The all-loop perturbative derivation of the NSVZ $\beta$-function and the NSVZ scheme in the non-Abelian case by summing singular contributions

The perturbative all-loop derivation of the NSVZ $\beta$-function for ${\cal N}=1$ supersymmetric gauge theories regularized by higher covariant derivatives is finalized by calculating the sum of singularities produced by quantum superfields. These singularities originate from integrals of double total derivatives and determine all contributions to the $\beta$-function starting from the two-loop approximation. Their sum is expressed in terms of the anomalous dimensions of the quantum gauge superfield, of the Faddeev--Popov ghosts, and of the matter superfields. This allows obtaining the NSVZ equation in the form of a relation between the $\beta$-function and these anomalous dimensions for the renormalization group functions defined in terms of the bare couplings. It holds for an arbitrary renormalization prescription supplementing the higher covariant derivative regularization. For the renormalization group functions defined in terms of the renormalized couplings we prove that in all loops one of the NSVZ schemes is given by the HD+MSL prescription.


Introduction
Ultraviolet divergences in supersymmetric theories are resticted by some non-renormalization theorems. In particular, in N = 1 supersymmetric theories the superpotential does not receive divergent quantum corrections [1]. However, even in the case of N = 1 supersymmetry there are some other similar statements. For example, the triple gauge-ghost vertices (which have two external lines of the Faddeev-Popov ghosts and one external line of the quantum gauge superfield) are also finite in all orders [2]. Moreover, according to Refs. [3,4,5,6], the βfunction of N = 1 supersymmetric gauge theories is related to the anomalous dimension of the matter superfields (γ φ ) j i by an equation which is usually called "the exact NSVZ β-function". For the theory with a simple gauge group G and chiral matter superfields in its representation R this equation is written as Note that we do not so far specify the definitions of the renormalization group functions (RGFs) and use the notations r ≡ dim G, where T A are the generators of the group G in the representation R such that [T A , T B ] = if ABC T C . Also we assume that the generators of the fundamental representation t A are normalized by the condition tr(t A t B ) = δ AB /2. Eq. (1) implies that the all-order β-function of the N = 1 supersymmetric Yang-Mills (SYM) theory without matter superfields is given by the geometric series. Moreover, if N = 2 supersymmetric gauge theories are considered as a special case of the N = 1 supersymmetric theories, then Eq. (1) leads to the finiteness beyond the one-loop approximation provided the quantization is made in N = 2 supersymmetric way [7,8]. This implies that the NSVZ βfunction is closely related to the N = 2 non-renormalization theorem derived in [9,10,11]. However, for its derivation N = 2 supersymmetry should be manifest at all steps of calculating quantum corrections. This can be achieved with the help of the harmonic superspace [12,13] and the invariant regularization [14]. The finiteness of N = 4 SYM theory [9,10,15,16] follows from the N = 2 non-renormalization theorem and, therefore, from the NSVZ β-function.
Equations analogous to NSVZ are also known for the Adler D-function in N = 1 SQCD [17,18] and for the renormalization of the gaugino mass in gauge theories with softly broken supersymmetry [19,20,21]. Recently an NSVZ-like equation was constructed for the renormalization of the Fayet-Iliopoulos term in D = 2, N = (2, 0) supersymmetric theories [22].
It is important that the NSVZ and NSVZ-like equations are valid only for certain renormalization prescriptions. In particular, for theories regularized by dimensional reduction [23] supplemented by modified minimal subtraction [24] (this scheme is usually called DR) Eq. (1) is not valid starting from the order O(α 2 , αλ 2 , λ 4 ), which corresponds to the three-loop β-function and the two-loop anomalous dimension [25,26,27,28,29]. 1 However, the detailed analysis made in these papers demonstrated that by a specially tuned finite renormalization of the gauge coupling constant one can restore the NSVZ equation (1), at least, in the three-and four-loop approximations. It should be noted that the possibility of this tuning is non-trivial, because of the presence of various group invariants (like C 2 , C(R) i j , etc.). If one considers only finite renormalizations polynomial in these invariants, then the NSVZ equation leads to some schemeindependent consequences [31,32]. This implies that, although the calculations made in the DR-scheme do not produce the NSVZ equation, they nevertheless confirm that it is valid in a certain (NSVZ) subtraction scheme. Using the general equations describing how the NSVZ equation changes under finite renormalizations [32,33], it is possible to construct an infinite set of the NSVZ schemes [34].
For a long time it was unknown, how to construct an all-order renormalization prescription giving the NSVZ scheme. However, recently it was understood that the solution can be found in the case of using the higher covariant derivative method [35,36] for regularizing supersymmetric theories. Unlike dimensional reduction [37], this regularization is mathematically consistent and can be formulated in a manifestly supersymmetric way in N = 1 superspace [38,39]. Although the calculations in higher derivative theories are rather complicated, some of them have been done in the last decades. For instance, a number of calculations of the one-loop effective potential for N = 1 higher derivative supersymmetric theories were made in [40,41,42,43,44,45,46]. In theories regularized by higher derivatives quantum corrections are obtained in a similar way. Some higher order calculations made with this regularization (see, e.g., [47,48,49,50]) demonstrated that the NSVZ equation is valid for RGFs defined in terms of the bare couplings. (In the Abelian case this has been proved in all orders [51,52]. Similar proofs of the NSVZ-like equations have also been constructed for the Adler D-function in N = 1 SQCD [17,18] and for the renormalization of the photino mass in N = 1 SQED with softly broken supersymmetry [53].) RGFs defined in terms of the bare couplings are scheme-independent for a fixed regularization [54], but depend on a regularization, so that Eq. (1) is valid for them for an arbitrary renormalization prescription supplementing the higher derivative regularization. The above mentioned calculations confirmed this in such an approximation where the dependence on a regularization is essential. Note that, according to [30,55], RGFs defined in terms of the bare couplings do not satisfy the NSVZ and NSVZ-like equations in the case of using dimensional reduction, again, for an arbitrary renormalization prescription supplementing it.
The important statement which allows constructing the NSVZ scheme is that RGFs defined in terms of the bare couplings and RGFs standardly defined in terms of the renormalized couplings up to the renaming of arguments coincide in the HD+MSL scheme [54], when the divergences of the theory regularized by Higher covariant Derivatives are removed by Minimal Subtractions of Logarithms [56,57]. This implies that the renormalization constants are constructed in such a way that they include only powers of ln Λ/µ, where Λ is a regularization parameter, analogous to the ultraviolet cut-off, and µ is a renormalization point. All finite constants in this case, by definition, are set to 0. The coincidence of various definitions of RGFs implies that the HD+MSL scheme appears to be one of the NSVZ schemes. In the Abelian case this has been proved in all loops [54]. Note that another all-loop NSVZ scheme in the Abelian case is the on-shell scheme [58]. For the renormalization of the photino mass in softly broken N = 1 SQED and for the Adler D-function in N = 1 SQCD the NSVZ-like schemes are also given by the HD+MSL prescription, see Refs. [59] and [60], respectively. In the non-Abelian case there are strong indications [2] that HD+MSL = NSVZ. Also there are some nontrivial explicit multiloop calculations confirming this fact [48,49,50], but, nevertheless, this statement has not yet been proved in all orders. This proof (started in Refs. [2,61,62]) will be finalized in this paper.
The main observation used for the derivation of the NSVZ and NSVZ-like equations for RGFs defined in terms of the bare couplings in theories regularized by higher derivatives is that the loop integrals giving the β-function are integrals of double total derivatives with respect to loop momenta. 2 Certainly, at least one of the loop integrals can be calculated analytically using equations like where f (Q 2 ) is a nonsingular function of the Euclidean momentum Q µ which rapidly tends to 0 at infinity. Note that the result does not vanish due to the nontrivial surface integration over the small sphere S 3 ε surrounding the point Q µ = 0. If we consider an L-loop contribution to the β-function, then Eq. (3) allows calculating one of the loop integrals, so that the resulting expression will contain only (L − 1) loop integrations. Therefore, it is possible to suggest that the result is a certain (L − 1)-loop quantum correction. According to [51,52], in the Abelian case it is really proportional to the (L − 1)-loop contribution to anomalous dimension of the matter superfields, so that the Abelian NSVZ equation [66,67] β(α 0 ) where N f is a number of flavors, is naturally produced for RGFs defined in terms of the bare couplings.
In the non-Abelian case the situation is more complicated, because Eq. (1) relates the βfunction to the anomalous dimension of the matter superfields in all previous orders due to the gauge coupling constant dependent denominator in the right hand side. However, according to Ref. [2], using the all-loop finiteness of triple gauge-ghost vertices the NSVZ equation (1) can be presented in an equivalent form which does not contain the coupling constant dependent denominator in the right hand side similarly to the Abelian case. (Note that in Eq. (5) we again do not specify the definitions of RGFs and omit some other possible arguments.) Similarly to the Abelian case, Eq. (5) relates an L-loop contribution only to (L − 1)-loop contributions to the anomalous dimensions of the quantum gauge superfield, of the Faddeev-Popov ghosts, and of the matter superfields, denoted by γ V , γ c , and (γ φ ) j i , respectively. This fact has a simple graphical interpretation analogous to the N = 1 SQED case considered in [64,68]. Namely, given a supergraph without external lines, the corresponding superdiagrams contributing to the β-function are constructed by attaching two external lines of the background gauge superfield in all possible ways, while the superdiagrams contributing to the anomalous dimensions are obtained by all possible cuts of internal lines. This qualitative picture is related to the structure of loop integrals giving the β-function defined in terms of the bare couplings. According to [61] they are integrals of double total derivatives in all orders in agreement with numerous calculations made with the higher covariant derivative regularization, see, e.g., [48,49,69,70,71]. Each cut of a certain propagator corresponds to taking an integral of a double total derivative analogous to (3). The sums of singularities generated by the cuts of various propagators produce the corresponding anomalous dimensions in Eq. (5) even at the level of loop integrals. 3 In the lowest orders this was explicitly demonstrated in [48,49,50,56,71]. In all loops the sums of singularities obtained by cutting the matter and Faddeev-Popov ghost propagators have been found in [62] using the method based on the Schwinger-Dyson effective superdiagrams proposed in [74]. These sums coincide with the terms containing (γ φ ) j i and γ c defined in terms of the bare couplings, respectively. Therefore, to complete the derivation of the NSVZ equation, one should find a sum of all singularities produced by cuts of the gauge propagators and demonstrate that it gives the term containing γ V in Eq. (5). 4 This is a purpose of the present paper.
The paper is organized as follows. Section 2 describes the quantization of N = 1 supersymmetric gauge theories regularized by higher covariant derivatives. In section 3 we rewrite the NSVZ β-function in the form of a relation between some two-point Green functions, which will be proved below. The proof is based on the method for constructing the integrals of double total derivatives giving the function β/α 2 0 which was proposed in [61], see also [50]. This method is described in section 4.1. It is slightly modified in section 4.2. Using this modification the sums of singularities produced by the matter superfields, by the Faddeev-Popov ghosts, and by the quantum gauge superfield are calculated exactly in all loops in section 5. In particular, in subsection 5.5 we find the sum of singularities produced by the quantum gauge superfield propagators, which is needed for finalizing the proof of the NSVZ β-function. Combining the results we derive Eq. (5) for RGFs defined in terms of the bare couplings. In the next section 6 we demonstrate that the HD+MSL prescription really gives one of the NSVZ schemes in all orders for RGFs defined in terms of the renormalized couplings. The last section 7 is devoted to the explicit perturbative verification of the results in the lowest orders of the perturbation theory. Doing the corresponding calculations we pay especial attention to checking the modification of the method for constructing integrals giving the function β/α 2 0 proposed in section 4.2.
2 N = 1 renormalizable supersymmetric gauge theories regularized by higher covariant derivatives We will consider a general renormalizable N = 1 SYM theory with a simple gauge group G and chiral matter superfields φ i in a representation R. In the massless limit in terms of N = 1 superfields [76,77,78] this theory is described by the action where for the superspace integration measures we use the brief notations In Eq. (6) V is the Hermitian background gauge superfield and V is the quantum gauge superfield, which satisfies the constraint V + = e −2V V e 2V . Note that in the first term of Eq. (6) the quantum and background gauge superfields are expanded in the generators of the fundamental representation as where T A are the generators of the gauge group in the representation R.
The gauge and Yukawa couplings are denoted by e 0 and λ ijk 0 , respectively, where the subscript 0 points the bare values. The function F(V ) = e 0 F(V ) A t A is needed, because the quantum gauge superfield is renormalized in the non-linear way [79,80,81]. In the lowest approximation this function was found in Refs. [82,83] and is written as where y 0 is the first bare constant in an infinite set of parameters describing the nonlinearity, and G ABCD ≡ (f AKL f BLM f CM N f DN K + permutations of B, C, and D)/6. The necessity of introducing the function F(V ) was also confirmed by the calculation of Ref. [84], where it was demonstrated that the renormalization group equations are satisfied only if the renormalization of the parameter y is taken into account. The gauge superfield strength is described by the chiral superfield which is a right spinor with respect to the Lorentz group. If the Yukawa couplings satisfy the equation the theory (6) is invariant under the background gauge transformations parameterized by the chiral superfield A which takes values in the Lie algebra of the gauge group G. The action (6) is also invariant under the quantum gauge transformations, but this invariance is broken by the gauge fixing procedure. The regularization is implemented by inserting into the classical action the higher derivative regulator functions R and F . They should rapidly increase at infinity and satisfy the condition R(0) = F (0) = 1. Then the action of the regularized theory can be written as where for f ( and the explicit expressions for the covariant derivatives have the form The modification of the action S → S reg allows regularizing all divergences beyond the one-loop approximation [85]. The generating functional for the regularized theory should also include a gauge fixing term, ghost actions, sources, and Pauli-Villars determinants needed for removing the remaining oneloop divergences [86], The source term is given by the expression where J A and j i are the real and chiral superfields, respectively. The anticommuting chiral superfields j A c , andj A c are the sources for the anticommuting chiral Faddeev-Popov ghost and antighost superfields denoted by c A andc A , respectively.
We will use the gauge fixing term invariant under the background gauge transformations (11) due to the presence of the background covariant derivatives The corresponding actions for the Faddeev-Popov and Nielsen-Kallosh ghosts (S FP and S NK , respectively) can be found in Refs. [61,62]. In this paper we will not use the explicit expressions for them. The bare gauge parameter ξ 0 and the parameters present in the function F(V ) (i.e., y 0 , etc.) can conveniently be included into a single infinite set Y 0 ≡ (ξ 0 , y 0 , . . .). Following Refs. [75,87], for regularizing the one-loop divergences we use two sets of the Pauli-Villars superfields. The first one consists of three commuting chiral superfields ϕ a in the adjoint representation of the gauge group with the mass M ϕ = a ϕ Λ. They cancel the divergences coming from the gauge and ghost loops. The corresponding determinant present in Eq. (15) is denoted by Det(P V, M ϕ ). The second set giving the determinant Det(P V, M ) consists of the commuting chiral superfields Φ i in a certain representation R PV that admits the gauge invariant mass term such that M ij M jk = M 2 δ i k with M = aΛ. 5 These superfields cancel the one-loop divergences produced by the matter superfields φ i if in Eq. (15) the degree of the corresponding determinant is c = T (R)/T (R PV ). Again we will not use explicit expressions of the Pauli-Villars determinants which can also be found in Refs. [61,62]. It should be only mentioned that the constants a ϕ and a must be independent of couplings. Starting from the generating functional for the connected Green functions W ≡ −i ln Z, one can construct the effective action where it is necessary to express the sources J A , j i , . . . in terms of quantum superfields V A , φ i , . . . using the equations 3 The NSVZ equation as a relation between two-point Green functions The NSVZ equation (5) for RGFs defined in terms of the bare couplings in the case of using the higher covariant derivative regularization can be written as a certain equation relating twopoint Green functions of the theory. The terms in the effective action corresponding to these Green functions can be presented in the form where ∂ 2 Π 1/2 ≡ −D aD2 D a /8 is the supersymmetric transversal projection operator. Writing Eq. (21) we took into account that the two-point Green function of the background gauge superfield is transversal due to the manifest background gauge invariance of the effective action. Similarly, in Eq. (22) we used the fact that quantum corrections to the two-point Green function of the quantum gauge superfield are also transversal due to the Slavnov-Taylor identity [72,73].
In our notation the renormalization constants are introduced with the help of the equations 6 where the subscript R stands for the renormalized superfields. The renormalization constants are constructed from the requirement that the functions d −1 , Z 2 V G V , (Z φ G φ ) i j , and Z c G c expressed in terms of the renormalized couplings α, λ ijk , and Y should be finite in the limit Λ → ∞. Note that due to the non-renormalization of the superpotential [1] the renormalized Yukawa couplings are given by We will always assume that no finite constants corresponding to finite renormalizations appear in this equation, so that Eq. (26) partially fixes the subtraction scheme. Similarly, due to the all-loop finiteness of the triple gauge-ghost vertices [2] 7 it is possible to choose a subtraction scheme in which the renormalization constants satisfy the condition Note that this equation is compatible with minimal subtractions of logarithms, because in the HD+MSL scheme all renormalization constants Z α , Z c , and Z V are equal to 1 for µ = Λ. The equation (27) allows to present the NSVZ equation in an equivalent form relating the β-function to the anomalous dimensions of the quantum gauge superfield, the Faddeev-Popov ghosts, and the matter superfields, where we take into account that RGFs (at least, γ V and γ c ) may in general depend on Y 0 . RGFs defined in terms of the bare couplings entering this equation can be related to the corresponding two-point Green functions by the equations That is why Eq. (28) can also be rewritten as an equation relating these Green functions which was first proposed in Ref. [2]. In this paper we will prove it in all orders of the perturbation theory.

Integrals of total derivatives 4.1 How to construct and calculate integrals of total derivatives
The key observation needed for deriving the NSVZ β-function is the factorization of loop integrals which give the function β(α 0 , λ 0 , Y 0 )/α 2 0 into integrals of double total derivatives in the case of using the higher covariant derivative regularization. The all-loop proof of this fact was done in Ref. [61]. The ideas used in this proof allowed constructing a method for obtaining these integrals in each order of the perturbation theory. For this purpose one should calculate only specially modified vacuum supergraphs 8 instead of a much larger number of superdiagrams with two external legs of the background gauge superfield. Some higher-order calculations made with the help of this method in Refs. [50,71,90] demonstrated that it works correctly and reproduces all known results. Moreover, the structure of quantum corrections which were first obtained by this method is in excellent agreement with some general theorems. Say, the β-function in the considered theories appeared to be gauge independent and satisfies the NSVZ equation in the HD+MSL scheme. Here using this method we will find an exact all-order expression for the function where is the one-loop β-function defined in terms of the bare couplings. (For the higher covariant derivative regularization considered in this paper it was calculated in Ref. [75].) Now, let us formulate the algorithm for constructing contributions to the function (34) following Refs. [50,61].
1. As a starting point we consider a vacuum supergraph with L loops and construct an expression for it using the superspace Feynman rules.
2. Next, one needs to find a point with the integration over d 4 θ (or convert the integration over d 2 θ into the integration over d 4 θ) and insert the expression θ 4 (v B ) 2 to the corresponding point of this supergraph. Here v B denotes a function slowly decreasing at a very large scale R → ∞. (Note that without this insertion any vacuum supergraph vanishes due to the integration over the anticommuting variables θ.) 3. We calculate the supergraph modified by the above insertion and omit terms suppressed by powers of 1/(RΛ). As a result we obtain an expression proportional to 4. At the next step it is necessary to choose L propagators with independent momenta denoted by Q µ i , where i = 1, . . . , L. (In our notation the capital letters denote Euclidean momenta which appear after the Wick rotation.) Let the gauge group indices corresponding to the beginnings and endings of these propagators be a i and b i , respectively. Then the product of the marked propagators contains the factor L i=1 δ b i a i . 5. The product L i=1 δ b i a i coming from the marked propagators in the integrand of the (Euclidean) loop integral should formally be replaced by the operator 6. At the last step, it is necessary to apply the operator to the resulting expression, where the derivative with respect to ln Λ should be calculated at fixed vales of the renormalized couplings prior to the integration over loop momenta.
After the above described procedure we obtain a contribution to the function (34) corresponding to the considered supergraph. It is produced by the sum of all two-point superdiagrams which are constructed from this supergraph by attaching two external lines of the background gauge superfield V in all possible ways. By construction, the result is given by a certain integral of double total derivatives with respect to the loop momenta.
The integrals of total derivatives do not vanish due to the singularities of the integrands, which appear when two momentum derivatives act on an inverse squared momentum, We will also need a modification of this identity which is obtained when the derivatives are taken with respect to the different momenta, where a 1 and a 2 are some constants. Note that in calculating the integrals of double total derivatives we should take these singular contributions with the opposite sign. Really, if f (Q 2 ) is a non-singular function which rapidly tends to 0 at infinity, then Note that terms in which double total derivatives act on Q −4 are not well-defined and cannot appear in the final expression for the function (34), although they can be present in expressions for separate supergraphs. This statement has been confirmed by some explicit two-and threeloop calculations in Refs. [49,71].

Graphs and total derivatives
Let consider a vacuum supergraph with L loops, V vertices, and P internal lines and set directions of all internal lines in an arbitrary way. These directions will be pointed by arrows. Also we denote momenta of all internal lines by certain letters. In our conventions an incoming momentum has the sign "minus", while an outcoming one has the sign "plus".
Let us construct a (V −1)×P matrix M corresponding to the supergraph under consideration. For this purpose we numerate the vertices in an arbitrary order and write the energy-momentum conservation laws in all vertices, except for the last one. (Evidently, the last equation is a linear combination of the others.) The resulting system of equations can be written in the form The matrix with the elements M i,J , where i = 1, . . . , V − 1 and J = 1, . . . , P , is the required matrix matched to the considered vacuum supergraph. From the above equations, which can briefly be written as we conclude that only P − V + 1 momenta are independent. According to the topological identity L = P − V + 1, this implies that (as well known) there are L independent momenta Q µ,α , α = 1, . . . , L in the considered supergraph, and the others can be expressed in terms of them, As a simple example we can consider a supergraph presented in Fig. 1 on the left. In this case the equation (43), the matrix M , and the equation (44) read as (45) For a more complicated three-loop graph presented in Fig. 1 on the right as another example, the matrix M takes the form Evidently, in this case there are three independent momenta. For example, it is possible to choose Q µ,1 , Q µ,2 , and Q µ,3 as independent variables. The corresponding propagators are denoted in Fig. 1 by the bold lines. Then the equation (44) takes the form An important observation made in Ref. [61] is that the equations which reflect the gauge invariance of vertices are very similar to Eq. (43). Really, let us consider a vertex with n outcoming lines corresponding to the term whereV is an operator acting on the product of various superfields of the theory denoted by ϕ I . Then the energy-momentum conservation in the considered vertex is expressed by the equation which is very similar to the equation which follows from the gauge invariance of theory where (T A ) I J are generators of the gauge group in a relevant representation. For the incoming lines the signs of the corresponding momenta in equations like (49) should be changed. Similarly, in equations like (50) one should replace the generators by the transposed generators of the conjugated representation taken with an opposite sign, although this replacement does not change the equation, because Let us consider a vacuum supergraph and replace a δ-symbol δ b J a J coming from the propagator with a number J by a certain matrix A a J b J . We will denote the resulting modified supergraph by [A] J . If we replace the i-th vertex operatorV I 1 I 2 ...In in a certain vacuum supergraph by the left hand side of Eq. (50), then, using this notation, we obtain the equation The system of these equations is analogous to Eq. (43) and contains the same matrix M i,J . Evidently, the solution of these equations can be written in a form similar to Eq. (44), where N I,α are exactly the same coefficients as in Eq. (44). According to the algorithm described in section 4.1 we need to insert θ 4 (v B ) 2 into the vacuum supergraph and make the replacement described in the item 5. We will denote supergraphs obtained in this way as L α,β=1 This means that θ 4 (v B ) 2 is inserted into an arbitrary point of the supergraph, two δ-symbols are replaced by the generators, and the double total derivatives are applied to the integrand of the loop integral.
Then the sum of all perturbative contributions to the function (34) can be presented in the form (55) By construction, this expression is an integral of double total derivatives. However, it does not vanish because the double total derivatives produce singular contributions acting on 1/Q 2 I . It is important to note that expressions for separate supergraphs (presented in the form of scalar integrals) can contain not only 1/Q 2 I , but also 1/(Q 2 I ) n ≡ 1/Q 2n I with n ≥ 2 if Q I is a momentum of a quantum gauge superfield propagator. Acting on 1/Q 2n I with n ≥ 2 the double total derivatives produce an expression which is not well-defined due to infrared divergences. However, we know that the left hand side of Eq. (55) is well-defined. Therefore, all bad terms coming from the 1/Q 2n I singularities should cancel each other. The calculations in the lowest orders [49,71] exactly confirm this statement.
Because all bad terms should cancel each other, it is necessary to consider only the 1/Q 2 I singularities. Then in the expression for a supergraph we need to consider only the product of all different inverse squared momenta multiplied by a nonsingular function f , where Q µ,I = Q µ,J for I = J. The prime indicates that the product includes only different momenta, a number of them being denoted by P ′ (which is evidently less or equal to P ). This is essential for supergraphs which contain coinciding momenta. The structure of such supergraphs is illustrated in Fig. 2. In this figure the gray circles denote 1PI subdiagrams, which are connected by propagators with coinciding momenta. By construction, the expression (56) should include only one such momentum.
k µ Figure 2: The structure of vacuum supergraphs with coinciding momenta (which are equal to k µ ). The gray circles denote 1PI subdiagrams.
δ-singularities appear when double total derivatives act on various factors in the product (56). Note that, according to Eq. (40), such singularities can appear even if derivatives with respect to different momenta act on the same inverse squared momentum, where we also took Eq. (44) into account. This implies that acting on the product (56) the double total derivatives give the singular contribution Taking into account that the bad terms containing 1/Q 2n I with n ≥ 2 should cancel each other, we can rewrite the expression (55) as a sum of supergraphs in which they are omitted. According to Eq. (58), in the remaining terms we replace one of 1/Q 2 I by 4π 2 δ 4 (Q I ) and sum up all expressions obtained after these substitutions. (The sign is "plus", because the singular contribution of the form (39) should be taken with the minus sign, see Eq. (41).) The result can be presented as Next, using Eq. (53) we calculate the sums over the indices α and β, Note that for the gauge and ghost propagators ( Basing on Eq. (60) it is possible to modify the step 5 of the algorithm described in section 4.1. Note that in what follows we formulate this step before the Wick rotation.
5. We construct a set of supergraphs in which one of the propagators is marked. If the original supergraph contains some propagators with the same momentum, then one can mark only one propagator with this momentum. After calculating the D-algebra, one should omit all terms proportional to 1/q 2n I with n ≥ 2 in scalar integrals and for any marked propagator with the momentum q µ in the remaining terms make the replacement for gauge and ghost propagators; After these replacements it is necessary to sum up the expressions for all supergraphs in the constructed set. Note that this prescription is the same for all vacuum supergraphs, unlike the one described in section 4.1. Really, the form of the operator (37) depends on the topology of the supergraph, while the replacement (61) does not. That is why it is possible, first, to find a sum of all vacuum supergraphs and, only after this, make these replacements. Below we will demonstrate that the sum of vacuum supergraphs really does not contain any bad terms which include 1/Q 2n I with n ≥ 2. This confirms the above argumentation, although the explicit mechanism of cancelling the bad contributions to the function (34) should be analyzed separately.
In the next section using the modification of the algorithm described above we will find the sum of all contributions to the function (34) which come from δ-singularities of the form (40). Various perturbative verifications of the above described modification will be made in section 7.
5 The all-loop sum of singular contributions 5

.1 The starting point and main idea
In this section we will calculate the sum of all singular contributions. They appear when the double total derivatives act on matter, ghost and gauge propagators and effectively cut the corresponding internal line. In [62] the sums of singularities produced by the cuts of matter and ghost propagators have been calculated in all orders. However, the method which is used in this paper is different. That is why it is expedient to recalculate the sums of matter and ghost singularities and verify that both approaches give the same result. Moreover, in this section we will find the all-loop expression for a sum of singularities produced by the cuts of quantum gauge superfield propagators. The method of summation generalizes the one proposed in Ref. [51] for the Abelian case.
First, we introduce the auxiliary action depending on a real constant g, where S total ≡ S reg + S gf + S FP + S NK and For g = 1 the action S g evidently coincides with S total . Moreover, all vertices generated by S g coincide with the ones produced by S total . However, the parts of S g quadratic in the quantum gauge superfield, the Faddeev-Popov ghosts, and the matter superfields differ from the corresponding quadratic parts of S total by the factor 1/g. This implies that, in comparison with the original theory, all propagators of the above mentioned quantum superfields obtained from the generating functional with c = T (R)/T (R PV ), acquire the factor g. Note that the propagators of the massive Pauli-Villars superfields and of the Nielsen-Kallosh ghosts remain the same. We need not modify them, because the Pauli-Villars propagators do not produce singularities, and the Nielsen-Kallosh ghosts are essential only in the one-loop approximation which is considered separately. Certainly, the functional Z g is treated formally, and we do not care about regularizing divergences in the corresponding supergraphs. Now, let us find the derivative of the functional Γ g (defined as a Legendre transform of W g = −i ln Z g ) with respect to the parameter g at vanishing superfields and g = 1. It is convenient to present the result in the form The angular brackets entering this equation are defined in the standard way where the sources should be expressed in terms of fields using Eq. (20). Evidently, the expression (65) can be rewritten in terms of the inverse two-point Green functions of the quantum superfields, (Certainly, in the last expression we also assume that g = 1 and the (super)fields are set to 0.) The original effective action Γ at vanishing fields is given by the sum of vacuum supergraphs. So far we do not calculate these supergraphs and work only with the formal expressions constructed with the help of supersymmetric Feynman rules. The derivative (65) is contributed by the same vacuum supergraphs, but each of them is multiplied by the number of propagators. Really, as we have already mentioned, each propagator of a quantum superfield is proportional to the parameter g, while the vertices do not contain g. Evidently, the first and second terms in Eq. (67) are contributed by vacuum supergraphs multiplied by the number of gauge and matter propagators, respectively. The last two terms give a sum of vacuum supergraphs multiplied by the number of ghost propagators.
Note that if a supergraph has some coinciding momenta as in Fig. 2, all propagators with the coinciding momenta should be taken into account when calculating the number of propagators. The contribution of these graphs to ∂Γ g /∂g is given by a relevant term in Eq. (67) and is proportional to a certain inverse Green function. If a usual two-point Green function is proportional to the function depending on the momentum, then the inverse Green function will include Evidently, a term containing (∆G) n corresponds to supergraphs of the structure presented in Fig. 2 with n coinciding momenta k µ . The gray circles in this supergraph certainly give various ∆G.
As we discussed in the previous section, for constructing the β-function one should cut only one propagator with the momentum k µ in a supergraph which have the structure presented in Fig. 2. (This cut is made by δ 4 (k) which appears after the replacement (61).) From the other side, the corresponding contribution to the expression (65) obtained by counting the propagators with the momentum k µ is proportional to n. Therefore, for calculating the contribution to the β-function it is necessary to replace G −1 by (The first term in Eq. (69) corresponds to the one-loop approximation, which is considered separately.) After this we proceed according to the algorithm described in section 4. Details of this calculation will be considered below, separately for the cuts of gauge, ghost, and matter propagators.

The all-loop sum of matter singularities
Let us start with calculating a contribution to the function (34) coming from the cuts of matter propagators. According to Eq. (67), the corresponding contribution to ∂Γ g /∂g is given by the expression which encodes the sum of vacuum supergraphs with a marked matter propagator. (By definition, the marking of a propagator does not change the expression for a supergraph, but supergraphs in which marked propagators are different are considered as different.) The two-point Green function of the matter superfields obtained by differentiating Eq. (23) with respect to φ and φ * reads as where The function (G φ ) i j present in Eq. (72) is normalized in such a way that in the tree approximation it is equal to δ j i F (∂ 2 /Λ 2 ). (In the limit Λ → ∞ this expression gives δ j i .) Therefore, it can be presented as where the sum of quantum corrections is denoted by (∆G φ ) i j . The function inverse to (72) entering Eq. (71) is given by the expression This can easily be verified using the first of the identities To obtain an expression which encodes the sum of vacuum supergraphs in which only one of propagators with coinciding momenta is marked, one should multiply terms with n insertions of ∆G φ by 1/n. (As we have already mentioned, the first term, which corresponds to the one-loop approximation, should be omitted.) After this we obtain the function Therefore, the first step for calculating the matter contribution to the function (34) is to make the replacement in Eq. (71). We see that no bad terms proportional to 1/q 2k with k ≥ 2 appear. Next, following the algorithm described in section 4, we should insert θ 4 (v B ) 2 to an arbitrary point of the supergraph. Evidently, in this case it is expedient to insert this expression to the point x.
Moreover, we need to make a replacement where q µ is the Minkowski momentum of the matter line which is cut. Note that (as we discussed in section 4.2) the operator ∂ 2 /∂q µ ∂q µ should act only on 1/q 2 . (After the Wick rotation it gives −∂ 2 /∂Q 2 µ , where Q µ is the corresponding Euclidean momentum.) The matter line is cut in the point x, so that q µ is the momentum of the propagator coming out of this point. Finally the result should be multiplied by −2π/(rV 4 ) and differentiated with respect to ln Λ at fixed values of renormalized couplings. This implies that the matter contribution to the function (34) can be written in the form where the function (ln G φ ) i j depends on the momentum q µ and The momentum integral is calculated in the Euclidean space after the Wick rotation, where we took into account that F (0) = 1 and wrote the result in terms of the anomalous dimension (γ φ ) j i using Eq. (32). Therefore, the expression (80) takes the form and calculate the integral over d 8 y with the help of the δ-function, Taking into account that we obtain the final expression for the sum of matter singularities It exactly coincides with the term containing the anomalous dimension of the matter superfields in the NSVZ equation written in the form (28) and agrees with another all-loop calculation made in [62] by a different method. This can be considered as a correctness test of the method proposed in this paper.

The all-loop sum of ghost singularities
A contribution of ghost singularities to the function (34) is calculated similarly to the matter contribution, but it is necessary to take into account that the ghost superfields are anticommuting. As a starting point we consider a part of the expression (67) containing the ghost Green functions, The two-point Green functions of the Faddeev-Popov ghosts are obtained by differentiating Eq. , In the tree approximation the function G c is equal to 1, so that it is convenient to present it in the form G c = 1 + ∆G c , where ∆G c is a sum of quantum corrections coming from 1PI supergraphs with two external ghost legs. By definition, the functions inverse to (89) satisfy the equations and are explicitly written as Repeating the transformations similar to the ones made in the previous section for the matter singularities, we conclude that for constructing a ghost contribution to the function (34), first, it is necessary to perform a substitution in the expression (88). Next, according to the algorithm described in section 4, θ 4 (v B ) 2 should be inserted to the point x. After this, it is necessary to make the replacement 9 and apply the operator (38) to the resulting expression. Note that the differentiation with respect to ln Λ should be made at fixed values of the renormalized couplings. As a result of this algorithm, we obtain the ghost contribution to the function (34) (Note that again no bad terms appear in this expression.) Similarly to the case of matter singularities, the momentum integral is calculated after the Wick rotation in the Euclidean space, With the help of this equation the considered contribution to the function (34) can be brought to the form Then we use the identity (84) and calculate the integral over d 8 y with the help of resulting δ 8 xy , According to Eq. (86), the remaining integral present in this equation is equal to 4V 4 , so that the sum of singular contributions produced by cuts of the Faddeev-Popov ghost propagators takes the form This result agrees with the one found in Ref. [62] by a different method and coincides with the term containing the anomalous dimension γ c in Eq. (28). Note that the sign of the expression (98) is different from the sign of the matter contribution (87) due to the anticommutation of the ghost superfields. The factor 2 in the ghost contribution appears because there are two sets of the chiral ghost superfields (the ghost c and the antighostc) in the adjoint representation of the gauge group.

Effective propagator of the quantum gauge superfield
Taking into account that (for g = 1) due to the Slavnov-Taylor identity quantum corrections to the two-point Green function of the quantum gauge superfield are transversal, it is possible to write the corresponding part of the effective action in the form (22). Substituting the explicit expression for the gauge fixing action, we equivalently present this equation as where the regulator function K depends on ∂ 2 /Λ 2 , and the argument q 2 of function G V should be replaced by −∂ 2 . Differentiating this expression with respect to the quantum gauge superfield, we obtain the corresponding two-point Green function By definition, the inverse function satisfies the equation Using the identity (see, e.g., [77]) which can be proved by (anti)commuting the covariant derivatives, it is possible to demonstrate that the explicit expression for the inverse two-point Green function of the quantum gauge superfield is written as Thus, we see that the effective propagator is proportional to Q −4 , where Q µ is the corresponding Euclidean momentum.

The all-loop sum of singularities produced by the quantum gauge superfield
A part of the expression (67) contributed by supergraphs in which one gauge propagator is marked is written as where the exact propagator of the quantum gauge superfield is given by Eq. (103). All quantum corrections inside this exact propagator are encoded in the function G V , which is equal to R(∂ 2 /Λ 2 ) in the tree approximation. That is why it is convenient to present it in the form where ∆G V is a sum of relevant quantum corrections. Exactly as for the matter and ghost contributions, we rewrite the inverse two-point Green function in Eq. (104) in the form of a series using the identities and (76). To simplify the resulting expression, we introduce the notation This expression is proportional to the usual (tree) propagator of the quantum gauge superfield. Then the inverse Green function present in Eq. (104) can equivalently be rewritten as The factors ∂ 2 Π 1/2 G V in this expression correspond to the 1PI subdiagrams (denoted in Fig.  2 by the gray circles), which are evidently transversal. As we discussed in section 5.1, for constructing a contribution to the function (34), one should first divide a term with the n-th power of ∆G V by n. (Certainly, the first term corresponding to the one-loop approximation should be omitted.) Then the function (108) will be replaced by the expression Calculating the sum of this series, we obtain that, at the first step, it is necessary to make in Eq. (104) the formal substitution After this, with the help of Eq. (106) we obtain It should be noted that all possible bad terms proportional to 1/q 2k with k ≥ 2 vanish, and the above expression really contains only the 1/q 2 singularity. To construct a contribution to the function (34) coming from the gauge singularities, we need to modify the expression (111) in a special way. Namely, it is necessary to insert θ 4 (v B ) 2 to the point x, apply the operator C 2 ∂ 2 /∂q µ ∂q µ to 1/q 2 coming from the gauge propagator which is cut in the point x, multiply the result by −2π/(rV 4 ), and differentiate it with respect to ln Λ. The sum of the gauge singularities constructed according to this procedure is given by As earlier, we calculate the momentum integral in the Euclidean space after the Wick rotation taking into account that R(0) = 1, For calculating the remaining integral we use the identity Therefore, Substituting this expression into Eq. (114) we obtain the final result for the sum of singular contributions to the function (34) produced by the quantum gauge superfield, We see that it coincides with the corresponding term in Eq. (28) in exact agreement with the guess made in Ref. [2].

The β-function defined in terms of the bare couplings
The overall result for the function (34) is obtained by summing the contributions produced by cuts of matter, ghost, and gauge propagators, which are given by Eqs. (87), (98), and (117), respectively, Taking into account that the one-loop contribution to the β-function is given by Eq. (35), we see that the resulting expression for the β-function defined in terms of the bare couplings coincides with the NSVZ equation written in the form (28). Certainly, it is highly important that the theory is regularized by higher covariant derivatives. As we have already mentioned, for RGFs defined in terms of the bare couplings the NSVZ equation does not hold in the case of using the regularization by dimensional reduction [55] due to a different structure of loop integrals [65]. Note that, in fact, in the previous sections we have also proved the equation (33) relating the two-point Green function of the considered theory in the limit of the vanishing external momenta.
To obtain the NSVZ relation in the usual form, one needs to involve the non-renormalization theorem for the triple gauge-ghost vertices. Namely, following Ref. [2], we differentiate Eq. (27) (which is a consequence of this theorem) with respect to ln Λ at fixed values of renormalized couplings and obtain the equation Excluding the sum γ c + γ V from Eqs. (28) and (119) we obtain the NSVZ relation in the original form of the relation between the β-function and the anomalous dimension of the matter superfields, Thus, we have proved that it is valid in all loops for RGFs defined in terms of the bare couplings in the case of using the regularization by higher covariant derivatives. The subtraction scheme in which the NSVZ equations hold for RGFs defined in terms of the renormalized couplings will be constructed in the next section.

NSVZ scheme for RGFs defined in terms of the renormalized couplings
We have proved that the NSVZ relations (1) and (5) are satisfied by RGFs defined in terms of the bare couplings in the case of using the higher covariant derivative regularization described in section 2. Note that these RGFs are scheme-independent [54], so that the NSVZ relations for them are valid independently of a renormalization prescription. However, the standard RGFs are defined in terms of the renormalized couplings and depend on both a regularization and a subtraction scheme. It is known that the NSVZ relation for RGFs defined in terms of the renormalized couplings is satisfied only in certain subtraction schemes, which are called the NSVZ schemes. In particular, the DR-scheme is not the NSVZ scheme [25,26,27]. However, using the results described above it is possible to demonstrate that an NSVZ scheme can be obtained in all orders with the help of the HD+MSL prescription [2]. For completeness, here we briefly explain, how this statement can be obtained.
In terms of the renormalized couplings RGFs are defined by the equations The key observation made in [54] is that both definitions of RGFs give the same functions (of different arguments) if certain boundary conditions are imposed on the renormalization constants. In the non-Abelian case considered here these boundary conditions (in the point x 0 which is a fixed value of ln Λ/µ) take the form We also assume that the renormalization of the Yukawa couplings is related to the renormalization of the matter superfields by Eq. (26). Note that, due to the boundary conditions (122) and the nonrenormalization of the triple ghost-gauge vertices, Eq. (27) is satisfied automatically, because in this case the equation has the only solution Z Similarly, it is easy to see that the condition Z ξ Z 2 V = 1 is also satisfied automatically.
The conditions (122) define a class of the HD+MSL-like schemes. Here HD means that the higher covariant derivative method is used for regularizing the theory under consideration. This is always assumed in this paper. The HD+MSL scheme is obtained for x 0 = 0. In this case all renormalization constants include only powers of ln Λ/µ and all finite constants which define a subtraction scheme vanish. The schemes corresponding to x 0 = 0 are related to the HD+MSL scheme by the redefinition of the renormalization point µ = exp(−x 0 )µ HD+MSL exactly as in the Abelian case considered in [34].
Under the boundary conditions (122) RGFs (121) can be related to RGFs (29) -(32) by the equations where the arrows point that it is necessary to make a formal change of the argument, say, to write α instead of α 0 , etc. All these equations can be proved repeating the argumentation of [54] (in which similar equations were derived in the Abelian case). For example, taking into account that α = α 0 Z α , we obtain Here the derivative d/d ln µ acts on both ln Λ/µ inside α, λ, Y and the explicitly written ln Λ/µ. The partial derivative ∂/∂ ln µ, by contrast, acts only on the explicitly written ln Λ/µ and does not act on ln Λ/µ inside α, λ, and Y . Therefore, Next, we consider Eq. (125) in the point ln Λ/µ = x 0 . Due to the boundary conditions (122) ln Z α (α, λ, Y, ln Λ/µ → x 0 ) = 0 for all values of α, λ, and Y . This implies that However, according to Eq. (122), in the point ln Λ/µ = x 0 the values of the bare and renormalized couplings coincide, (For the Yukawa couplings it is also necessary to use Eq. (26), which relates the renormalization of the Yukawa couplings to the renormalization of the chiral matter superfields.) Taking into account that the left hand side of Eq. (125) is considered as a function of the renormalized couplings, the right hand side should also be expressed in terms of them using Eq. (129). This implies that we should formally replace the bare couplings by the renormalized ones, so that The other equations in (124) can be proved in a similar way. For example, with the help of the chain rule for the derivative with respect to ln µ we can present the anomalous dimension of the quantum gauge superfield in the form Again we consider this equation in the point ln Λ/µ = x 0 and take into account that due to Eq. (122) ln Z V (α, λ, Y, ln Λ/µ → x 0 ) = 0. Then, repeating the above argumentation, we obtain Earlier we have demonstrated that RGFs defined in terms of the bare couplings satisfy the NSVZ relations (28) and (120) for theories regularized by higher covariant derivatives independently of a renormalization prescription. (Let us recall that these RGFs are scheme-independent for a fixed regularization.) After the formal change of arguments α 0 → α, λ 0 → λ, Y 0 → Y these equalities certainly remain valid. Therefore, from Eq. (124) we conclude that in the case of using the higher covariant derivative regularization supplemented by the renormalization prescription (122) RGFs defined in terms of the renormalized couplings also satisfy the equations This implies that in the non-Abelian case the prescription (122) also provides the NSVZ scheme in all loops. The HD+MSL prescription is obtained by imposing the boundary conditions (122) with x 0 = 0. For other values of x 0 we obtain a family of schemes which differ from HD+MSL by redefinitions of the normalization point µ. Certainly, in these schemes the NSVZ relation is also valid in all loops. These statements agree with the explicit calculations (of some scheme dependent terms in RGFs) made in Refs. [48,49].

Verifications in the lowest orders
In this section we verify the general argumentation discussed above by explicit calculations in the lowest orders of the perturbation theory made with the help of the higher covariant derivative regularization. For this purpose we will use the results for various groups of supergraphs obtained earlier. Namely, using the method described in section 4.1 the two-loop β-function for a general N = 1 supersymmetric gauge theory with a simple gauge group has been calculated in Ref. [71]. Also using this method the parts of the three-loop β-function containing the Yukawa couplings and ghost loops have been found in Refs. [61] and [50], respectively. Note that before this the part of the three-loop β-function containing the Yukawa couplings has also been calculated with the help of the standard technique in Refs. [48,49]. For N = 1 SQED with N f flavors the complete three-loop β-function in a general ξ-gauge has been calculated by the algorithm of section 4.1 in Ref. [90]. Here all these calculations are used for checking the exact results described in the previous sections. Note that we will not verify that the equations (28), (120), and (133) really hold, because this has already been done in Refs. [48,49,50,71,90]. The main purpose of this section is to test the argumentation which was used for the all-loop derivation of the NSVZ equations at intermediate steps.

The two-loop approximation
First, we reanalyse the two-loop calculation made in Ref. [71]. The corresponding contribution to the β-function is generated by the vacuum supergraphs presented in Fig. 3. The standard technique requires calculating all superdiagrams contributing to the two-point Green function of the background gauge superfield. They are obtained from the ones presented in Fig.  3 by attaching two external V -legs in all possible ways. The method of Ref. [61] considerably simplifies the calculation, because it deals only with vacuum supergraphs. In this paper we have made some modifications, which will be verified here. This allows checking the general argumentation used for deriving the NSVZ equation and illustrating it by explicit calculations.
Let us start with the supergraph B1. The expression for it has been calculated in Ref. [61]. If we include the factor −2π/(rV 4 ) · d/d ln Λ, the result can be written as B1 B2 B3 B4 B5 B6 B7 Figure 3: Vacuum supergraphs generating the two-loop contribution to the β-function.
where F K ≡ F (K 2 /Λ 2 ). The integrand here contains three inverse squared momenta. Therefore, we should sum up three (equal) expressions obtained by replacing one of these inverse squared momenta 1/P 2 (together with the corresponding δ-symbol δ n m ) by 4π 2 C(R) m n δ 4 (P ). This gives the contribution to the function (34) of the form This result exactly agrees with the one obtained in Ref. [48] by the direct calculation of two-point superdiagrams with two external V -legs. The supergraphs B2 and B3 containing a matter loop produce two different contributions. If a matter loop corresponds to the superfields φ i and to the Pauli-Villars superfields Φ i , then (for an arbitrary value of ξ 0 ) the result for the (properly modified) vacuum supergraphs multiplied by the operator (38) is given by where R K ≡ R(K 2 /Λ 2 ). Unlike the separate supergraphs B2 and B3, it does not contain bad singularities proportional to the inverse momenta to the fourth power and terms which depend on ξ 0 . Next, we should find a sum of the expressions obtained from (136) either by replacing 1/K 2 (which comes from the gauge propagator) by 4π 2 C 2 δ 4 (K) or by replacing δ j i /Q 2 or δ j i /(Q + K) 2 (coming from the propagators of the superfields φ i ) by 4π 2 C(R) i j δ 4 (Q) or 4π 2 C(R) i j δ 4 (Q + K), respectively. Note that no replacement should be made for the nonsingular Pauli-Villars propagators. The above procedure gives the result ∆ matter After calculating the integrals of the δ-functions and some transformations we obtain the expression which exactly coincides with the one found in Ref. [71] with the help of a different prescription. The solid lines in the supergraphs B2 and B3 can also stand for the propagators of the Pauli-Villars superfields ϕ a . In this case the calculation of the vacuum superdiagrams with an insertion of θ 4 (v B ) 2 multiplied by the operator (38) for an arbitrary value of ξ 0 gives In this case all matter propagators are massive and do not produce singularities. Therefore, it is only necessary to make the replacement 1/K 2 → 4π 2 C 2 δ 4 (K), after which the contribution of the Pauli-Villars superfields ϕ a to the function (34) takes the form This expression again coincides with the one found in Ref. [71] with the help of a different technique. 10 The contribution of the vacuum supergraphs B4, B5, B6, and B7 is written in the form We see that all bad terms containing inverse momenta to the fourth power really cancel each other, although they are present in expressions for the separate supergraphs. This agrees with the general argumentation of section 5. Also all terms dependent on the gauge parameter ξ 0 cancel each other. To find a contribution to the function (34) we should replace one of the gauge or ghost inverse squared momenta by the corresponding δ-function multiplied by 4π 2 C 2 and sum up all expression thus obtained. The result is After calculating the integrals of the δ-functions this expression can be rewritten as Again we see that the result coincides with the one obtained by the method of Ref. [61] in Ref. [71]. This confirms the correctness of the argumentation used in this paper for deriving the NSVZ equation.
The overall two-loop expression for the β-function reads as and coincides with the one obtained in Ref. [71]. Certainly, this expression satisfies the NSVZ equations (28) and (120) and agrees with the previous calculations (first made in Ref. [91]).

The three-loop approximation: supergraphs with Yukawa vertices
Next, we will verify the argumentation of this paper on the example of the three-loop contribution to β-function which contains the Yukawa couplings. It is generated by the supergraphs B8 -B11 presented in Fig. 4. The direct calculation of the two-point superdiagrams obtained from them by attaching two external lines of the background gauge superfield V has been made in Refs. [48,49] in the Feynman gauge ξ 0 = 1. Subsequently, the result has been reobtained by the method described in section 4.1 in Ref. [61]. Now we will demonstrate how the modification proposed in this paper works in this case. First, we consider the supergraph B8 which is quartic in the Yukawa couplings. The result for the modified vacuum supergraph is given by the expression which does not contain bad terms proportional to the inverse momenta to the fourth power.
As earlier, to construct the corresponding contribution to the function (34), we should sum all expressions obtained from (145) by replacing squared inverse momenta multiplied by δ-symbols coming from the corresponding propagator by momentum δ-functions multiplied by 4π 2 C(R). The tensor structure of the resulting factors can easily be viewed from the structure of the graph under consideration. After constructing and calculating the integrals of the δ-functions we obtain the required contribution which exactly coincides with the corresponding result of Refs. [48,61] found by different methods.
According to [61] the properly modified vacuum supergraph B9 is given by the expression where K µ corresponds to the propagator of the quantum gauge superfield and, following Ref.
[49], we use the notation Again, Eq. (147) does not contain bad terms. Next, we proceed according to the algorithm of section 4.2. In this case the relevant replacements are 1/K 2 → 4π 2 C 2 δ 4 (K) and δ j i /P 2 → 4π 2 C(R) i j δ 4 (P ), where P µ stands for momenta of the matter superfields, namely, Q µ , (Q + K) µ , (Q − L) µ , and (Q + K − L) µ . After some transformations involving the identities which follow from Eq. (10), we get the result which agrees with the one found in [49,61] by different methods. The supergraph B10 is given by the expression where K µ denotes the momentum of the quantum gauge superfield propagator and, again following Ref. [49], In this case the algorithm of section 4.2 produces the expression which also coincides with the results of Refs. [49,61]. The expression for the last supergraph B11 reads as where the momentum of the gauge superfield propagator is denoted by K µ and Replacing the squared inverse momenta and the corresponding δ-symbols according to the prescription given in section 4.2 we obtain the part of the function (34), Again, it agrees with the calculations made previously.

The three-loop approximation: supergraphs with ghost loops
All three-loop vacuum supergraphs containing loops of the Faddeev-Popov ghosts have been calculated in Ref. [50]. These supergraphs are presented in Fig. 5. Note that in this approximation the first nonlinear term in the function F(V ) (see Eq. (8)) is essential. It generates the vertex present in the supergraph B21, see Ref. [84] for details. This vertex is denoted by a cross. Expressions for the supergraphs in Fig. 5 contain 1/K 4 singularities, where K µ is the Euclidean momentum of the quantum gauge superfield. As in the two-loop approximation, these singularities should disappear after adding the purely gauge vacuum supergraphs. However, the sum of the three-loop diagrams containing only the gauge propagators has not yet been calculated with the higher covariant derivative regularization. Nevertheless, the sum of the supergraphs with ghost loops does not contain any singularities proportional to the inverse ghost or matter momenta to the fourth power. This implies that it is possible to compare the sums of singularities coming from the cuts of ghost and matter lines with the corresponding anomalous dimensions γ c (α 0 , λ 0 , Y 0 ) and (γ φ ) i j (α 0 , λ 0 , Y 0 ). With the help of the method described in section 4.1 this has been done in Ref. [50]. Also it is possible to use these results for checking the modification of the algorithm discussed in section 4.2. This is made as follows: 1. First, we consider a vacuum supergraph containing a ghost loop (or loops) with an insertion of θ 4 (v B ) 2 and calculate it using the D-algebra.
2. Next, we replace one of inverse squared momenta (say, 1/Q 2 ) of a ghost propagator by 4π 2 C 2 δ 4 (Q). If matter propagators are present in the considered supergraph, then δ j i /Q 2 coming from a matter propagator is replaced by 4π 2 C(R) i j δ 4 (Q). All expressions obtained after the replacements of one propagator should be summed.
3. After calculating the integrals of the δ-functions the results are compared with the corresponding contributions to the ghost (or matter) anomalous dimensions. These contributions are given by the sums of all two-point supergraphs which are produced by all possible cuts of ghost (or matter) propagators in the considered supergraph, see Fig. 5.
We have done this for all three-loop vacuum supergraphs presented in Fig. 5. As a result, we have obtained in agreement with Eqs. (87) and (98). Note that this verification is different from the one made in Ref. [50], because the algorithm used for constructing the contribution to the function (34) is different.

N = 1 SQED in the three-loop approximation
As one more example for checking the method proposed in section 4.2 we can consider N = 1 SQED with N f flavors. In this case the gauge group is U (1), r = 1, C 2 = 0, and C(R) i j → δ αβ ·1 2 , where α, β = 1, . . . , N f and 1 2 denotes an identity matrix of the size 2 × 2. This theory does not contain Yukawa terms triple in chiral superfields. Therefore, the calculations in this case have been done for F (x) = 1. Although N = 1 SQED is a particular case of N = 1 supersymmetric gauge theories considered earlier, this example is not trivial, because for this theory the complete expression for the three-loop β-function has explicitly been calculated with the higher derivative regularization. That is why it is possible to perform one more nontrivial test of the method proposed in this paper. The gray circles denote insertions of the one-loop polarization operator of the quantum gauge superfield, see Ref. [87] for details.
Using the results of Ref. [90] the sum of two-and three-loop vacuum supergraphs modified by an insertion of θ 4 (v B ) 2 and multiplied by −2π/V 4 · d/d ln Λ in the general ξ 0 -gauge is written as and does not contain inverse momenta to the fourth power. According to the method discussed in section 4.2, to find a contribution to the function (34), we should sum the expressions obtained from Eq. (162) by replacing one of 1/P 2 by 4π 2 δ 4 (P ), where P µ is a momentum of the (massless) matter superfields. Note that it is not necessary to take into account the gauge superfield propagators, because now C 2 = 0, and nonsingular propagators of the massive Pauli-Villars superfields. In the first two terms of Eq. (162) the momentum of the quantum gauge superfield is denoted by K µ , while in the last term the momenta of the quantum gauge superfield propagators are K µ and L µ . Then the above described procedure gives Again, this expression correctly reproduces the result obtained earlier by different methods, see Refs. [51,90]. Therefore, for all considered supergraphs the results obtained with the help of the technique discussed in section 4.2 (which is a certain modification of the one proposed in [61]) coincided with the expressions found earlier by different methods. This confirms the correctness of the method which was used in this paper for the all-loop perturbative derivation of the exact NSVZ β-function.

Conclusion
In this paper we have finished the all-order perturbative derivation of the exact NSVZ βfunction (1) for non-Abelian N = 1 supersymmetric gauge theories which was started in Refs. [2,61,62]. Its main ingredient is the higher covariant derivative regularization, which allows revealing some interesting features of quantum corrections in these theories. For instance, according to Ref. [62], with this regularization all loop integrals giving the β-function defined in terms of the bare couplings are integrals of double total derivatives with respect to loop momenta. These integrals do not vanish due to singularities, which appear when double total derivatives act on massless propagators. The sum of the singularities produces all contributions to the β-function starting from the two-loop approximation. (The one-loop quantum corrections should be considered separately. This has been done in Ref. [75].) It is possible to divide singular contributions into three groups depending on a propagator on which the double total derivatives act. Namely, they can act on the matter superfield propagators, on the propagators of the Faddeev-Popov ghosts, and on the propagators of the quantum gauge superfield. Qualitatively, this can be interpreted as a cutting of the corresponding internal line in a certain vacuum supergraph. All such cuts give a set of superdiagrams contributing to the anomalous dimension of the corresponding quantum superfield. From the other side, attaching two external gauge lines of the background gauge superfield to the considered supergraph we obtain a set of superdiagrams contributing to the β-function. The NSVZ equation in the form (5) relates this contribution to the β-function to the parts of the anomalous dimensions of quantum superfields coming from the superdiagrams produced by cuts of internal lines. Note that the factorization into double total derivatives allows to calculate analytically one of loop integrals, so that the β-function in a certain order is really related to the anomalous dimensions in the previous order. The proof of this fact has been done in this paper. We have demonstrated that the sums of singularities coming from the cuts of certain propagators are really equal to the corresponding terms in Eq. (5). (Note that in Ref. [62] this has been done for the matter and Faddeev-Popov ghost superfields. The results of this paper obtained by a different method are the same. However, in the present paper we have also calculated the all-loop sum of singularities produced by cuts of the gauge propagators.) Thus, Eq. (5) for RGFs defined in terms of the bare couplings is proved in all orders in the case of using the higher covariant derivative regularization. 11 The original NSVZ β-function (1) for these RGFs can be obtained with the help of the non-renormalization theorem for the triple gauge-ghost vertices proved in [2]. Note that RGFs defined in terms of the bare couplings are scheme-independent for a fixed regularization, so that both these equations are valid for any renormalization prescription supplementing the higher covariant derivative regularization.
Taking into account that in the HD+MSL scheme RGFs defined in terms of the renormalized couplings coincide with the ones defined in terms of the bare couplings up to the renaming of arguments, we conclude that for standardly defined RGFs one of the NSVZ schemes is given by the HD+MSL prescription in all orders. By other words, to obtain the NSVZ equation in all loops, one should regularize a theory by higher covariant derivatives and include into renormalization constants only powers of ln Λ/µ (or, equivalently, set all finite constants to 0).
As we have already mentioned, explicit calculations exactly confirm the statements discussed above even in the approximations, where the dependence on a regularization and a renormalization prescription becomes essential, see, e.g., [48,49,50,90]. Also the results of this paper confirm the correctness of the expression for the three-loop β-function of a general N = 1 supersymmetric gauge theory with matter superfields and a simple gauge group derived in Ref. [92] from the NSVZ equation.