The (cid:12) -function of N = 1 supersymmetric gauge theories regularized by higher covariant derivatives as an integral of double total derivatives

,


Introduction
Ultraviolet divergences in supersymmetric theories are restricted by some nonrenormalization theorems. According to one of them, N = 4 supersymmetric Yang-Mills (SYM) theory is finite in all orders [1][2][3][4]. Divergencies in N = 2 theories exist only in the one-loop approximation [1,4,5], so that it is even possible to construct finite N = 2 supersymmetric theories by choosing a gauge group and a matter representation in such a way that the one-loop divergencies cancel [6]. All these non-renormalization theorems can be derived [7,8] from the equation which relates the β-function of N = 1 supersymmetric gauge theories with the anomalous dimension of the matter superfields [9][10][11][12] β(α, λ) = − where α is the gauge coupling constant and λ denotes the Yukawa couplings. Note that so far we do not specify the definitions of the renormalization group functions (RGFs) and what couplings are considered as their arguments. Eq. (1.1) called the exact NSVZ β-function can also be considered as a non-renormalization theorem in addition to the JHEP10(2019)011 well-known statement that the superpotential in N = 1 supersymmetric theories is not renormalized [13]. According to one more non-renormalization theorem derived in [14], the triple ghost-gauge vertices in N = 1 supersymmetric gauge theories are finite in all orders. 1 With the help of this non-renormalization theorem the exact NSVZ β-function can be equivalently rewritten in a new form [14], which relates the β-function to the anomalous dimensions of the quantum gauge superfield (γ V ), of the Faddeev-Popov ghosts (γ c ), and of the matter superfields ( γ φ i j ).
Some NSVZ-like relations can be written for other theories. For example, in theories with softly broken supersymmetry an analogous equation describes the renormalization of the gaugino mass [18][19][20]. Also it is possible to construct the NSVZ-like equations for the Adler D-function in N = 1 SQCD [21,22] and even for the renormalization of the Fayet-Iliopoulos term in two-dimensional N = (0, 2) supersymmetric models [23].
Various derivations of the exact NSVZ β-function involve general arguments based on the analysis of the instanton contributions [7,9], anomalies [10,12,24], and nonrenormalization of the topological term [25]. However, a direct perturbative verification of eq. (1.1) in all orders appeared to be a highly non-trivial problem. Even to start solving this problem, one should first pay attention to some important subtleties related to the regularization, quantization, and renormalization.
Really, the calculations of quantum corrections made in the DR-scheme (that is with the help of dimensional reduction [26] supplemented by the modified minimal subtractions [27]) in refs. [28][29][30][31][32] demonstrate that the NSVZ relation is not valid for this renormalization prescription. However, the difference can be explained by the scheme dependence of the NSVZ relation, which is described by the general equations derived in [33,34]. Namely, it is possible to tune the renormalization scheme in such a way that the NSVZ equation will take place [28][29][30]. 2 It is important that this possibility is highly non-trivial due to some scheme-independent equations following from the NSVZ relation [34,36]. Nevertheless, at present there is no general all-loop prescription giving the NSVZ scheme in the case of using the regularization by dimensional reduction.
The NSVZ renormalization prescription can be naturally formulated in all loops if N = 1 supersymmetric gauge theories are regularized by the higher covariant derivative method [37,38] in the supersymmetric version [39,40]. The matter is that using of this regularization reveals the underlying structure of the loop integrals responsible for appearing the NSVZ relation. Namely, in this case the integrals giving the β-function defined in terms of the bare couplings appear to be integrals of double total derivatives with respect to loop momenta. 3 This was first noted in calculating quantum corrections for N = 1 su-
The integrals of double total derivatives do not vanish due to the identity where Q is an Euclidean momentum. The δ-function reduces the number of loop integrations by 1, so that in the Abelian case an L-loop contribution to the β-function appears to be related to an (L − 1)-loop contribution to the anomalous dimension of the matter superfields. The sum of singularities in the Abelian case was calculated in [54,55], where it was expressed in terms of the anomalous dimension of the matter superfields. The relation between the β-function and the anomalous dimension obtained in this way is nothing else than the NSVZ equation for RGFs defined in terms of the bare couplings. Thus, at least in the Abelian case, it naturally appears in the case of using the higher derivative regularization. Note that the RGFs defined in terms of the bare couplings are scheme independent if a regularization is fixed (see, e.g., [57]), so that the NSVZ equation for these RGFs is valid for an arbitrary renormalization prescription. 4 In the non-Abelian case the situation is much more complicated. Eq. (1.1) relates an L-loop contribution to the β-function to the anomalous dimension of the matter superfields in all previous orders. That is why it is more probable that it is eq. (1.2) that originally appears in the perturbative calculations. Moreover, unlike eq. (1.1), eq. (1.2) can be visualized in the same way as in the Abelian case (see refs. [44,50]). Namely, starting from a supergraph without external lines, it is possible to obtain a contribution to the β-function by attaching two external lines of the background gauge superfield and contributions to the anomalous dimensions by cutting internal lines. Thus obtained contributions are related by eq. (1.2).
The similarity between eq. (1.2) and the Abelian NSVZ equation [58,59] allows suggesting that the factorization of integrals into double total derivatives also produces the NSVZ equation in the non-Abelian case. This guess was confirmed by numerous calculations in the lowest loops, see, e.g., [47,51,53,60]. This implies that all higher order corrections to the β-function (starting from the two-loop approximation) appear from the δ-singularities. Therefore, to derive the NSVZ relation in the non-Abelian case (for RGFs defined in terms of the bare couplings with the higher covariant derivative regularization), it is necessary only to sum singular contributions and to prove that they give the sum of the anomalous dimensions in the right hand side of eq. (1.2). If this is really so, then the NSVZ scheme for

JHEP10(2019)011
RGFs defined in terms of the renormalized couplings is given by the so-called HD+MSL prescription [14] exactly as in the Abelian case [34,36,57]. 5 This means that the theory is regularized by higher covariant derivatives supplemented by the minimal subtractions of logarithms, when only powers of ln Λ/µ are included into the renormalization constants. 6 The paper is organized as follows: in section 2 we formulate the theory under consideration in N = 1 superspace, regularize it by higher covariant derivatives, and describe the quantization. Also in this section we introduce some auxiliary constructions, which will be needed for the investigation of the loop integrals giving the β-function. RGFs defined in terms of the bare couplings are introduced in section 3. In this section we also present the β-function and the NSVZ relation for it in the form which is mostly convenient for the analysis. In section 4 we demonstrate that the β-function defined in terms of the bare couplings is given by integrals of double total derivatives with respect to loop momenta. Here we also describe the method which allows to construct these integrals in a simple way. This method is applied for calculating the three-loop contribution to the β-function containing the Yukawa couplings in section 5. In particular, we demonstrate that the result exactly coincides with the one obtained in ref. [53] with the help of the standard supergraph calculation.
2 N = 1 supersymmetric gauge theories: regularization, quantization, and auxiliary parameters It is convenient to describe N = 1 supersymmetric gauge theories using N = 1 superspace with the coordinates (x µ , θ), where θ is an auxiliary anticommuting Majorana spinor. In this case N = 1 supersymmetry of the theory is manifest. Moreover, it becomes possible to perform the quantization and calculate quantum corrections in a manifestly N = 1 supersymmetric way [64][65][66]. At the classical level the considered theory in the massless limit is described by the action where V is the Hermitian gauge superfield and φ i are the chiral matter superfields in a representation R of a gauge group G which is assumed to be simple. In the classical theory (2.1) the supersymmetric gauge superfield strength is defined as W a ≡D 2 e −2V D a e 2V /8. The gauge coupling constant is defined as α = e 2 /4π, and the Yukawa couplings are denoted by λ ijk . Note that at the classical level we do not distinguish between bare and renormalized couplings. This difference is essential in the quantum theory. Below, considering the quantum theory, we will denote the bare couplings by α 0 = e 2 0 /4π and λ ijk 0 , while the renormalized couplings will be denoted by α and λ ijk . 5 HD+MSL prescription also gives the NSVZ-like schemes for the Adler D-function [52] and for the renormalization of the photino mass in softly broken N = 1 SQED [61]. 6 This NSVZ scheme is not unique [62]. For example, in N = 1 SQED the on-shell scheme is also NSVZ [63].

JHEP10(2019)011
Below t A and T A are the generators of the fundamental representation and the representation R, respectively. These sets of generators satisfy the conditions We will always assume that tr(T A ) = 0. Also we will use the notation (The generators of the adjoint representation are expressed in terms of the structure constants as ( Under the condition parameterized by a Lie algebra valued chiral superfield A. To quantize the theory (2.1), it is also necessary to take into account that the quantum gauge superfield is renormalized in a nonlinear way [67][68][69] (see also refs. [70,71]). The necessity of this nonlinear renormalization has been demonstrated by explicit calculations in refs. [72,73]. Moreover, the two-loop calculation of the Faddeev-Popov ghost anomalous dimension in [74] showed that without this nonlinear renormalization the renormalization group equations are not satisfied. Thus, it is really needed for quantum calculations. To take into account the nonlinear renormalization, following ref. [68], we substitute the gauge superfield V by the function F (V ) in the action functional. Moreover, it is necessary to replace e and λ by the bare couplings e 0 and λ 0 , respectively.
For obtaining a manifestly gauge invariant effective action we will use the background field method [75][76][77] formulated in N = 1 superspace [1,64]. A distinctive feature of the background field method in the supersymmetric case is the nonlinear background-quantum splitting which in the considered case can be implemented by the substitution where in the right hand side V and V are the quantum and background gauge superfields, respectively. 7 In this case the quantum gauge superfield satisfies the constrain V + = e −2V V e 2V . Due to the background-quantum splitting the gauge invariance produces two different types of gauge transformations. Under the background gauge symmetry the superfields of the theory change as The standard form of the background quantum splitting is e 2F (V ) → e Ω + e 2F (V ) e Ω , the background gauge superfield being defined by the equation e 2V = e Ω + e Ω . However, after the change of variables V → e −Ω + V e Ω + in the generating functional we arrive to eq. (2.6).

JHEP10(2019)011
This invariance remains unbroken at the quantum level and becomes a manifest symmetry of the effective action. Alternatively, the quantum gauge invariance is broken by the gauge fixing procedure. It is convenient to introduce the background supersymmetric covariant derivatives ∇ a and∇ȧ and the gauge supersymmetric covariant derivatives ∇ a and∇ȧ defined by the equations Note that for the purposes of this paper it is more convenient to use a different representation for them in comparison with refs. [74,78]. In the representation (2.9) the covariant derivatives ∇ a and∇ȧ should act on a function X which changes as X → e −A + X. In this case they transform in the same way under both background and quantum transformations. This is also valid for the background covariant derivatives ∇ a and∇ȧ, but only in the case of the background gauge transformations. If we use the background field method and take into account the nonlinear renormalization of the quantum gauge superfield, then the gauge superfield strength is defined as (2.11) Below we will also need some auxiliary parameters. The coordinate-independent complex parameter g describes the continuous deformation of the original theory (corresponding to g = 1) into the theory in which quantum superfields interact only with the background gauge superfield (corresponding to g → 0). This parameter is introduced by making the substitutions Then, it is easy to see that an L-loop contribution to the two-point Green function of the background gauge superfield is proportional to (gg * ) L−1 . Also we introduce the auxiliary chiral superfield 8 g(x, θ). It is added to g in such a way that all quantum corrections containing g will actually depend on the (coordinatedependent) combination

JHEP10(2019)011
Now, let us include the parameters g and g into the classical action. For this purpose we write all terms containing the quantum gauge superfield as integrals over d 4 x d 4 θ ≡ d 8 x with the help of eq. (2.10). After this we modify the result by introducing the auxiliary parameters in the following way: where the integration measures are Note that we do not include the superfield g in the first term of eq. (2.14), which does not contain the quantum gauge superfield V . This allows to avoid breaking of the background gauge invariance (2.7). However, the action (2.14) is invariant under the quantum gauge transformations (2.8) only if g = 0 (but for an arbitrary value of the coordinate independent parameter g). Nevertheless, it is not important, because the parameter g is auxiliary and actually we are interested only in the cases when g = 0, 1 and g = 0.
The most important ingredient needed for deriving the NSVZ β-function for RGFs defined in terms of the bare couplings is the higher covariant derivative regularization [37,38]. In this paper we will use the version similar to the one considered in ref. [78] with some modifications appearing due to the presence of the auxiliary parameters and the function F (V ). To regularize a theory by higher covariant derivatives, at the first step, it is necessary to add a higher derivative term S Λ to its action. As a result, propagators will contain higher degrees of momenta that, in turn, leads to the finiteness of the regularized theory beyond the one-loop approximation [82]. In the case g = 0 the regularized action S reg = S + S Λ invariant under both background and quantum gauge transformations can be constructed as where the higher derivative regulators R(x) and F (x) are functions rapidly growing at infinity which satisfy the conditions R(0) = F (0) = 1. In eq. (2.16) and below the subscript Adj means that remains unbroken. This can be done similarly to constructing the action (2.14). However, it is more difficult due to the presence of the function R(x). We present this function in the form Then the regularized action can be written as It is important that this action is invariant under the background gauge transformations, but the quantum gauge invariance exists only for g = 0. In this case the action (2.19) is reduced to eq. (2.16). Moreover, all terms containing the quantum superfields depend on auxiliary parameters only in the combination g = g + g. (The first term, which depends on the constant g and does not depend on the superfield g, contains only the background gauge superfield.) To obtain a manifestly gauge invariant effective action, it is necessary to use a gauge fixing term invariant under the background transformations (2.7). Taking into account that a higher derivative regulator should be also inserted into this term [78], the gauge fixing action can be chosen as Certainly, the quantization procedure also requires to introduce the Faddeev-Popov action. The Faddeev-Popov ghosts and the corresponding antighosts in the supersymmetric case are described by the chiral superfields c A andc A , respectively. The action for them obtained in a standard way takes the form (2.21) In the case of using the background superfield method it is also necessary to take into account the Nielsen-Kallosh ghost action

JHEP10(2019)011
Here the Nielsen-Kallosh ghosts b are chiral anticommuting superfields in the adjoint representation, which interact only with the background gauge superfield. The arrow points out that the parameters g and e 0 can be excluded from the Nielsen-Kallosh action by the change of variables b → e 0 gb; b + → e 0 g * b + in the generating functional. (It is easy to see that the corresponding determinant is equal to 1.) After the gauge fixing procedure the quantum gauge transformations (2.8) are no longer a symmetry of the total action (that, in particular includes the gauge fixing term and ghosts). The total action is invariant under the BRST transformations [83,84]. In N = 1 superspace the BRST transformations have been formulated in ref. [67]. For the theory considered in this paper the BRST invariance is a symmetry of the action only in the case g = 0, but for an arbitrary value of the coordinate independent parameter g.
As we mentioned above, the one-loop divergences cannot be regularized by adding the higher derivative term to the action. For this purpose it is necessary to supplement the higher derivative method by the Pauli-Villars regularization which is introduced by inserting the Pauli-Villars determinants into the generating functional [85]. According to refs. [78,86], to cancel the one-loop divergences appearing in supersymmetric gauge theories, one should introduce three chiral Pauli-Villars superfields ϕ a with a = 1, 2, 3 in the adjoint representation of the gauge group, and chiral superfields Φ i in a certain representation R PV which admits a gauge invariant mass term. The superfields ϕ a cancel one-loop divergences coming from the loops of the quantum gauge superfield, of the Faddeev-Popov ghosts and of the Nielsen-Kallosh ghosts. The superfields Φ i cancel the one-loop divergences coming from the matter loop. This occurs if the generating functional is defined as where Dµ denotes the measure of the functional integration and c = T (R)/T (R PV ). The sources are included into 9 The Pauli-Villars determinants are constructed as

JHEP10(2019)011
and M jk M * ki = M 2 δ j i . (We assume that the representation R PV is chosen in such a way that this condition can be satisfied. For example, it is possible to use the adjoint representation.) To obtain a regularized theory with a single dimensionful parameter, it is necessary to require that the Pauli-Villars masses M ϕ and M should be proportional to the parameter Λ, It is important that we consider a regularization for which a ϕ and a do not depend on couplings. The effective action is standardly defined as the Legendre transform of the generating functional W = −i ln Z for connected Green functions, where the sources should be expressed in terms of (super)fields from the equations (2.30)

Renormalization and RGFs defined in terms of the bare couplings
In this section we present the β-function defined in terms of the bare couplings in a form which is the most convenient for proving the factorization of the corresponding loop integrals into integrals of double total derivatives. This factorization is an important step towards constructing the all-loop perturbative derivation of the exact NSVZ β-function.
That is why in this section we also rewrite the NSVZ relation (1.2) in such a form that can be used as a starting point of this derivation. To find the β-function defined in terms of the bare couplings, we consider the twopoint Green function of the background gauge superfield. Note that in our conventions the term "two-point" in particular means that the auxiliary superfield g is set to 0, but the dependence on the parameter g is kept. It is easy to see that the considered Green function depends on g, α 0 , λ 0 , and λ * 0 only via the combinations gg * α 0 and gg * λ ijk 0 λ * 0mnp . (For simplicity, below we will denote the latter one by gg * λ 0 λ * 0 .) Really, in the case g = 0 the total action depends on gg * α 0 , gλ 0 and g * λ * 0 . However, the numbers of λ 0 and λ * 0 in any supergraph contributing to the considered Green function are equal. Therefore, the Yukawa couplings enter it only in the combination gg * λ 0 λ * 0 . Similar arguments also work for the two-point Green functions of the quantum gauge superfield, of the Faddeev-Popov ghosts, and for the two-point Green function φ * i φ j of the matter superfields. Below we will use the notation ρ ≡ |g| 2 = gg * , (3.1) so that the above mentioned two-point Green functions actually depend on ρα 0 and ρλ 0 λ * 0 . Due to the background gauge invariance the two-point Green function of the background gauge superfield is transversal and (in the massless limit) can be written as

JHEP10(2019)011
where the supersymmetric transversal projection operator is defined by the equation With the help of the Slavnov-Taylor identities [87,88] (and some other similar equations) it is possible to prove that quantum corrections to the two-point Green function of the quantum gauge superfield are also transversal, Also we will need the two-point Green functions of the Faddeev-Popov ghosts and of the matter superfields, Renormalized couplings α, λ and the renormalization constants pressed in terms of α and λ in the limit Λ → ∞. Note that due to the non-renormalization of the superpotential [13] the renormalized Yukawa couplings are related to the bare ones by the equation Similarly, due to the non-renormalization of the triple ghost-gauge vertices [14] the renormalization constants can be chosen in such a way that We will always assume that the renormalization constants satisfy eqs. (3.7) and (3.8).
(Certainly the renormalization constants are not uniquely defined [89], and these constrains partially fix an arbitrariness in choosing a subtraction scheme.) It is important that in the non-Abelian case the quantum gauge superfield is renormalized in a nonlinear way [67][68][69]. The non-linear renormalization can be realized as a linear renormalization of an infinite set of parameters. For example, in the lowest approximation it is possible to present the function F (V ) in the form where y 0 is a new bare parameter and (3.10)

JHEP10(2019)011
Then, the result for the nonlinear renormalization obtained in [72,73] can be equivalently written in the form where ξ is the renormalized gauge parameter and k 1 is a finite constant which appears due to the arbitrariness in choosing a subtraction scheme. The explicit calculation of ref. [74] demonstrated that the renormalization group equations cannot be satisfied without introducing the parameter y 0 (or, possibly, implementing the nonlinear renormalization by some different way). Certainly, in higher orders an infinite set of parameters similar to y 0 is needed. All these parameters are similar to the gauge fixing parameter ξ 0 , because by a proper change of variables in the generating functional it is possible to prove that a nonlinear renormalization is equivalent to a nonlinear change of a gauge [67]. That is why below we will include the gauge fixing parameter and the parameters of the nonlinear renormalization inside the function F (V ) into a single set The corresponding renormalized values will be denoted by Y = (ξ, y, . . .).
We believe that the NSVZ relation is valid for RGFs defined in terms of the bare couplings in the case of using the higher covariant derivative regularization. These RGFs are defined by the equations and do not depend on a renormalization prescription for a fixed regularization [57]. It is easy to see that RGFs defined in terms of the bare couplings can be obtained by differentiating the corresponding Green functions. For example, the β-function defined in terms of the bare couplings can be constructed by differentiating the quantum corrections in the twopoint Green function of the background gauge superfield in the limit of the vanishing external momentum, (3.14) Note that the term 1/(gg * α 0 ) appears in the function d −1 in the tree approximation and corresponds to the first term in eq. (2.19). The limit p → 0 is needed for removing terms proportional to (p/Λ) k , where k is a positive integer. The equality follows from the finiteness of the function d −1 expressed in terms of the renormalized couplings.

JHEP10(2019)011
It is well known that for g = 1 the β-function can be presented as the series where the (Y 0 -independent) coefficient is obtained by calculating the one-loop contribution to the β-function. (For the considered regularization the details of this calculation can be found in [78].) For g = 1 it is easy to see that the L-loop contribution to the β-function is proportional to gg * L+1 = ρ L+1 . Therefore, the dependence of the expression β(ρα 0 , ρλ 0 λ * . . If we consider g and g * as independent variables, then Consequently, where +0 means that ρ = 0, but ρ → 0. Taking into account that the limit ρ → 0 corresponds to the theory in which quantum superfields interact only with the background gauge superfield, so that nontrivial quantum corrections exist only in the one-loop approximation, we obtain Therefore, the β-function defined in terms of the bare couplings (for the original theory which corresponds to g = 1) can be calculated with the help of the equation Due to the finiteness of the functions dimensions of the quantum superfields can also be related to the corresponding Green functions by the equations In the one-loop order these anomalous dimensions contain terms proportional to α 0 and λ 0 λ * 0 (the latter ones appear only in γ φ i j ), (3.24) and the terms corresponding to the L-loop approximation are proportional to gg * L = ρ L . Using this fact, from the identity (3.18) we obtain This implies that for deriving the NSVZ relation (1.2) it is sufficient to prove that to this equation with the help of eqs. (3.14) and (3.21) -(3.23). In eq. (3.26) the derivative with respect to ln Λ is very important, because it removes infrared divergences which could appear in the limit of the vanishing external momentum. Explicit loop calculations (e.g., in refs. [51,53]) demonstrate that loop integrals written without d/d ln Λ are not well defined, while after the differentiation all bad terms disappear.
The derivatives with respect to g and g * are not so important and can be excluded from eq. (3.26). Certainly, in this case it is necessary to add the constant corresponding to the one-loop contribution, For g = 1 this identity was first suggested in ref. [14]. However, for deriving the NSVZ relation in all loops it is more preferable to use eq. (3.26).
The left hand side of eq. (3.26) can be constructed starting from the expression for the two-point Green function of background gauge superfield (3.2). To extract the function d −1 , it is convenient to make the formal substitution where θ 4 ≡ θ a θ aθȧθȧ .
where v A 0 = const and X µ = (x i , ix 0 ) are the Euclidean coordinates. The corresponding Euclidean momenta are denoted by From eq. (3.31) we see that v A (P ) is essentially different from 0 only in a small region of the size 1/R → 0. This implies that substituting the functions (3.30) into eq. (3.2) we automatically obtain the limit P → 0 (or, equivalently, p → 0), which is needed for constructing RGFs defined in terms of the bare couplings. Let us consider quantum corrections encoded in the expression where S total includes the usual action, the gauge fixing term, and the ghost actions. (Certainly, the terms proportional to Λ −k , where k is a positive integer, should be omitted).
Then we consider a part of ∆Γ corresponding to the two-point Green function of the background gauge superfield. Performing the Wick rotation and making the substitution (3.29), after some transformations, in the limit R → ∞ we obtain where we have introduced the notation Thus, we see that the substitution (3.29) allows extracting the β-function defined in terms of the bare couplings from the considered part of the effective action in the case of using the higher covariant derivative regularization. (In the case of using the dimensional reduction one should be much more careful, see [41,42] for details.) Differentiating eq. (3.33) with respect to the parameters g and g * and multiplying the result by the factor 2π/V 4 , we obtain the left hand side of eq. (3.26). In turn, the derivatives with respect to the coordinate-independent parameters g and g * can be expressed in terms of the derivatives with respect to the chiral superfield g and the antichiral superfield g * , respectively. Really, all terms in the action containing quantum superfields depend only on the combinations g and g * , see eqs. which depends on g and g * in a different way is the first term in eq. (2.19), but it does not affect quantum corrections and does not enter ∆Γ. Therefore, it is possible to relate the derivatives of ∆Γ with respect to g and g * to the derivatives with respect to g and g * , Thus, to derive the NSVZ relation, it is sufficient to prove the identity V denotes a part of Γ which is quadratic in the background gauge superfield and does not contain the other superfields except for g. Note that writing eq. (3.37) we took into account that S Note that here we do not set the auxiliary external superfields g and g * to 0, because eq. (3.37) contains the derivatives with respect to these superfields. In this paper we will consider only N = 1 supersymmetric gauge theories with a simple gauge group. In this case it is easy to see that any invariant tensor I AB should be proportional to δ AB . 10 Therefore, for simple gauge groups With the help of eqs. (3.38) and (3.39) for a simple gauge group it is possible to rewrite eq. (3.37) in the form mostly convenient for proving, namely, According to the above discussion, for the theory regularized by higher covariant derivatives this equation is equivalent to the NSVZ relations (1.1) and (1.2) for RGFs defined in terms of the bare couplings. Below we will prove that the left hand side of eq. (3.40) is given by integrals of double total derivatives. 10 The considered invariant tensor satisfies the equation [T A Adj , I] = 0, so that it commutes with all generators of the adjoint representation. For a simple group the adjoint representation is irreducible. Therefore, IAB should be proportional to δAB.

JHEP10(2019)011
4 The β-function as an integral of double total derivatives

The Slavnov-Taylor identity for the background gauge invariance
The background gauge invariance is a manifest symmetry of the theory under consideration (even in the presence of the auxiliary superfield g). At the quantum level symmetries are encoded in the Slavnov-Taylor identities [87,88]. The Slavnov-Taylor identity corresponding to the background gauge transformations constructed in this section is a very important ingredient for the all-loop proof of the factorization into double total derivatives. This identity is derived by standard methods, namely, it is necessary to make the change of variables in the functional integral (2.23), which does not change the generating functional Z. This change of variables coincides with the background gauge transformations of the quantum superfields. Due to the background gauge invariance, the total gauge fixed action and the Pauli-Villars determinants remain unchanged if the background gauge superfield is also modified as However, the source term S sources transforms nontrivially. This implies that in the linear order in A the invariance of the generating functional W = −i ln Z under the change of variables (4.1) can be expressed by the equation where the variations of various superfields under the infinitesimal background gauge transformations are written as 11 where B is a function(al) depending on the superfields of the theory. 11 The expression for δV = δV B t B is obtained in the standard way from the identity 0 = δ[V , e 2V ].

JHEP10(2019)011
Rewriting eq. (4.4) in terms of (super)fields, we obtain the equation which expresses the manifest background gauge invariance of the effective action, (4.7) It is important that in this equation (super)fields are not set to 0, so that this equation encodes an infinite set of identities relating Green functions of the theory. That is why we will call it the generating Slavnov-Taylor identity.
Considering A and A + as independent variables and differentiating eq. (4.7) with respect to A A we obtain where the matrix [f (X) Adj ] AB is defined by the equation Expressing the generators of the adjoint representation in terms of the structure constants it is possible to rewrite the generating Slavnov-Taylor identity (4.7) corresponding to the background gauge symmetry in the form where the operatorÔ A is given by the expression To verify eq. (4.10), it is necessary to take into account that a derivative with respect to a chiral superfield is also chiral and use the identity valid for an arbitrary chiral superfield φ.
It is important that due to eq. (4.10) the effective action satisfies the equation whereŌ Ȧ a ≡ −Dȧ This can be verified with the help of the equality Therefore, taking into account that f AAC = 0, after the differentiation we see that Note that here all fields (including the background gauge superfield V ) should be set to 0, but the auxiliary superfield parameter g remains arbitrary. To derive the last equality, it is necessary to use eq. (4.13) and the identity This expression can be presented as a sum of certain one particle irreducible (1PI) supergraphs, because the effective action is the generating functional for 1PI Green functions (see, e.g., [90]). Therefore, it can be calculated using the tools of the perturbation theory, which include standard rules for working with supergraphs. Note that the external lines in the superdiagrams contributing to the expression (4.19) are attached to the points x, y, z 1 , and z 2 and correspond to θ 2θȧ v B x , θ 2θḃ v B y , 1, and 1, respectively.

JHEP10(2019)011
Evidently, any two points of an 1PI graph can be connected by a chain of vertices and propagators. This allows to shift v B in an arbitrary point of the supergraph, because additional terms produced by such shifts are suppressed by powers of 1/R. Really, propagators contain derivatives with respect to the superspace coordinates acting on δ 8 xy . Certainly, v B commutes with ∂/∂θ a and ∂/∂θȧ due to the independence of θ. As for the derivatives with respect to the space-time coordinates x µ , the shifting of v B from the superspace point 1 to the point 2 is made according to the procedure where we took into account that the space-time derivatives of v B are proportional to powers of 1/R, see, e.g., eq. (3.30). (To be exact, the dimensionless parameter in this case is 1/(ΛR).) Certainly, the terms proportional to 1/R can be omitted in the limit R → ∞, which is actually equivalent to the limit p → 0 in equations like eq. (3.14). Below we will always ignore them.
With the help of equations like (4.20) we can shift v B to an arbitrary point of the supergraph. Let us shift both v B in eq. (4.19) to the point z 1 , Note that in this case the usual coordinates x µ on which v B depends should be replaced by the chiral coordinates y µ = x µ + iθȧ(γ µ )ȧ b θ b to obtain a manifestly supersymmetric expression. Certainly, this is possible, because the difference is proportional to powers of 1/R and vanishes in the limit R → ∞. Also it is possible to prove thatθȧ andθ˙b in eq. (4.19) can be shifted in an arbitrary point. Really, let us consider a supergraph contributing to the expression (4.19). It is calculated according to the well-known algorithm (see, e.g., [65]), the result being given by an integral over the full superspace. 12 The integral over the full superspace includes integration over d 4 θ and does not vanish only if the integrand contains θ 4 = θ 2θ2 . Note that new θ-s cannot be produced in calculating the supergraphs, in spite of their presence inside the supersymmetric covariant derivatives. Therefore, any supergraph with θ-s on external lines does not vanish only if it contains at least two right components θ a and two left componentsθȧ. The expression (4.19) is quadratic inθ, which can be shifted along a pass consisting of vertices and propagators using equations like Here O(1) denotes terms which do not containθ. They appear when the covariant derivatives are commuted withθ-s with the help of the identity {θȧ,D˙b} = δȧ b . The arrow in eq. (4.22) points that we omit them, because these terms do not contribute to eq. (4.19). Really, the original expression is quadratic inθ, so that the contributions of O(1) terms are no more than linear inθ-s. This implies that they are removed by the final integration over d 4 θ.

JHEP10(2019)011
Thus, we see thatθ-s in supergraphs contributing to eq. (4.19) can be shifted in an arbitrary way using equations like (4.22). This allows shiftingθȧ andθ˙b from the points x and y to the point z 2 , θȧ After this, we use the identity (4.24) (Here we essentially use that bothθ-s are placed into a single point z 2 .) As a result, we obtain that after the shifts (4.21) and (4.23) the considered expression is written as . (4.25) Note that due to the antichirality ofθ 2 this expression remains manifestly supersymmetric. The right components θ cannot be shifted in an arbitrary way, because the considered expression is quartic in θ a (here we count only the degree of the right components). However, in this case it is possible to use a special identity derived in ref. [55]. Let us consider an 1PI supergraph contributing to the expression (4.25) and construct two passes connecting the point x with z 1 and the point z 1 with y, see figure 1. The corresponding sequences of vertices and propagators we will denote by A and B, respectively. Actually, A and B are products of the expressions in which various derivatives (namely, ∂ µ , D a ,Dȧ, and 1/∂ 2 ) act on superspace δ-functions. Then according to ref. [55] θ 2 ABθ 2 + 2(−1) P A +P B θ a Aθ 2 Bθ a − θ 2 Aθ 2 B − Aθ 2 Bθ 2 = O(θ), (4.26) where (−1) P X is the Grassmannian parity of an expression X, and O(θ) denotes terms which are no more than linear in θ. For completeness, we also present the proof of this identity in appendix A. (The point x is on the left of each term, the point y is on the right, and the point z 1 is between A and B.) Evidently, the O(θ) terms in eq. (4.26) do not contribute to eq. (4.25), because the integral over d 4 θ which remains after the calculation of the supergraph removes them. Therefore, with the help of eq. (4.26) the left hand side of eq. (3.40) can be rewritten in the form where we take into account that all propagators are Grassmannian even. This expression can be equivalently expressed in terms of the operatorÔ A as

JHEP10(2019)011
x y z 1 A B d d © Figure 1. The points x, z 1 , and y of a supergraph can be connected by a pass which consists of the gauge, matter, and ghost propagators. A corresponds to its part connecting the points x and z 1 , and B corresponds to the part connecting the points z 1 and y.
To see this, it is necessary to use the identity where we also took the identity into account. Eq. (4.30) is a convenient starting point for presenting the left hand side of eq. (3.40) in the form of an integral of double total derivatives. This will be made in the next section.

Formal calculation
Numerous explicit calculations of the β-function reveal that it is given by integrals of double total derivatives in the momentum space for both the Abelian [44,50] and non-Abelian [47-49, 51, 53] N = 1 supersymmetric theories regularized by higher covariant derivatives. In the Abelian case this factorization into integrals of double total derivatives has been proved in all orders in refs. [54,55]. For generalizing this result to the non-Abelian case we consider the left hand side of eq. (where ρ = gg * ) and present it in the form (4.30). Below we will demonstrate that it is given by integrals of double total derivatives in the momentum space in all orders.

JHEP10(2019)011
An important observation is that the expression (4.30) formally vanishes as a consequence of the Slavnov-Taylor identity (4.7). In fact, it is not true because of singular contributions, which will be discussed in section 4.5. However, first, we describe the formal calculation.
As a starting point we consider the Slavnov-Taylor identity (4.7) in which we set the superfields V , φ i , c A , andc A to 0. However, the auxiliary superfields remain arbitrary. This gives the equation Its left hand side is a functional of the background gauge superfield V and the auxiliary external superfields g and g * . Next, we differentiate eq. (4.33) with respect to V B y and, after this, set the background gauge superfield to 0. Then using eq. (4.5) we obtain where we also took into account that (even for g = 0) .

1) and (4.3) in the form
where ε aB is a coordinate independent anticommuting parameter. This implies that A B = ε aB θ a . Substituting these parameters into eq. (4.34) and differentiating with respect tō εȧ B , we obtain the equation the left hand side of which being a functional of the auxiliary superfield g. Therefore, it is possible to differentiate with respect to g and g * , so that the part of eq. (4.30) obtained

JHEP10(2019)011
from the second term in the round brackets vanishes. The part obtained from the first term vanishes due to the same reason. This implies that LHS of eq. (3.40) = 2 d 8 x d 8 y d 6 z 1 d 6z The similar arguments can be used for this expression (which corresponds to the third term in the round brackets in eq. (4.30)). In this case it is necessary to choose the superfield A as where a B µ are real coordinate-independent parameters. Therefore, A B = ia B µ y µ , where the chiral coordinates y µ and the antichiral coordinates (y µ ) * are defined as respectively. In this case from eq. (4.34) for arbitrary g we formally 13 obtain the identity Consequently, the expression (4.39) seems to vanish. This implies (see eqs. (3.19) and (4.32)) that all higher order corrections to the β-function vanish and the β-function is completely defined by the one-loop approximation. Certainly, it is not true. The matter is that the above calculation was made formally and something very important was missed. The origin of the incorrect result can be found analyzing the explicit calculations made with the higher covariant derivative regularization [45][46][47][50][51][52][53]. They demonstrate that all integrals giving the β-function are integrals of double total derivatives in the momentum space, and that all loop corrections come from δ-singularities. Below in section 4.4 we will see that the integrals of (double) total derivatives appear due to the presence of x µ in eq. (4.40). These total derivatives produce singular contributions which were ignored in the formal calculation. Note that eq. (4.37) does not contain x µ , so that the momentum total derivatives do not appear in the first two terms of eq. (4.30). This implies that the higher (L ≥ 2) loop corrections to the β-function are completely determined by the third term inside the round brackets in eq. (4.30). It is this term that produces the double total derivatives in the momentum space. To derive this fact in section 4.4, here we relate this term with the second variation of the functional integral giving the effective action under the change of variables corresponding to the background gauge transformations.
Let us set all quantum superfields to 0. Then the effective action will depend only on the external superfields V and g. Taking into account that (at least, in the perturbation theory) the vanishing of the quantum (super)fields corresponds to the vanishing of the sources, we obtain Γ quantum fields=0 where Z is given by the functional integral (2.23). 13 This identity is not actually valid, because the parameter A too rapidly grows at infinity.

JHEP10(2019)011
Similarly to the derivation of the Slavnov-Taylor identity in section 4.1, we perform the change of variables (4.1) in this functional integral, but the parameter A will be chosen in the form (4.40). Let us denote the variation of the effective action under the background gauge transformations of the quantum superfields byδ a . (This variation does not include the transformation of the background gauge superfield V .) Taking into account that the generating functional (4.43) remains the same after the considered change of variables, while the total action is invariant under the background gauge transformation, we obtain the equation similar to eq. (4.7), which is certainly a mere consequence of the Slavnov-Taylor identity. (Note that the background superfield V and the external superfield g are not so far set to 0.) Differentiating eq. (4.44) with respect to a B µ gives The derivative of the effective action with respect to V A entering this equation can be presented as the functional integral where the angular brackets are defined by eq. (4.6) and we also introduced the notation In this functional integral it is possible to perform again the change of variables (4.1) with the parameter A = ib B µ t B y µ . After this change of variables we set the background gauge superfield V to 0. As a result, we obtain the identity As usual, the subscript "fields = 0" means that the superfields V , φ i , c,c, and V are set to 0, while the chiral superfield g can take arbitrary values. The symbolδ b denotes the variation under the transformations (4.1) of the quantum superfields parameterized by A = ib A µ t A y µ , the background gauge superfield V being fixed. Let us transform the right hand side of this expression taking into account that the total action (4.2) and the Pauli-Villars actions S ϕ and S Φ (given by eqs.  (4.49) where δ b V is given by eq. (4.5). From eq. (4.49) it is possible to obtain the identities They can be derived by commuting the derivative with respect to V B y to the left, if we take into account that it commutes withδ b and use the equation which is valid because f AAC = 0. The operatorδ b in eq. (4.48) acts on the expression inside the angular brackets and on the actions S total , S ϕ , and S Φ in the exponents. Eqs. (4.49) and (4.50) allow expressing the result in terms of the derivatives with respect to the background gauge superfield. From the other side, the derivative of the angular brackets with respect to V also acts on the expression inside these brackets and on the actions in the exponents. This implies that The expression δ b V entering this equation is given by eq. (4.5). Differentiating it with respect to b B µ and setting the background gauge superfield to 0, we obtain Therefore, taking into account eq. (4.44), we see that the formal calculation gives (Note that in this expression we do not set the external superfield g to 0.) However, in what follows we will see that the first equality is not true, because doing the formal calculation we ignore singular contributions. These singular contributions will be discussed below.

JHEP10(2019)011
If we apply the operator to the left hand side of eq. (4.54) and, after this, set the auxiliary external superfield g to 0, then we obtain the expression (4.39), According to this equation all higher order corrections to the β-function vanish. Certainly, it is not true. As we have already mentioned above, such a result appears, because singular contributions were missed in the formal calculation described above.
Although from eq. (4.56) we obtain the same (incorrect) formal result as from eq. (4.42), eq. (4.56) will be very useful below, because it allows explaining the factorization of the loop integrals giving the β-function into integrals of double total derivatives.

Integrals of double total derivatives
Although the calculation described in the previous section is formal, it allows explaining why the β-function (defined in terms of the bare couplings with the higher derivative regularization) is given by integral of double total derivatives in the momentum space. This can be done starting from eq. (4.56). Its left hand side is related to the β-function by eq. (4.32). In this section we present the right hand side of eq. (4.56) as a sum of integrals of double total derivatives and formulate a prescription for constructing these integrals.
Let ϕ I denotes the whole set of superfields of the theory, where the index I corresponds to quantum numbers with respect to the gauge group, and j I are the corresponding sources. In the momentum representation the propagators can be presented in the form (4.57) where Z 0 is the generating functional for the free theory.
Let us make the change of the integration variables (4.1) with the parameter A given by eq. (4.40) in the generating functional Z with the sources and the background gauge superfield set to 0. Although under this change of variables the generating functional remains invariant, the propagators and vertices transform nontrivially. Really, if S 2 and S int are the quadratic part of the action and the interaction, respectively, then , where Z 0 ≡ Dϕ exp iS 2 ϕ (ϕ) +iϕ·j , (4.59) are also different from the old ones. Now, let us try to understand how the evident equality Z = Z appears at the level of superdiagrams. For this purpose we write the transformation (4.1) with the parameter (4.40) and concentrate on the terms linear in x µ , where (T A ) I J are the generators of the gauge group in a relevant representation, and the terms which do not explicitly depend on x µ are denoted by dots. 14 Then the propagator changes as Using this equation it is possible to demonstrate that in the momentum representation the change of the propagator (4.61) is related to its derivative with respect to the momentum, Next, let us proceed to the interaction vertices. An n-point vertex can be formally written in the form In (x 1 , x 2 , . . . , x n ; θ 1 , θ 2 , . . . , θ n ) 14 Note that if the sources are not set to 0, then . . In this case the arguments of the effective action change as ϕI = δW/δj I → ϕ I = ϕI + ia A µ x µ (T A )I J ϕJ + . . . This implies that the considered change of the integration variables actually generates the transformationδa.
. . , k n ; θ 1 , θ 2 , . . . , θ n ) + . . . (4.72) Next, it is necessary to note a resemblance between eq. (4.69) and eq. (4.65). In eq. (4.65) each generator actually corresponds to a propagator coming from the considered vertex exactly as momenta in eq. (4.69). This implies that such equations appear in pairs. Say, if the considered vertex is placed inside a certain graph in which the momentum k µ 2 can be expressed in terms of k µ 3 , . . . , k µ n , then where c 3 , . . . c n are some numerical coefficients. In this caseδ aV I 1 I 2 ...In will be proportional to Thus, the variationsδ a of vertices inside a supergraph contain only derivatives with respect to independent momenta. It is well known that due to the momentum conservation in each vertex (encoded in equations like eq. (4.69)) in an L loop graph without external lines only L momenta are independent. (In our case this is also true, because the momenta of all external lines vanish.) Therefore, we can mark L propagators whose momenta are considered as independent parameters, see figure 2 (which corresponds to the case L = 3). Then, using the resemblance between eq. (4.69) and eq. (4.65), it is possible to construct L independent structures in which the generators correspond to certain propagators, e.g., to the propagators whose momenta we consider as independent parameters. Any graph in which T A stands on a certain propagator can be expressed in terms of these structures.
Let us consider a closed loop, consisting of vertices and propagators, which includes one of the independent momenta, say, k µ . Then according to eqs. (4.63), (4.72) and (4.75), from the terms containing the derivative ∂/∂k µ we obtain the contribution to the first variation of the considered supergraph given by an integral of a total derivative  where the generator T A should be inserted on the propagator with the momentum k µ . This is graphically illustrated in figure 2.
The second variation is calculated similarly. Thus, we have a prescription, how to find integrals of double total derivatives which contribute to the β-function. The starting point is the expression (4.77) First, we consider a certain L loop supergraph contributing to it and (in an arbitrary way) mark L propagators with the (Euclidean) momenta Q µ i considered as independent. Let a i be the indices corresponding to their beginnings. Next, it is necessary to calculate the supergraph using the standard rules. The result includes a coefficient which contains couplings and some group factors. This coefficient should be replaced by a certain differential operator which is obtained by calculating the "second variation" of the expression i δ b i a i , where δ b i a i comes from the marked propagators, formally setting In other words, we make the replacement Next, one should multiply the result by the factor where the sign "−" appears, because

JHEP10(2019)011
Finally, it is necessary to rewrite the result in terms of ρ = gg * and perform the integration dρ. (4.82) The expression obtained according to the algorithm described above coincides with a contribution to β/α 2 0 coming from the sum of all superdiagrams which are obtained from the original vacuum supergraphs by attaching two external lines of the background gauge superfield in all possible ways.
Below in section 5 we will verify this algorithm for some particular examples.

The role of singularities
From the discussion of the previous section we can conclude that in the case of using the higher derivative regularization the integrals giving the β-function are integrals of double total derivatives. This agrees with the results of explicit calculations which also reveal that all higher order corrections to the β-function originate from singularities of the momentum integrals. Actually it is the contributions of the singularities that have been missed in the formal calculation of section 4.3. Let us demonstrate, how they appear, by considering the integral as a simple example. In eq. (4.83) Q µ denotes the Euclidean momentum, and f (Q 2 ) is a nonsingular function which rapidly tends to 0 in the limit Q 2 → ∞.
If we calculate the integral (4.83) formally, then it vanishes, because it is an integral of a total derivative. Actually, using the divergence theorem, we reduce the integral under consideration to the integral over the infinitely large sphere S 3 ∞ in the momentum space. Evidently, the result is equal to 0, because the function f vanishes on this sphere, where dS µ is the integration measure on S 3 ∞ . Actually, in section 4.3 we made a similar calculation. However, the result obtained in eq. (4.84) is evidently incorrect due to a singularity of the integrand at Q µ = 0.
To correct the above calculation, it is necessary to surround the singularity by a sphere S 3 ε of an infinitely small radius ε (with the inward-pointing normal) and take into account the integral over this sphere,

JHEP10(2019)011
Let us visualize this result by reobtaining it in a different way. First, we note that defining the integral I we actually do not distinguish between the expression (4.83) and the integral However, it is possible to introduce the operator ∂/∂Q µ which is similar to ∂/∂Q µ , but, by definition, the integral of it is always reduced to the integral over the sphere S 3 ∞ only. Moreover, we assume that this operator is commuted with Q µ /Q 4 in the integrand with the help of the identity (4.87) In terms of the operator ∂/∂Q µ the considered integral is defined as Then, if we integrate by parts taking into account vanishing of the integral of a total derivative and eq. (4.87), we obtain From this equation we see that the integral I is determined by a contribution of the δ-singularity. Note that in the coordinate representation where a is a certain function, while Such a structure of loop integrals appears in the Abelian case (see, e.g., [54]). In the non-Abelian case the structure analogous to (4.90) is the right hand side of eq. (4.56), while its left hand side is an analog of the expression (4.91). Therefore, it becomes clear that making the calculations formally in the previous section we ignored the δ-singularities. Thus, to make the calculation properly, it is necessary to take into account singular contributions, which generate all terms containing the anomalous dimensions in the NSVZ equation (1.2) for RGFs defined in terms of the bare couplings. We hope to describe how to sum these singularities in a future publications. Figure 3. Graphs generating terms containing the Yukawa couplings in the three-loop β-function.

JHEP10(2019)011
We point out independent momenta and indices corresponding to beginnings of the respective propagators using the same notations as in the calculation described in the text.

Verification in the lowest orders
To confirm the correctness of the general arguments presented above, it is desirable to verify them by explicit calculations in the lowest orders. In section 4.4 we have formulated the prescription, how to construct integrals of double total derivatives which appear in calculating the β-function in the case of using the higher covariant derivative regularization. For obtaining these integrals one usually calculates a set of superdiagrams which are obtained from a given graph by attaching two external lines of the background gauge superfield in all possible ways. For example, in ref. [51] this has been done for the three-loop contributions quartic in the Yukawa couplings. All three-loop terms containing the Yukawa couplings have been subsequently found in ref. [53]. (Both these calculations were made in the Feynman gauge ξ = 1 for the higher derivative regulator K = R.) Unfortunately, at present no other three-loop contributions to the β-function are known in the case of using the higher covariant derivative regularization. Nevertheless, the results of refs. [51,53] allow verifying the general argumentation of the present paper by comparing the algorithm described in section 4.4 with the result of the standard calculation.
A part of the three-loop β-function which contains the Yukawa couplings originates from the supergraphs presented in figure 3. Within the standard technique used in refs. [51,53] they generate large sets of superdiagrams with two external lines corresponding to the background gauge superfield which have to be calculated. However, now it is possible to derive the result for their sums by a different (and much simpler) way. Namely, we

JHEP10(2019)011
should calculate the (specially modified) superdiagrams without external lines and, after this, follow the algorithm described in section 4.4. Here we describe this calculation for the graph (1) in details and present the similar results for the remaining graphs (2) - (5).
As a starting point we find the contribution of the graph (1) to the expression (4.77). Due to the derivatives with respect to the superfields g and g * and subsequent integrations, two vertices in this graph take the form Then, after some standard calculations, for the contribution of the supergraph (1) (in the Euclidean space after the Wick rotation) we obtain Note that although here the superfield g is set to 0, the coordinate independent parameter g can in general be present in the Yukawa vertices and gauge propagators. However, the graph (1) appears to be independent on g and, therefore, on ρ = gg * . According to the prescription described in section 4.4 for obtaining the contribution to the β-function, at the first step, it is necessary to replace the factor λ ijk 0 λ * 0ijk (which in the original graph comes from the expression λ ijk 0 λ * 0pmn δ p i δ m j δ n k ) by a certain differential operator acting on the integrand in eq. (5.2). To construct this operator, we consider the propagators with the independent momenta K µ and Q µ . Let they are proportional to δ m j and δ n k , respectively. Then, we construct the second "variation" formally replacing This operation changes the Yukawa coupling dependent factor in eq. (5.2) as Replacing the factor λ ijk 0 λ * 0ijk in eq. (5.2) by this operator and taking into account that the Euclidean momenta K µ and Q µ enter the integrand of eq. (5.2) symmetrically, we obtain the expression

JHEP10(2019)011
To simplify it, we use two identities. The first one, follows from eq. (2.4), while the second one, can be verified by direct differentiating after some changes of integration variables in the resulting integrals. Then the expression under consideration takes the form To find the contribution to the function β(α 0 , λ 0 λ * 0 , Y 0 )/α 2 0 , it is necessary to multiply this expression by −2π/rV 4 and apply the operator to the result. For the graph (1) this integration gives the factor 1, because the expression for this graph does not depend on ρ. Therefore, This result exactly coincides with the one derived in ref. [51] by direct summation of the superdiagrams contributing to the two-point Green function of the background gauge superfield. Certainly, the calculation described here is much simpler, because we had to calculate the only superdiagram without external lines. The agreement of the results confirms the correctness of the general arguments presented in this paper. However, it is desirable to verify also the three-loop results corresponding to the graphs (2) -(5) in figure 3. As in refs. [51,53] we will use the Feynman gauge, so that in what follows the parameter ξ 0 is set to 1 and the higher derivative regulator K is chosen equal to R. Calculating the supergraph (2) in figure 3 we should take into account that θ 2 andθ 2 can appear in different points. This produces a set of subgraphs presented in the curly brackets in figure 4. However, all these subgraphs differ only in the numeric coefficients. Really, they are quartic in θ-s, so that these θ-s can be shifted to an arbitrary point of the supergraph. (Terms with lower degrees of θ, which can appear after such shifts, evidently vanish due to the integration over d 4 θ.) For example, it is possible to shift θ-s as it is shown in the right hand side of figure 4. 15 JHEP10(2019)011 Figure 4. Subgraphs of the supergraph (2) correspond to different positions of θ 2 andθ 2 . However, the sum of them is effectively reduced to a single supergraph in which θ 4 can be placed in an arbitrary point and g = g * = 1.
The result for their sum (in the Euclidean space after the Wick rotation) can be written as where, following ref.
[53], we use the notation As earlier, we should replace the factor λ ijk 0 λ * 0imn (T B ) j m (T B ) k n by a relevant differential operator. For constructing this differential operator we again mark the propagators with the independent momenta Q µ , L µ , and K µ , see figure 3. The beginnings of the lines which denote them correspond to the indices m, i, and B. They refer to the representations R (in which the matter superfields lie),R, and Adj, respectively. Then, the calculation of the first "variation" gives where we take into account that T Ā R = − T A t (with T A being the generators of the representation R) and T A Adj BC = −if ABC . The second "variation" is calculated in a JHEP10(2019)011 similar way. After some (rather non-trivial) transformations involving eq. (2.4) we obtain that the differential operator for the considered graph has the form Then it is necessary to repeat the same algorithm as for the graph (1), namely, n by the operator (5.14); 2. multiply the result by −2π/rV 4 ; 3. apply the operator (5.9).
The three-loop supergraphs are proportional to gg * = ρ, so that in the considered case the integration gives 16 Thus, the contribution of the graph (2) to the function β/α 2 0 takes the form We see that this result coincides with the one obtained in ref. [53] by the straightforward calculation of superdiagrams with two external legs of the background gauge superfield. The expression for the next graph (3) has the form 17) 16 In general, an L-loop supergraph is proportional to ρ L−2 , and the integration gives the factor (L − 1) −2 .
This implies that in the general case to find a contribution to eq. (4.77), it is possible to start with a vacuum supergraph contributing to the effective action with g = g * = 1 and simply insert θ 4 to an arbitrary point which contains integration over the full superspace. (Note that the integrations over d 6 x or d 6x in the Yukawa terms can always be converted to the integrals over the full superspace.)

JHEP10(2019)011
where Similar to the previous supergraphs, we replace the factor λ ijk 0 λ * 0ijl T B k m T B m l by a differential operator. To obtain this operator, we begin with calculating the first "variation" of the considered factor, The second "variation" is constructed by a similar procedure. The result can be written in the form (5.20) Proceeding according to the above described algorithm, we find the contribution of the supergraph (3) to the function β/α 2 0 , Note that the last term in eq. (5.20) is not essential, because the corresponding contribution to β/α 2 0 vanishes. (It changes the sign under the sequence of the variable changes L µ → L µ − Q µ ; Q µ → −Q µ ; K µ → −K µ .) The result (5.21) also coincides with the one obtained in ref. [53].
The expression for the supergraph (4) is Here we use the same notation as in ref. [53],

JHEP10(2019)011
where the prime and the subscript Q denote the derivative with respect to Q 2 /Λ 2 . The corresponding operator is exactly the same as for the supergraph (3) and is given by eq. (5.20). Similarly to the case of the supergraph (3), the last term in this expression does not contribute to β/α 2 0 , so that This result also agrees with the calculation of ref. [53].
The last supergraph (5) is given by the expression The first "variation" of the factor λ ijk 0 λ * 0ijl λ mnl 0 λ * 0mnk is written as The second "variation" can be found by a similar method, but, to simplify the resulting expression, it is necessary to involve the identities which follow from eq. (2.4). Using these identities and taking into account that the integrand of eq. (5.25) is symmetric in Q and L, we find the required replacement

JHEP10(2019)011
Constructing the contribution of the graph (5) to the function β/α 2 0 with the help of this operator and using the equations we obtain This expression also agrees with refs. [51,53]. Thus, we see that the algorithm described in this paper allows reproducing all results obtained earlier by the direct summation of the superdiagrams with two external lines of the background gauge superfield. Certainly, this fact can be viewed as an evidence in favour of the correctness of the general consideration made in this paper.

Conclusion
We have proved that for N = 1 supersymmetric gauge theories the integrals giving the βfunction defined in terms of the bare couplings are integrals of double total derivatives with respect to the loop momenta in all orders in the case of using the regularization by higher covariant derivatives. This fact agrees with the results of numerous explicit calculations in the lowest orders and generalizes the similar statement for the Abelian case [54,55]. The proof of the factorization into double total derivatives is a very important step towards the all-loop perturbative derivation of the exact NSVZ β-function. This derivation consists of the following main steps: 1. Using the finiteness of the triple ghost-gauge vertices (which has been demonstrated in ref. [14]) we rewrite the NSVZ equation in the equivalent form (1.2).
2. The β-function defined in terms of the bare couplings is extracted from the difference between the effective action and the classical action by the formal substitution (3.29). Then, using the identity (4.26) and the background gauge invariance, the result is presented as an integral of a double total derivative in the momentum space. This integral is reduced to the sum of singular contributions which are given by integrals of the momentum δ-functions. (This has been done in this paper.)

JHEP10(2019)011
3. The remaining step is to sum the singular contributions and to prove that they produce the anomalous dimensions of the quantum superfields in eq. (1.2). Now this work is in progress.
As a result, we presumably obtain eqs. (1.1) and (1.2) for RGFs defined in terms of the bare couplings in the case of using the higher covariant derivative regularization (in agreement with the results of explicit multiloop calculations). Due to scheme independence of these RGFs (for a fixed regularization) this statement is valid for all renormalization prescriptions.
If the NSVZ relation is really valid for RGFs defined in terms of the bare couplings for theories regularized by higher covariant derivatives, then the all-order prescription for constructing the NSVZ scheme for RGFs defined in terms of the renormalized couplings is HD+MSL. This means using of the higher covariant derivative regularization supplemented by minimal subtractions of logarithms, when only powers of ln Λ/µ are included into renormalization constants.
As a by-product of the proof presented in this paper we have obtained a simple method for constructing the loop integrals contributing to the β-function defined in terms of the bare couplings. Actually, it is necessary to calculate (a specially modified) supergraphs without external lines and replace the products of couplings and group factors by a certain differential operator specially constructed for each supergraph. The result is equal to the sum of a large number of superdiagrams which are obtained from the original supergraph by attaching two external lines of the background gauge superfield in all possible ways. Certainly, this drastically simplifies the calculations.
As an illustration of this method we considered all three-loop contributions containing the Yukawa couplings and compared the result with the one found by the standard calculation in refs. [51,53]. The coincidence of the expressions obtained by both these methods confirms the correctness of the algorithm proposed in this paper.
Anticommuting θ a with supersymmetric covariant derivatives inside A and B we obtain expressions which do not explicitly depend on θ. This implies that the right hand side of eq. (A.1) is proportional to the second degree of (explicitly written) θ. After commuting the remaining θ 2 to the left, the expression (A.1) can be presented as Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.