1 Introduction

Ultraviolet divergences in supersymmetric theories are resticted by some non-renormalization theorems. In particular, in \({{\mathcal {N}}}=1\) supersymmetric theories the superpotential does not receive divergent quantum corrections [1]. However, even in the case of \({{\mathcal {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,7,8], the \(\beta \)-function of \({{\mathcal {N}}}=1\) supersymmetric gauge theories is related to the anomalous dimension of the matter superfields \((\gamma _\phi )_j{}^i\) by an equation which is usually called “the exact NSVZ \(\beta \)-function”. For the theory with a simple gauge group G and chiral matter superfields in its representation R this equation is written as

$$\begin{aligned} \frac{\beta (\alpha ,\lambda )}{\alpha ^2} = - \frac{3 C_2 - T(R) + C(R)_i{}^j \big (\gamma _\phi \big )_j{}^i(\alpha ,\lambda )/r}{2\pi (1- C_2\alpha /2\pi )}. \end{aligned}$$
(1)

Note that we do not so far specify the definitions of the renormalization group functions (RGFs) and use the notations \(r\equiv \dim \, G\),

$$\begin{aligned}&\text {tr}(T^A T^B) \equiv T(R)\delta ^{AB};\qquad (T^A T^A)_i{}^j \equiv C(R)_i{}^j; \nonumber \\&f^{ACD} f^{BCD} \equiv C_2\delta ^{AB}, \end{aligned}$$
(2)

where \(T^A\) are the generators of the group G in the representation R such that \([T^A, T^B] = i f^{ABC} T^C\). Also we assume that the generators of the fundamental representation \(t^A\) are normalized by the condition \(\text {tr}(t^A t^B) = \delta ^{AB}/2\).

The equation (1) implies that the all-order \(\beta \)-function of the \({{\mathcal {N}}}=1\) supersymmetric Yang–Mills (SYM) theory without matter superfields is given by the geometric series. Moreover, if \({{\mathcal {N}}}=2\) supersymmetric gauge theories are considered as a special case of the \({{\mathcal {N}}}=1\) supersymmetric theories, then Eq. (1) leads to the finiteness beyond the one-loop approximation provided the quantization is made in \({{\mathcal {N}}}=2\) supersymmetric way [9, 10]. This implies that the NSVZ \(\beta \)-function is closely related to the \({{\mathcal {N}}}=2\) non-renormalization theorem derived in [11,12,13]. However, for its derivation \({{\mathcal {N}}}=2\) supersymmetry should be manifest at all steps of calculating quantum corrections. This can be achieved with the help of the harmonic superspace [14, 15] and the invariant regularization [16]. The finiteness of \({{\mathcal {N}}}=4\) SYM theory [11, 12, 17, 18] follows from the \({{\mathcal {N}}}=2\) non-renormalization theorem and, therefore, from the NSVZ \(\beta \)-function.

Equations analogous to NSVZ are also known for the Adler D-function in \({{\mathcal {N}}}=1\) SQCD [19, 20] and for the renormalization of the gaugino mass in gauge theories with softly broken supersymmetry [21,22,23]. Recently an NSVZ-like equation was constructed for the renormalization of the Fayet–Iliopoulos term in \(D=2\), \({{\mathcal {N}}}=(2,0)\) supersymmetric theories [24].

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 [25] supplemented by modified minimal subtraction [26] (this scheme is usually called \(\overline{\text {DR}}\)) Eq. (1) is not valid starting from the order \(O(\alpha ^2,\alpha \lambda ^2,\lambda ^4)\), which corresponds to the three-loop \(\beta \)-function and the two-loop anomalous dimension [27,28,29,30,31].Footnote 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 scheme-independent consequences [33, 34]. This implies that, although the calculations made in the \(\overline{\text {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 [34, 35], it is possible to construct an infinite set of the NSVZ schemes [36].

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 [37, 38] for regularizing supersymmetric theories. Unlike dimensional reduction [39], this regularization is mathematically consistent and can be formulated in a manifestly supersymmetric way in \({{\mathcal {N}}}=1\) superspace [40, 41]. 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 \({{\mathcal {N}}}=1\) higher derivative supersymmetric theories were made in [42,43,44,45,46,47,48]. 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., [49,50,51,52]) 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 [53, 54]. Similar proofs of the NSVZ-like equations have also been constructed for the Adler D-function in \({{\mathcal {N}}}=1\) SQCD [19, 20] and for the renormalization of the photino mass in \({{\mathcal {N}}}=1\) SQED with softly broken supersymmetry [55].) RGFs defined in terms of the bare couplings are scheme-independent for a fixed regularization [56], 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 [32, 57], 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 [56], when the divergences of the theory regularized by Higher covariant Derivatives are removed by Minimal Subtractions of Logarithms [58, 59]. This implies that the renormalization constants are constructed in such a way that they include only powers of \(\ln \Lambda /\mu \), where \(\Lambda \) is a regularization parameter, analogous to the ultraviolet cut-off, and \(\mu \) 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 [56]. Note that another all-loop NSVZ scheme in the Abelian case is the on-shell scheme [60]. For the renormalization of the photino mass in softly broken \({{\mathcal {N}}}=1\) SQED and for the Adler D-function in \({{\mathcal {N}}}=1\) SQCD the NSVZ-like schemes are also given by the HD+MSL prescription, see Refs. [61] and [62], 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 [50,51,52], but, nevertheless, this statement has not yet been proved in all orders. This proof (started in Refs. [2, 63, 64]) 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 \(\beta \)-function are integrals of double total derivatives with respect to loop momenta.Footnote 2 Certainly, at least one of the loop integrals can be calculated analytically using equations like

$$\begin{aligned} \int \frac{d^4Q}{(2\pi )^4} \frac{\partial ^2}{\partial Q_\mu \partial Q^\mu }\bigg (\frac{f(Q^2)}{Q^2}\bigg ) = \frac{1}{4\pi ^2} f(0), \end{aligned}$$
(3)

where \(f(Q^2)\) is a nonsingular function of the Euclidean momentum \(Q_\mu \) 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_\varepsilon \) surrounding the point \(Q_\mu =0\). If we consider an L-loop contribution to the \(\beta \)-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 [53, 54], 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 [68, 69]

$$\begin{aligned} \frac{\beta (\alpha _0)}{\alpha _0^2} = \frac{N_f}{\pi }\bigg (1-\gamma (\alpha _0)\bigg ), \end{aligned}$$
(4)

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 \(\beta \)-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

$$\begin{aligned} \frac{\beta (\alpha ,\lambda )}{\alpha ^2}= & {} - \frac{1}{2\pi }\bigg (3 C_2 - T(R) - 2C_2 \gamma _c(\alpha ,\lambda ) - 2C_2 \gamma _V(\alpha ,\lambda ) \nonumber \\&\quad + \frac{1}{r} C(R)_i{}^j \big (\gamma _\phi \big )_j{}^i(\alpha ,\lambda )\bigg ), \end{aligned}$$
(5)

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 \(\gamma _V\), \(\gamma _c\), and \((\gamma _\phi )_j{}^i\), respectively. This fact has a simple graphical interpretation analogous to the \({{\mathcal {N}}}=1\) SQED case considered in [66, 70, 71]. Namely, given a supergraph without external lines, the corresponding superdiagrams contributing to the \(\beta \)-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 \(\beta \)-function defined in terms of the bare couplings. According to [63] 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., [50, 51, 72,73,74]. 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.Footnote 3 In the lowest orders this was explicitly demonstrated in [50,51,52, 58, 74]. In all loops the sums of singularities obtained by cutting the matter and Faddeev–Popov ghost propagators have been found in [64] using the method based on the Schwinger–Dyson effective superdiagrams proposed in [77]. These sums coincide with the terms containing \((\gamma _\phi )_j{}^i\) and \(\gamma _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 \(\gamma _V\) in Eq. (5).Footnote 4 This is a purpose of the present paper.

The paper is organized as follows. Section 2 describes the quantization of \({{\mathcal {N}}}=1\) supersymmetric gauge theories regularized by higher covariant derivatives. In Sect. 3 we rewrite the NSVZ \(\beta \)-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 \(\beta /\alpha _0^2\) which was proposed in [63], see also [52]. This method is described in Sect. 4.1. It is slightly modified in Sect. 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 Sect. 5. In particular, in Sect. 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 \(\beta \)-function. Combining the results we derive Eq. (5) for RGFs defined in terms of the bare couplings. In the next Sect. 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 Sect. 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 \(\beta /\alpha _0^2\) proposed in Sect. 4.2.

2 \({{\mathcal {N}}}=1\) renormalizable supersymmetric gauge theories regularized by higher covariant derivatives

We will consider a general renormalizable \({{\mathcal {N}}}=1\) SYM theory with a simple gauge group G and chiral matter superfields \(\phi _i\) in a representation R. In the massless limit in terms of \({{\mathcal {N}}}=1\) superfields [79,80,81] this theory is described by the action

$$\begin{aligned} S= & {} \frac{1}{2e_0^2}\, \text {Re}\,\text {tr} \int d^6x\, W^a W_a \nonumber \\&+ \frac{1}{4} \int d^8x\, \phi ^{*i} (e^{2{{\mathcal {F}}}(V)} e^{2{\varvec{V}}})_i{}^j \phi _j \nonumber \\&+ \bigg (\frac{1}{6} \lambda _0^{ijk} \int d^6x\, \phi _i \phi _j \phi _k + \text {c.c.}\bigg ), \end{aligned}$$
(6)

where for the superspace integration measures we use the brief notations

$$\begin{aligned}&d^8x \equiv d^4x\, d^4\theta _x;\qquad d^6x \equiv d^4x\, d^2\theta _x; \nonumber \\&d^6{\bar{x}} \equiv d^4x\, d^2{\bar{\theta }}_x. \end{aligned}$$
(7)

In Eq. (6) \({\varvec{V}}\) is the Hermitian background gauge superfield and V is the quantum gauge superfield, which satisfies the constraint \(V^+ = e^{-2{\varvec{V}}} V e^{2{\varvec{V}}}\). 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 \(V = e_0 V^A t^A\) and \({\varvec{V}} = e_0 {\varvec{V}}^A t^A\), while in the second term \(V = e_0 V^A T^A\) and \({\varvec{V}} = e_0 {\varvec{V}}^A T^A\), 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 \(\lambda _0^{ijk}\), respectively, where the subscript 0 points the bare values. The function \({{\mathcal {F}}}(V) = e_0 {{\mathcal {F}}}(V)^A t^A\) is needed, because the quantum gauge superfield is renormalized in a non-linear way [82,83,84]. In the lowest approximation this function was found in Refs. [85, 86] and is written as

$$\begin{aligned} {{\mathcal {F}}}(V)^A = V^A + y_0\, e_0^2\, G^{ABCD}\, V^B V^C V^D + \cdots , \end{aligned}$$
(8)

where \(y_0\) is the first bare constant in an infinite set of parameters describing the nonlinearity, and \(G^{ABCD} \equiv \bigg (f^{AKL} f^{BLM} f^{CMN} f^{DNK} + \text {permutations of} \, B, C, \text {and} \, D\bigg )/6\). The necessity of introducing the function \({{\mathcal {F}}}(V)\) was also confirmed by the calculation of Ref. [87], 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

$$\begin{aligned} W_a \equiv \frac{1}{8} {\bar{D}}^2 \bigg [ e^{-2{\varvec{V}}} e^{-2{{\mathcal {F}}}(V)}\, D_a \bigg (e^{2{{\mathcal {F}}}(V)}e^{2{\varvec{V}}}\bigg )\bigg ], \end{aligned}$$
(9)

which is a right spinor with respect to the Lorentz group.

If the Yukawa couplings satisfy the equation

$$\begin{aligned} \lambda _0^{ijm} (T^A)_m{}^k + \lambda _0^{imk} (T^A)_m{}^j + \lambda _0^{mjk} (T^A)_m{}^i = 0, \end{aligned}$$
(10)

the theory (6) is invariant under the background gauge transformations

$$\begin{aligned}&e^{2{\varvec{V}}} \rightarrow e^{-A^+} e^{2{\varvec{V}}} e^{-A};\qquad V \rightarrow e^{-A^+} V e^{A^+}; \nonumber \\&\phi _i \rightarrow (e^A)_i{}^j \phi _j \end{aligned}$$
(11)

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

$$\begin{aligned} S_{\mathrm{reg}}= & {} \frac{1}{2 e_0^2}\,\text {Re}\, \text {tr} \int d^6x\, W^a \nonumber \\&\times \bigg [e^{-2{\varvec{V}}} e^{-2{{\mathcal {F}}}(V)}\, R\bigg (-\frac{{\bar{\nabla }}^2 \nabla ^2}{16\Lambda ^2}\bigg )\, e^{2{{\mathcal {F}}}(V)}e^{2{\varvec{V}}}\bigg ]_{Adj} W_a \qquad \nonumber \\&+ \frac{1}{4} \int d^8x\, \phi ^{*i} \bigg [\, F\bigg (-\frac{{\bar{\nabla }}^2 \nabla ^2}{16\Lambda ^2}\bigg ) e^{2{{\mathcal {F}}}(V)}e^{2{\varvec{V}}}\bigg ]_i{}^j \phi _j \nonumber \\&+ \bigg (\frac{1}{6} \lambda _0^{ijk} \int d^6x\, \phi _i \phi _j \phi _k + \text {c.c.} \bigg ), \end{aligned}$$
(12)

where for \(f(x) = 1 + f_1 x + f_2 x^2 +\cdots \)

$$\begin{aligned} f(X)_{Adj} Y \equiv Y + f_1 [X, Y] + f_2 [X,[X,Y]] + \cdots , \end{aligned}$$
(13)

and the explicit expressions for the covariant derivatives have the form

$$\begin{aligned} \nabla _a = D_a;\qquad {\bar{\nabla }}_{\dot{a}} = e^{2{{\mathcal {F}}}(V)} e^{2{\varvec{V}}} {\bar{D}}_{\dot{a}} e^{-2{\varvec{V}}} e^{-2{{\mathcal {F}}}(V)}. \end{aligned}$$
(14)

The modification of the action \(S \rightarrow S_{\mathrm{reg}}\) allows regularizing all divergences beyond the one-loop approximation [88].

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 one-loop divergences [89],

$$\begin{aligned} Z= & {} \int D\mu \, \text {Det}(PV,M_\varphi )^{-1}\text {Det}(PV,M)^c \nonumber \\&\exp \bigg \{i\bigg (S_{\mathrm{reg}} + S_{\mathrm{gf}} + S_{\mathrm{FP}} + S_{\mathrm{NK}} + S_{\mathrm{sources}}\bigg )\bigg \}. \end{aligned}$$
(15)

The source term is given by the expression

$$\begin{aligned} S_{\mathrm{sources}}= & {} \int d^8x\, J^A V^A \nonumber \\&+ \bigg (\int d^6x\, \bigg (j^i \phi _i + j_c^A c^A + {\bar{j}}_c^A {\bar{c}}^A\bigg ) + \text {c.c.}\bigg ),\nonumber \\ \end{aligned}$$
(16)

where \(J^A\) and \(j^i\) are the real and chiral superfields, respectively. The anticommuting chiral superfields \(j_c^A\) and \({\bar{j}}_c^A\) are the sources for the anticommuting chiral Faddeev–Popov ghost and antighost superfields denoted by \(c^A\) and \({\bar{c}}^A\), respectively.

We will use the gauge fixing term

$$\begin{aligned} S_{\mathrm{gf}} = -\frac{1}{16\xi _0 e_0^2}\, \text {tr} \int d^8x\, {\varvec{\nabla }}^2 V K\bigg (-\frac{{\varvec{{\bar{\nabla }}}}^2 {\varvec{\nabla }}^2}{16\Lambda ^2}\bigg )_{Adj} {\varvec{{\bar{\nabla }}}}^2 V\nonumber \\ \end{aligned}$$
(17)

invariant under the background gauge transformations (11) due to the presence of the background covariant derivatives

$$\begin{aligned} {\varvec{\nabla }}_a \equiv D_a;\qquad \quad {\varvec{{\bar{\nabla }}}}_{\dot{a}} \equiv e^{2{\varvec{V}}} {\bar{D}}_{\dot{a}} e^{-2{\varvec{V}}}. \end{aligned}$$
(18)

The corresponding actions for the Faddeev–Popov and Nielsen–Kallosh ghosts (\(S_{\mathrm{FP}}\) and \(S_{\mathrm{NK}}\), respectively) can be found in Refs. [63, 64]. In this paper we will not use the explicit expressions for them. The bare gauge parameter \(\xi _0\) and the parameters present in the function \({{\mathcal {F}}}(V)\) (i.e., \(y_0\), etc.) can conveniently be included into a single infinite set \(Y_0 \equiv (\xi _0,\,y_0,\ldots )\).

Following Refs. [78, 90], for regularizing the one-loop divergences we use two sets of the Pauli–Villars superfields. The first one consists of three commuting chiral superfields \(\varphi _a\) in the adjoint representation of the gauge group with the mass \(M_\varphi = a_\varphi \Lambda \). They cancel the divergences coming from the gauge and ghost loops. The corresponding determinant present in Eq. (15) is denoted by \( \text {Det}(PV,M_\varphi )\). The second set giving the determinant \(\text {Det}(PV,M)\) consists of the commuting chiral superfields \(\Phi _i\) in a certain representation \(R_{\mathrm{PV}}\) that admits the gauge invariant mass term such that \(M^{ij} M_{jk} = M^2 \delta ^i_k\) with \(M = a\Lambda \).Footnote 5 These superfields cancel the one-loop divergences produced by the matter superfields \(\phi _i\) if in Eq. (15) the degree of the corresponding determinant is \(c=T(R)/T(R_{\mathrm{PV}})\). Again we will not use explicit expressions of the Pauli–Villars determinants which can also be found in Refs. [63, 64]. It should be only mentioned that the constants \(a_\varphi \) and a must be independent of couplings.

Starting from the generating functional for the connected Green functions \(W \equiv -i\ln Z\), one can construct the effective action

$$\begin{aligned} \Gamma [{\varvec{V}}, V, \phi _i,\ldots ] \equiv W - S_{\mathrm{sources}}\bigg |_{\mathrm{sources}\, \rightarrow \, \mathrm{fields}}, \end{aligned}$$
(19)

where it is necessary to express the sources \(J^A, j^i,\ldots \) in terms of quantum superfields \(V^A,\phi _i,\ldots \) using the equations

$$\begin{aligned} V^A \equiv \frac{\delta W}{\delta J^A},\qquad \phi _i \equiv \frac{\delta W}{\delta j_i},\qquad \text {etc}. \end{aligned}$$
(20)

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 two-point Green functions of the theory. The terms in the effective action corresponding to these Green functions can be presented in the form

$$\begin{aligned} \Gamma ^{(2)}_{{\varvec{V}}}= & {} - \frac{1}{8\pi } \text {tr} \int \frac{d^4p}{(2\pi )^4}\, d^4\theta \, {\varvec{V}}(-p,\theta ) \nonumber \\&\times \partial ^2 \Pi _{1/2} {\varvec{V}}(p,\theta )\, d^{-1}(\alpha _0, \lambda _0, Y_0, \Lambda /p); \end{aligned}$$
(21)
$$\begin{aligned} \Gamma ^{(2)}_V - S_{\mathrm{gf}}^{(2)}= & {} -\frac{1}{4} \int \frac{d^4q}{(2\pi )^4}\, d^4\theta \, V^A(-q,\theta )\nonumber \\&\times \partial ^2\Pi _{1/2} V^A(q,\theta )\, G_V(\alpha _0, \lambda _0, Y_0, \Lambda /q); \end{aligned}$$
(22)
$$\begin{aligned} \Gamma ^{(2)}_\phi= & {} \frac{1}{4}\int \frac{d^4q}{(2\pi )^4}\, d^4\theta \, \phi ^{*i}(-q,\theta ) \phi _j(q,\theta ) \nonumber \\&\times \big (G_\phi \big )_i{}^j(\alpha _0, \lambda _0, Y_0, \Lambda /q); \end{aligned}$$
(23)
$$\begin{aligned} \Gamma ^{(2)}_c= & {} \frac{1}{4}\int \frac{d^4q}{(2\pi )^4}\, d^4\theta \nonumber \\&\times \bigg ({\bar{c}}^{+A}(-q,\theta ) c^{A}(q,\theta ) - {\bar{c}}^A(-q,\theta ) c^{+A}(q,\theta )\bigg ) \nonumber \\&\times G_c(\alpha _0, \lambda _0, Y_0, \Lambda /q). \end{aligned}$$
(24)

where \(\partial ^2\Pi _{1/2} \equiv - D^a {\bar{D}}^2 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 [75, 76].

In our notation the renormalization constants are introduced with the help of the equationsFootnote 6

$$\begin{aligned} Z_\alpha\equiv & {} \frac{\alpha }{\alpha _0};\qquad \quad V^A \equiv Z_V (V_R)^A;\qquad \quad Z_y \equiv \frac{y_0}{y};\qquad \quad \cdots \nonumber \\ Z_\xi\equiv & {} \frac{\xi }{\xi _0};\qquad {\bar{c}}^A c^B = Z_c ({\bar{c}}_R)^A (c_R)^B; \nonumber \\ \phi _i= & {} \left( \sqrt{Z_\phi }\right) _i{}^j \left( \phi _R\right) _j,\qquad \end{aligned}$$
(25)

where the subscript R stands for the renormalized superfields. The renormalization constants are constructed from the requirement that the functions \(d^{-1}\), \(Z_V^2 G_V\), \((Z_\phi G_\phi )_i{}^j\), and \(Z_c G_c\) expressed in terms of the renormalized couplings \(\alpha \), \(\lambda ^{ijk}\), and Y should be finite in the limit \(\Lambda \rightarrow \infty \). Note that due to the non-renormalization of the superpotential [1] the renormalized Yukawa couplings are given by

$$\begin{aligned} \lambda ^{ijk} = \left( \sqrt{Z_\phi }\right) _m{}^i \left( \sqrt{Z_\phi }\right) _n{}^j \left( \sqrt{Z_\phi }\right) _p{}^k \lambda _0^{mnp}. \end{aligned}$$
(26)

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]Footnote 7 it is possible to choose a subtraction scheme in which the renormalization constants satisfy the condition

$$\begin{aligned} Z_\alpha ^{-1/2} Z_c Z_V =1. \end{aligned}$$
(27)

Note that this equation is compatible with minimal subtractions of logarithms, because in the HD+MSL scheme all renormalization constants \(Z_\alpha \), \(Z_c\), and \(Z_V\) are equal to 1 for \(\mu =\Lambda \).

The Eq. (27) allows to present the NSVZ equation in an equivalent form relating the \(\beta \)-function to the anomalous dimensions of the quantum gauge superfield, the Faddeev–Popov ghosts, and the matter superfields,

$$\begin{aligned} \frac{\beta (\alpha _0,\lambda _0,Y_0)}{\alpha _0^2}= & {} - \frac{1}{2\pi }\bigg (3 C_2 - T(R) - 2C_2 \gamma _c(\alpha _0,\lambda _0,Y_0)\nonumber \\&\quad - 2C_2 \gamma _V(\alpha _0,\lambda _0,Y_0) \nonumber \\&\quad + \frac{1}{r} C(R)_i{}^j \big (\gamma _\phi \big )_j{}^i(\alpha _0,\lambda _0,Y_0)\bigg ), \end{aligned}$$
(28)

where we take into account that RGFs (at least, \(\gamma _V\) and \(\gamma _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

$$\begin{aligned} \frac{\beta (\alpha _0,\lambda _0,Y_0)}{\alpha _0^2}= & {} -\frac{d}{d\ln \Lambda }\bigg (\frac{1}{\alpha _0}\bigg )\bigg |_{\alpha ,\lambda ,Y=\mathrm{const}} \nonumber \\= & {} \frac{d}{d\ln \Lambda } \bigg (d^{-1} - \alpha _0^{-1}\bigg )\bigg |_{\alpha ,\lambda ,Y=\mathrm{const};\, p\rightarrow 0};\qquad \end{aligned}$$
(29)
$$\begin{aligned} \gamma _V(\alpha _0,\lambda _0,Y_0)= & {} - \frac{d\ln Z_V}{d\ln \Lambda }\bigg |_{\alpha ,\lambda ,Y = \mathrm{const}} \nonumber \\= & {} \frac{1}{2}\,\frac{d\ln G_V}{d\ln \Lambda }\bigg |_{\alpha ,\lambda ,Y = \mathrm{const};\ q\rightarrow 0}; \end{aligned}$$
(30)
$$\begin{aligned} \gamma _c(\alpha _0,\lambda _0,Y_0)= & {} - \frac{d\ln Z_c}{d\ln \Lambda }\bigg |_{\alpha ,\lambda ,Y = \mathrm{const}} \nonumber \\= & {} \frac{d\ln G_c}{d\ln \Lambda }\bigg |_{\alpha ,\lambda ,Y = \mathrm{const};\ q\rightarrow 0}; \end{aligned}$$
(31)
$$\begin{aligned} \big (\gamma _\phi \big )_i{}^j(\alpha _0,\lambda _0,Y_0)= & {} - \frac{d\big (\ln Z_\phi \big )_i{}^j}{d\ln \Lambda }\bigg |_{\alpha ,\lambda ,Y = \mathrm{const}} \nonumber \\= & {} \frac{d\big (\ln G_\phi \big )_i{}^j}{d\ln \Lambda }\bigg |_{\alpha ,\lambda ,Y = \mathrm{const};\ q\rightarrow 0}. \end{aligned}$$
(32)

That is why Eq. (28) can also be rewritten as an equation relating these Green functions

$$\begin{aligned}&\left. \frac{d}{d\ln \Lambda } \bigg (d^{-1} -\alpha _0^{-1}\bigg ) \right| _{\alpha ,\lambda ,Y = \mathrm{const};\ p\rightarrow 0} \nonumber \\&\quad = -\frac{1}{2\pi } \bigg (3C_2-T(R)\bigg ) + \frac{1}{2\pi }\, \frac{d}{d\ln \Lambda } \bigg ( 2C_2 \ln G_c \qquad \nonumber \\&\qquad + C_2 \ln G_V - \frac{1}{r} C(R)_i{}^j \big (\ln G_\phi \big )_j{}^i\bigg )\bigg |_{\alpha ,\lambda ,Y = \mathrm{const};\ q\rightarrow 0},\nonumber \\ \end{aligned}$$
(33)

which was first proposed in Ref. [2]. In this paper we will prove it in all orders of the perturbation theory.

4 Integrals of total derivatives

4.1 How to construct and calculate integrals of total derivatives

The key observation needed for deriving the NSVZ \(\beta \)-function is the factorization of loop integrals which give the function \(\beta (\alpha _0,\lambda _0,Y_0)/\alpha _0^2\) 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. [63]. 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 supergraphsFootnote 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. [52, 74, 93] 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 \(\beta \)-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

$$\begin{aligned} \frac{\beta (\alpha _0,\lambda _0,Y_0)}{\alpha _0^2} - \frac{\beta _{\mathrm{1-loop}}(\alpha _0)}{\alpha _0^2}, \end{aligned}$$
(34)

where

$$\begin{aligned} \beta _{\mathrm{1-loop}}(\alpha _0) = -\frac{\alpha _0^2}{2\pi }\left( 3C_2-T(R)\right) \end{aligned}$$
(35)

is the one-loop \(\beta \)-function defined in terms of the bare couplings. (For the higher covariant derivative regularization considered in this paper it was calculated in Ref. [78].)

Now, let us formulate the algorithm for constructing contributions to the function (34) following Refs. [52, 63].

  1. 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. 2.

    Next, one needs to find a point with the integration over \(d^4\theta \) (or convert the integration over \(d^2\theta \) into the integration over \(d^4\theta \)) and insert the expression \(\theta ^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\rightarrow \infty \). (Note that without this insertion any vacuum supergraph vanishes due to the integration over the anticommuting variables \(\theta \).)

  3. 3.

    We calculate the supergraph modified by the above insertion and omit terms suppressed by powers of \(1/(R\Lambda )\). As a result we obtain an expression proportional to

    $$\begin{aligned} {{\mathcal {V}}}_4 \equiv \int d^4x\, \left( v^B\right) ^2 \rightarrow \infty . \end{aligned}$$
    (36)
  4. 4.

    At the next step it is necessary to choose L propagators with independent momenta denoted by \(Q_i^\mu \), where \(i=1,\ldots ,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 \(\prod _{i=1}^L \delta _{a_i}^{b_i}\).

  5. 5.

    The product \(\prod _{i=1}^L \delta _{a_i}^{b_i}\) coming from the marked propagators in the integrand of the (Euclidean) loop integral should formally be replaced by the operator

    $$\begin{aligned} \sum _{k,l=1}^L \bigg (\prod \limits _{i\ne k,l} \delta _{a_i}^{b_i}\bigg )\, (T^A)_{a_k}{}^{b_k} (T^A)_{a_l}{}^{b_l} \frac{\partial ^2}{\partial Q^\mu _k \partial Q^\mu _l}. \end{aligned}$$
    (37)
  6. 6.

    At the last step, it is necessary to apply the operator

    $$\begin{aligned} - \frac{2\pi }{r{{\mathcal {V}}}_4}\cdot \frac{d}{d\ln \Lambda } \end{aligned}$$
    (38)

to the resulting expression, where the derivative with respect to \(\ln \Lambda \) 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 \({\varvec{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,

$$\begin{aligned} \frac{\partial ^2}{\partial Q_\mu \partial Q_\mu } \bigg (\frac{1}{Q^2}\bigg ) = -4\pi ^2 \delta ^4(Q). \end{aligned}$$
(39)

We will also need a modification of this identity which is obtained when the derivatives are taken with respect to the different momenta,

$$\begin{aligned}&\frac{\partial ^2}{\partial Q_{\mu ,1} \partial Q_{\mu ,2}} \bigg (\frac{1}{(a_1 Q_{\mu ,1}+ a_2 Q_{\mu ,2}+ Q_{\mu ,3})^2}\bigg ) \nonumber \\&\quad = - 4\pi ^2 a_1 a_2\, \delta ^4\bigg (a_1 Q_{\mu ,1}+ a_2 Q_{\mu ,2}+ Q_{\mu ,3}\bigg ), \end{aligned}$$
(40)

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

$$\begin{aligned}&\int \frac{d^4Q}{(2\pi )^4}\, \frac{\partial ^2}{\partial Q_\mu \partial Q_\mu }\bigg (\frac{f(Q^2)}{Q^2}\bigg ) \nonumber \\&\quad = \frac{1}{4\pi ^2} f(0) = \int \frac{d^4Q}{(2\pi )^4}\, f(Q^2)\cdot 4\pi ^2 \delta ^4(Q). \end{aligned}$$
(41)

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 three-loop calculations in Refs. [51, 74].

4.2 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)\times 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

$$\begin{aligned} \left( \begin{array}{cccc} M_{1,1} &{} M_{1,2} &{} \cdots &{} M_{1,P}\\ M_{2,1} &{} M_{2,2} &{} \cdots &{} M_{2,P}\\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ M_{V-1,1} &{} M_{V-1,2} &{} \cdots &{} M_{V-1,P} \end{array} \right) \left( \begin{array}{c} Q_{\mu ,1}\\ Q_{\mu ,2}\\ \vdots \\ Q_{\mu ,P} \end{array} \right) =0. \end{aligned}$$
(42)

The matrix with the elements \(M_{i,J}\), where \(i=1,\ldots ,V-1\) and \(J=1,\ldots ,P\), is the required matrix matched to the considered vacuum supergraph. From the above equations, which can briefly be written as

$$\begin{aligned} \sum \limits _{J=1}^P M_{i,J} Q_{\mu ,J} = 0, \end{aligned}$$
(43)

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_{\mu ,\alpha }\), \(\alpha =1,\ldots , L\) in the considered supergraph, and the others can be expressed in terms of them,

$$\begin{aligned} Q_{\mu ,I} = \sum \limits _{\alpha =1}^L N_{I,\alpha } Q_{\mu ,\alpha }. \end{aligned}$$
(44)
Fig. 1
figure 1

Examples of vacuum supergraphs, for which we construct the matrices M. The bold lines denote propagators the momenta of which are considered as independent

As a simple example we can consider a supergraph presented in Fig. 1 on the left. In this case Eq. (43), the matrix M, and Eq. (44) read as

$$\begin{aligned}&\bigg (1\ 1 \ 1\bigg ) \left( \begin{array}{c} Q_{\mu ,1}\\ Q_{\mu ,2}\\ Q_{\mu ,3} \end{array} \right) = 0;\qquad M = \bigg (1\ 1 \ 1\bigg ); \nonumber \\&\left( \begin{array}{c} Q_{\mu ,1}\\ Q_{\mu ,2}\\ Q_{\mu ,3} \end{array} \right) = \left( \begin{array}{c} 1\\ 0\\ -1 \end{array} \right) Q_{\mu ,1} + \left( \begin{array}{c} 0\\ 1\\ -1 \end{array} \right) Q_{\mu ,2}. \end{aligned}$$
(45)

For a more complicated three-loop graph presented in Fig. 1 on the right as another example, the matrix M takes the form

$$\begin{aligned} M = \left( \begin{array}{cccccc} 0 &{} -1 &{} -1 &{} 1 &{} 0 &{} 0\\ 1 &{} 0 &{} 0 &{} -1 &{} 1 &{} 0\\ 0 &{} 1 &{} 0 &{} 0 &{} -1 &{} 1 \end{array} \right) . \end{aligned}$$
(46)

Evidently, in this case there are three independent momenta. For example, it is possible to choose \(Q_{\mu ,1}\), \(Q_{\mu ,2}\), and \(Q_{\mu ,3}\) as independent variables. The corresponding propagators are denoted in Fig. 1 by the bold lines. Then the Eq. (44) takes the form

$$\begin{aligned} \left( \begin{array}{c} Q_{\mu ,1} \\ Q_{\mu ,2} \\ Q_{\mu ,3} \\ Q_{\mu ,4} \\ Q_{\mu ,5} \\ Q_{\mu ,6} \end{array} \right) = \left( \begin{array}{c} 1 \\ 0 \\ 0 \\ 0 \\ -1 \\ -1 \end{array} \right) Q_{\mu ,1} + \left( \begin{array}{c} 0 \\ 1 \\ 0 \\ 1 \\ 1 \\ 0 \end{array} \right) Q_{\mu ,2} + \left( \begin{array}{c} 0 \\ 0 \\ 1 \\ 1 \\ 1 \\ 1 \end{array} \right) Q_{\mu ,3}.\nonumber \\ \end{aligned}$$
(47)

An important observation made in Ref. [63] 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

$$\begin{aligned} \int d^8x\, {\hat{V}}^{I_1 I_2 \ldots I_n}\, \varphi _{I_1} \varphi _{I_2} \ldots \varphi _{I_n}, \end{aligned}$$
(48)

where \({\hat{V}}\) is an operator acting on the product of various superfields of the theory denoted by \(\varphi _I\). Then the energy-momentum conservation in the considered vertex is expressed by the equation

$$\begin{aligned} Q_{\mu ,1} + Q_{\mu ,2} + \cdots + Q_{\mu ,n} = 0, \end{aligned}$$
(49)

which is very similar to the equation which follows from the gauge invariance of theory

$$\begin{aligned}&(T^A)_{K}{}^{I_1}\, {\hat{V}}^{K I_2 \ldots I_n} + (T^A)_{K}{}^{I_2}\, {\hat{V}}^{I_1 K \ldots I_n} + \cdots \nonumber \\&\quad + (T^A)_{K}{}^{I_n}\, {\hat{V}}^{I_1 I_2 \ldots K} = 0, \end{aligned}$$
(50)

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

$$\begin{aligned} \left( T^A_{R}\right) _K{}^I \ \rightarrow \ - \left( T^A_{{\bar{R}}}\right) _K{}^I = \left( T^A_R\right) _K{}^I. \end{aligned}$$
(51)

Let us consider a vacuum supergraph and replace a \(\delta \)-symbol \(\delta _{a_J}^{b_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 operator \({\hat{V}}^{I_1 I_2\ldots I_n}\) in a certain vacuum supergraph by the left hand side of Eq. (50), then, using this notation, we obtain the equation

$$\begin{aligned} \sum \limits _{J=1}^P M_{i,J} \left[ \,T^A\right] _J = 0. \end{aligned}$$
(52)

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),

$$\begin{aligned} \left[ \,T^A\right] _I = \sum \limits _{\alpha =1}^L N_{I,\alpha } \left[ \,T^A\right] _\alpha , \end{aligned}$$
(53)

where \(N_{I,\alpha }\) are exactly the same coefficients as in Eq. (44).

According to the algorithm described in Sect. 4.1 we need to insert \(\theta ^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

$$\begin{aligned} \sum \limits _{\alpha ,\beta =1}^L \bigg [\,\theta ^4 (v^B)^2; \frac{\partial ^2}{\partial Q_{\mu ,\alpha } \partial Q_{\mu ,\beta }} \left[ \, T^A\right] _\alpha \left[ \, T^A\right] _\beta \bigg ]. \end{aligned}$$
(54)

This means that \(\theta ^4 (v^B)^2\) is inserted into an arbitrary point of the supergraph, two \(\delta \)-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

$$\begin{aligned}&\frac{\beta (\alpha _0,\lambda _0,Y_0)}{\alpha _0^2} - \frac{\beta _{\mathrm{1-loop}}(\alpha _0)}{\alpha _0^2} \nonumber \\&\quad = -\frac{2\pi }{r{{\mathcal {V}}}_4} \frac{d}{d\ln \Lambda } \sum \limits _{\begin{array}{c} \mathrm{{vacuum}} \\ \mathrm{{supergraphs}} \end{array}} \sum \limits _{\alpha ,\beta =1}^L \nonumber \\&\qquad \times \bigg [\,\theta ^4 (v^B)^2;\, \frac{\partial ^2}{\partial Q_{\mu ,\alpha } \partial Q_{\mu ,\beta }} \left[ \, T^A\right] _\alpha \left[ \, T^A\right] _\beta \bigg ]. \end{aligned}$$
(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_I^2\). It is important to note that expressions for separate supergraphs (presented in the form of scalar integrals) can contain not only \(1/Q_I^2\), but also \(1/(Q_I^2)^n \equiv 1/Q_I^{2n}\) with \(n\ge 2\) if \(Q_I\) is a momentum of a quantum gauge superfield propagator. Acting on \(1/Q_I^{2n}\) with \(n\ge 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_I^{2n}\) singularities should cancel each other. The calculations in the lowest orders [51, 74] exactly confirm this statement.

Because all bad terms should cancel each other, it is necessary to consider only the \(1/Q_I^2\) 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,

(56)

where \(Q_{\mu ,I} \ne Q_{\mu ,J}\) for \(I\ne 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.

Fig. 2
figure 2

The structure of vacuum supergraphs with coinciding momenta (which are equal to \(k_\mu \)). The gray circles denote 1PI subdiagrams

\(\delta \)-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,

$$\begin{aligned} \frac{\partial ^2}{\partial Q_{\mu ,\alpha } \partial Q_{\mu ,\beta }} \bigg (\frac{1}{Q_I^2}\bigg ) = - 4\pi ^2 N_{I,\alpha } N_{I,\beta }\, \delta ^4(Q_I), \end{aligned}$$
(57)

where we also took Eq. (44) into account. This implies that acting on the product (56) the double total derivatives give the singular contribution

(58)

Taking into account that the bad terms containing \(1/Q_I^{2n}\) with \(n\ge 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_I^2\) by \(4\pi ^2\delta ^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

$$\begin{aligned}&-\frac{2\pi }{r{{\mathcal {V}}}_4} \frac{d}{d\ln \Lambda } \sum \limits _{\begin{array}{c} \mathrm{{vacuum}} \\ \mathrm{{supergraphs}} \end{array}} \sum \limits _{I=1}^{P'} \sum \limits _{\alpha ,\beta =1}^L N_{I,\alpha } N_{I,\beta }\nonumber \\&\qquad \times \bigg [\,\theta ^4 (v^B)^2;\, \left[ \, T^A\right] _\alpha \left[ \, T^A\right] _\beta ;\, \frac{1}{Q^{2n}}\bigg |_{n\ge 2} \nonumber \\&\quad \rightarrow 0;\, \frac{1}{Q_I^2} \rightarrow 4\pi ^2\delta ^4(Q_I) \bigg ]. \end{aligned}$$
(59)

Next, using Eq. (53) we calculate the sums over the indices \(\alpha \) and \(\beta \),

$$\begin{aligned}&-\frac{2\pi }{r{{\mathcal {V}}}_4} \frac{d}{d\ln \Lambda } \sum \limits _{\begin{array}{c} \mathrm{{vacuum}} \\ \mathrm{{supergraphs}} \end{array}} \sum \limits _{I=1}^{P'} \bigg [\,\theta ^4 (v^B)^2;\, \left[ \, T^A T^A\right] _I;\, \frac{1}{Q^{2n}}\bigg |_{n\ge 2} \nonumber \\&\quad \rightarrow 0;\, \frac{1}{Q_I^2} \rightarrow 4\pi ^2\delta ^4(Q_I) \bigg ]. \end{aligned}$$
(60)

Note that for the gauge and ghost propagators \((T^A T^A)_{BC} = C_2\delta _{BC}\), while for the matter propagators \((T^A T^A)_i{}^j = C(R)_i{}^j\).

Basing on Eq. (60) it is possible to modify the step 5 of the algorithm described in Sect. 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_I^{2n}\) with \(n\ge 2\) in scalar integrals and for any marked propagator with the momentum \(q_\mu \) in the remaining terms make the replacement

$$\begin{aligned}&\frac{1}{q^2}\quad \ \rightarrow \ \ C_2\, \frac{\partial ^2}{\partial q_\mu \partial q^\mu } \bigg (\frac{1}{q^2}\bigg )\qquad \qquad \text {for gauge and ghost propagators};\nonumber \\&\delta _i^j\, \frac{1}{q^2}\ \rightarrow \ \ C(R)_i{}^j\, \frac{\partial ^2}{\partial q_\mu \partial q^\mu } \bigg (\frac{1}{q^2}\bigg )\qquad \text {for matter propagators}. \end{aligned}$$
(61)

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 Sect. 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_I^{2n}\) with \(n\ge 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 \(\delta \)-singularities of the form (40). Various perturbative verifications of the above described modification will be made in Sect. 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 [64] 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. [53] for the Abelian case.

First, we introduce the auxiliary action

$$\begin{aligned} S_g \equiv S_{\mathrm{total}} +\Delta S_g \end{aligned}$$
(62)

depending on a real constant g, where \(S_{\mathrm{total}} \equiv S_{\mathrm{reg}} + S_{\mathrm{gf}} + S_{\mathrm{FP}} + S_{\mathrm{NK}}\) and

$$\begin{aligned} \Delta S_g\equiv & {} \frac{1}{4} \bigg (\frac{1}{g}-1\bigg ) \int d^8x \nonumber \\&\times \bigg [ - V^A R(\partial ^2/\Lambda ^2)\, \partial ^2 \Pi _{1/2} V^A -\frac{1}{8\xi _0} \nonumber \\&\times D^2 V^A K(\partial ^2/\Lambda ^2) {\bar{D}}^2 V^A + {\bar{c}}^{+A} c^A \nonumber \\&-{\bar{c}}^A c^{+A} + \phi ^{*i} F(\partial ^2/\Lambda ^2)\phi _i\bigg ]. \end{aligned}$$
(63)

For \(g=1\) the action \(S_g\) evidently coincides with \(S_{\mathrm{total}}\). Moreover, all vertices generated by \(S_g\) coincide with the ones produced by \(S_{\mathrm{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_{\mathrm{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

$$\begin{aligned} Z_g\equiv & {} \int D\mu \, \text {Det}(PV, M_\varphi )^{-1}\, \text {Det}(PV, M)^c \nonumber \\&\times \exp \bigg (i S_g + i S_{\mathrm{sources}}\bigg ), \end{aligned}$$
(64)

with \(c=T(R)/T(R_{\mathrm{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 \(\Gamma _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

$$\begin{aligned}&\frac{\partial \Gamma _g}{\partial g}\bigg |_{g=1;\, \mathrm{fields=0}} =\bigg \langle \frac{\partial S_g}{\partial g}\bigg \rangle \bigg |_{g=1,\, \mathrm{fields=0}} \nonumber \\&\quad = \bigg \langle \frac{1}{4} \int d^8x\, \bigg \{ V^A\bigg [ \partial ^2\Pi _{1/2} R(\partial ^2/\Lambda ^2) + \frac{1}{16\xi _0}\bigg ({\bar{D}}^2 D^2 \nonumber \\&\qquad + D^2 {\bar{D}}^2\bigg ) K(\partial ^2/\Lambda ^2) \bigg ]V^A - \phi ^{*i} F(\partial ^2/\Lambda ^2) \phi _i \nonumber \\&\qquad - {\bar{c}}^{+A} c^A + {\bar{c}}^A c^{+A} \bigg \} \bigg \rangle \bigg |_{g=1;\, \mathrm{fields=0}}. \end{aligned}$$
(65)

The angular brackets entering this equation are defined in the standard way

$$\begin{aligned} \langle A \rangle\equiv & {} \frac{1}{Z}\int D\mu \, A\, \text {Det}(PV, M_\varphi )^{-1}\, \text {Det}(PV, M)^c \nonumber \\&\times \exp \bigg (i S_{\mathrm{total}} + i S_{\mathrm{sources}}\bigg ), \end{aligned}$$
(66)

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,

$$\begin{aligned}&\frac{\partial \Gamma _g}{\partial g}\bigg |_{g=1;\, \mathrm{fields=0}} \nonumber \\&\quad = \frac{i}{4}\int d^8x\,\bigg \{\bigg [ \partial ^2\Pi _{1/2} R(\partial ^2/\Lambda ^2) + \frac{1}{16\xi _0} \nonumber \\&\qquad \times \bigg ({\bar{D}}^2 D^2 + D^2 {\bar{D}}^2\bigg ) K(\partial ^2/\Lambda ^2) \bigg ]_x \nonumber \\&\qquad \times \bigg (\frac{\delta ^2\Gamma }{\delta V^A_{x}\delta V^A_y}\bigg )^{-1} - F(\partial ^2/\Lambda ^2)_x \bigg (\frac{\delta ^2\Gamma }{\delta \phi ^{*i}_{\,,x} \delta \phi _{i,y}}\bigg )^{-1} \nonumber \\&\qquad - \bigg (\frac{\delta ^2\Gamma }{\delta c^A_{x}\,\delta {\bar{c}}^{+A}_{y}}\bigg )^{-1} - \bigg (\frac{\delta ^2\Gamma }{\delta {\bar{c}}^A_{x}\,\delta c^{+A}_{y}}\bigg )^{-1} \bigg \}\bigg |_{y=x}. \end{aligned}$$
(67)

(Certainly, in the last expression we also assume that \(g=1\) and the (super)fields are set to 0.)

The original effective action \(\Gamma \) 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 \(\partial \Gamma _g/\partial 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

$$\begin{aligned} G \equiv 1+ \Delta G \end{aligned}$$
(68)

depending on the momentum, then the inverse Green function will include

$$\begin{aligned} G^{-1} = 1 + \sum \limits _{n=1}^\infty (-1)^n \left( \Delta G\right) ^n. \end{aligned}$$
(69)

Evidently, a term containing \((\Delta G)^n\) corresponds to supergraphs of the structure presented in Fig. 2 with n coinciding momenta \(k_\mu \). The gray circles in this supergraph certainly give various \(\Delta G\).

As we discussed in the previous section, for constructing the \(\beta \)-function one should cut only one propagator with the momentum \(k_\mu \) in a supergraph which have the structure presented in Fig. 2. (This cut is made by \(\delta ^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_\mu \) is proportional to n. Therefore, for calculating the contribution to the \(\beta \)-function it is necessary to replace \(G^{-1}\) by

$$\begin{aligned} \sum \limits _{n=1}^\infty \frac{(-1)^{n}}{n} \left( \Delta G\right) ^n = -\ln G. \end{aligned}$$
(70)

(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 Sect. 4. Details of this calculation will be considered below, separately for the cuts of gauge, ghost, and matter propagators.

5.2 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 \(\partial \Gamma _g/\partial g\) is given by the expression

$$\begin{aligned} -\frac{i}{4}\int d^8x\, F(\partial ^2/\Lambda ^2)_x \bigg (\frac{\delta ^2\Gamma }{\delta \phi ^{*i}_{\,,x}\delta \phi _{i,y}}\bigg )^{-1}\bigg |_{y=x;\ \mathrm{fields}=0},\nonumber \\ \end{aligned}$$
(71)

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 \(\phi \) and \(\phi ^*\) reads as

$$\begin{aligned} \frac{\delta ^2\Gamma }{\delta \phi ^{*i}_{\,,x}\, \delta \phi _{j,y}}\bigg |_{\mathrm{fields}=0} = \frac{1}{16} D_x^2 {\bar{D}}_y^2 \left( G_\phi \right) _i{}^j \delta ^8_{xy}, \end{aligned}$$
(72)

where

$$\begin{aligned} \delta ^8_{xy} \equiv \delta ^4(x-y)\, \delta ^4(\theta _x-\theta _y). \end{aligned}$$
(73)

The function \((G_\phi )_i{}^j\) present in Eq. (72) is normalized in such a way that in the tree approximation it is equal to \(\delta _i^j F(\partial ^2/\Lambda ^2)\). (In the limit \(\Lambda \rightarrow \infty \) this expression gives \(\delta _i^j\).) Therefore, it can be presented as

$$\begin{aligned} \left( G_\phi \right) _i{}^j = \delta _i^j F + \left( \Delta G_\phi \right) _i{}^j, \end{aligned}$$
(74)

where the sum of quantum corrections is denoted by \((\Delta G_\phi )_i{}^j\). The function inverse to (72) entering Eq. (71) is given by the expression

$$\begin{aligned} \bigg (\frac{\delta ^2\Gamma }{\delta \phi ^{*i}_{\,,x}\, \delta \phi _{j,y}}\bigg )^{-1}= & {} -\frac{D_x^2 {\bar{D}}_y^2}{4\partial ^2} \left( G_\phi ^{-1}\right) _j{}^i \delta ^8_{xy} \nonumber \\= & {} -4 \bigg (\frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F} \delta _j^i + \frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F}\, \left( \Delta G_\phi \right) _j{}^i\, \nonumber \\&\times \frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F} + \frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F}\, \left( \Delta G_\phi \right) _j{}^k\, \frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F}\nonumber \\&\times \left( \Delta G_\phi \right) _k{}^i\, \frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F} +\cdots \bigg ) \delta ^8_{xy}. \end{aligned}$$
(75)

This can easily be verified using the first of the identities

$$\begin{aligned} {\bar{D}}^2 D^2 {\bar{D}}^2 = -16 {\bar{D}}^2\partial ^2;\qquad D^2 {\bar{D}}^2 D^2 = -16 D^2\partial ^2. \end{aligned}$$
(76)

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 \(\Delta G_\phi \) 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

$$\begin{aligned}&-4 \bigg (\frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F} \left( \Delta G_\phi \right) _j{}^i \frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F} + \frac{1}{2}\cdot \frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F} \nonumber \\&\quad \qquad \times \left( \Delta G_\phi \right) _j{}^k \frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F} \left( \Delta G_\phi \right) _k{}^i \frac{D_x^2 {\bar{D}}_x^2}{16\partial ^2 F} + \cdots \bigg ) \delta ^8_{xy} \nonumber \\&\quad = \frac{D_x^2 {\bar{D}}_x^2}{4\partial ^2 F} \nonumber \\&\qquad \times \bigg (\frac{\Delta G_\phi }{F} - \frac{(\Delta G_\phi )^2}{2 F^2} + \frac{(\Delta G_\phi )^3}{3 F^3} + \cdots \bigg )_j{}^i \delta ^8_{xy} \nonumber \\&\quad = \frac{D_x^2 {\bar{D}}_x^2}{4\partial ^2 F} \bigg (\ln \frac{G_\phi }{F}\bigg )_j{}^i \delta ^8_{xy}. \end{aligned}$$
(77)

Therefore, the first step for calculating the matter contribution to the function (34) is to make the replacement

$$\begin{aligned} \bigg (\frac{\delta ^2\Gamma }{\delta \phi ^{*i}_{\,,x}\, \delta \phi _{i,y}}\bigg )^{-1}\bigg |_{\mathrm{fields}=0}\ \rightarrow \ \frac{D_x^2 {\bar{D}}_x^2}{4\partial ^2 F} \bigg (\ln \frac{G_\phi }{F}\bigg )_i{}^i \delta ^8_{xy} \end{aligned}$$
(78)

in Eq. (71). We see that no bad terms proportional to \(1/q^{2k}\) with \(k\ge 2\) appear. Next, following the algorithm described in Sect. 4, we should insert \(\theta ^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

$$\begin{aligned}&\bigg (\ln \frac{G_\phi }{F}\bigg )_i{}^i = \bigg (\ln \frac{G_\phi }{F}\bigg )_j{}^i \delta _i^j \nonumber \\&\quad \rightarrow \ \bigg (\ln \frac{G_\phi }{F}\bigg )_j{}^i C(R)_i{}^j \frac{\partial ^2}{\partial q_\mu \partial q^\mu }, \end{aligned}$$
(79)

where \(q_\mu \) is the Minkowski momentum of the matter line which is cut. Note that (as we discussed in Sect. 4.2) the operator \(\partial ^2/\partial q_\mu \partial q^\mu \) should act only on \(1/q^2\). (After the Wick rotation it gives \(-\partial ^2/\partial Q_\mu ^2\), where \(Q_\mu \) is the corresponding Euclidean momentum.) The matter line is cut in the point x, so that \(q_\mu \) is the momentum of the propagator coming out of this point. Finally the result should be multiplied by \(-2\pi /(r{{\mathcal {V}}}_4)\) and differentiated with respect to \(\ln \Lambda \) at fixed values of renormalized couplings. This implies that the matter contribution to the function (34) can be written in the form

$$\begin{aligned} \Delta _{\mathrm{matter}}\bigg (\frac{\beta }{\alpha _0^2}\bigg )= & {} - \frac{i\pi }{8 r{{\mathcal {V}}}_4} C(R)_i{}^j \frac{d}{d\ln \Lambda } \int d^8x\,d^8y\, \left( \theta ^4\right) _x \left( v^B\right) ^2_x\nonumber \\&\times \int \frac{d^4q}{(2\pi )^4}\, \delta ^8_{xy}(q)\, \bigg (\ln \frac{G_\phi }{F}\bigg )_j{}^i \nonumber \\&\times \frac{\partial ^2}{\partial q_\mu \partial q^\mu } \bigg (\frac{1}{q^2}\bigg ) D_x^2 {\bar{D}}_x^2 \delta ^8_{xy}, \end{aligned}$$
(80)

where the function \((\ln G_\phi )_i{}^j\) depends on the momentum \(q_\mu \) and

$$\begin{aligned} \delta ^8_{xy}(q) \equiv \delta ^4(\theta _x-\theta _y) e^{iq_\mu (x^\mu - y^\mu )}. \end{aligned}$$
(81)

The momentum integral is calculated in the Euclidean space after the Wick rotation,

$$\begin{aligned}&\frac{d}{d\ln \Lambda } \int \frac{d^4q}{(2\pi )^4}\, \delta ^8_{xy}(q)\, \bigg (\ln \frac{G_\phi }{F}\bigg )_j{}^i \frac{\partial ^2}{\partial q_\mu \partial q^\mu } \bigg (\frac{1}{q^2}\bigg ) \nonumber \\&\quad \rightarrow \ - 4i\pi ^2 \delta ^4(\theta _x-\theta _y) \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4}\, \bigg (\ln \frac{G_\phi }{F}\bigg )_j{}^i\, \delta ^4(Q)\nonumber \\&\quad = - \frac{i}{4\pi ^2} \delta ^4(\theta _x-\theta _y) \frac{d}{d\ln \Lambda }\left( \ln G_\phi \right) _j{}^i\bigg |_{Q=0} \nonumber \\&\quad = - \frac{i}{4\pi ^2} \delta ^4(\theta _x-\theta _y) \left( \gamma _\phi \right) _j{}^i(\alpha _0,\lambda _0,Y_0),\qquad \end{aligned}$$
(82)

where we took into account that \(F(0)=1\) and wrote the result in terms of the anomalous dimension \((\gamma _\phi )_j{}^i\) using Eq. (32). Therefore, the expression (80) takes the form

$$\begin{aligned} \Delta _{\mathrm{matter}}\bigg (\frac{\beta }{\alpha _0^2}\bigg )&= - \frac{1}{32\pi r{{\mathcal {V}}}_4} C(R)_i{}^j \left( \gamma _\phi \right) _j{}^i(\alpha _0,\lambda _0,Y_0) \nonumber \\&\quad \times \int d^8x\,d^8y\,\left( \theta ^4\right) _x \left( v^B\right) ^2_x\nonumber \\&\quad \times \delta ^4(\theta _x-\theta _y) D_x^2 {\bar{D}}_x^2 \delta ^8_{xy}. \end{aligned}$$
(83)

Next, it is necessary to use the identity

$$\begin{aligned} \delta ^4(\theta _x-\theta _y) D_x^2 {\bar{D}}_x^2 \delta ^8_{xy} = 4\, \delta ^8_{xy} \end{aligned}$$
(84)

and calculate the integral over \(d^8y\) with the help of the \(\delta \)-function,

$$\begin{aligned} \Delta _{\mathrm{matter}}\bigg (\frac{\beta }{\alpha _0^2}\bigg )= & {} - \frac{1}{8\pi r{{\mathcal {V}}}_4} C(R)_i{}^j \left( \gamma _\phi \right) _j{}^i(\alpha _0,\lambda _0,Y_0) \nonumber \\&\times \int d^8x\,\left( \theta ^4\right) _x \left( v^B\right) ^2_x. \end{aligned}$$
(85)

Taking into account that

$$\begin{aligned} \int d^8x\, \left( \theta ^4\right) _x \left( v^B\right) ^2_x = 4{{\mathcal {V}}}_4, \end{aligned}$$
(86)

we obtain the final expression for the sum of matter singularities

$$\begin{aligned} \Delta _{\mathrm{matter}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) = - \frac{1}{2\pi r} C(R)_i{}^j \left( \gamma _\phi \right) _j{}^i(\alpha _0,\lambda _0,Y_0). \end{aligned}$$
(87)

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 [64] by a different method. This can be considered as a correctness test of the method proposed in this paper.

5.3 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,

$$\begin{aligned} -\frac{i}{4}\int d^8x\,\bigg \{ \bigg (\frac{\delta ^2\Gamma }{\delta c^A_{x}\,\delta {\bar{c}}^{+A}_{y}}\bigg )^{-1} + \bigg (\frac{\delta ^2\Gamma }{\delta {\bar{c}}^A_{x}\,\delta c^{+A}_{y}}\bigg )^{-1} \bigg \}\bigg |_{y=x;\ \mathrm{fields}=0}.\nonumber \\ \end{aligned}$$
(88)

The two-point Green functions of the Faddeev–Popov ghosts are obtained by differentiating Eq. (24),

$$\begin{aligned} \frac{\delta ^2\Gamma }{\delta {\bar{c}}^{A}_{x}\, \delta c^{+B}_{y}}\bigg |_{\mathrm{fields}=0}= & {} \frac{\delta ^2\Gamma }{\delta c^{A}_{x}\, \delta {\bar{c}}^{+B}_{y}}\bigg |_{\mathrm{fields}=0} \nonumber \\= & {} \delta _{AB} G_c \frac{{\bar{D}}_x^2 D^2_y}{16} \delta ^8_{xy}. \end{aligned}$$
(89)

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 + \Delta G_c\), where \(\Delta 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

$$\begin{aligned}&\int d^6z\, \bigg (\frac{\delta ^2\Gamma }{\delta c^{+A}_{x}\, \delta {\bar{c}}^{D}_{z}}\bigg )^{-1} \frac{\delta ^2\Gamma }{\delta {\bar{c}}^{D}_{z}\, \delta c^{+B}_{y}}\bigg |_{\mathrm{fields}=0} \nonumber \\&\quad = -\frac{1}{2}\delta _{AB} D^2_x \delta ^8_{xy};\nonumber \\&\int d^6z\, \bigg (\frac{\delta ^2\Gamma }{\delta {\bar{c}}^{+A}_{x}\, \delta c^{D}_{z}}\bigg )^{-1} \frac{\delta ^2\Gamma }{\delta c^{D}_{z}\, \delta {\bar{c}}^{+B}_{y}}\bigg |_{\mathrm{fields}=0} \nonumber \\&\quad = -\frac{1}{2}\delta _{AB} D^2_x \delta ^8_{xy} \end{aligned}$$
(90)

and are explicitly written as

$$\begin{aligned}&\bigg (\frac{\delta ^2\Gamma }{\delta {\bar{c}}^{A}_{x}\, \delta c^{+B}_{y}}\bigg )^{-1}\bigg |_{\mathrm{fields}=0} = \bigg (\frac{\delta ^2\Gamma }{\delta c^{A}_{x}\, \delta {\bar{c}}^{+B}_{y}}\bigg )^{-1}\bigg |_{\mathrm{fields}=0} \nonumber \\&\quad = \delta _{AB}\,\frac{{\bar{D}}^2_x D^2_y}{4\partial ^2 G_c} \delta ^8_{xy}. \end{aligned}$$
(91)

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

$$\begin{aligned}&\bigg (\frac{\delta ^2\Gamma }{\delta c^A_{x}\, \delta {\bar{c}}^{+A}_{y}}\bigg )^{-1}\bigg |_{\mathrm{fields}=0}\ \rightarrow \ - \delta _{AA} \frac{{\bar{D}}_x^2 D_y^2}{4\partial ^2} \ln G_c\, \delta ^8_{xy};\qquad \nonumber \\&\bigg (\frac{\delta ^2\Gamma }{\delta {\bar{c}}^A_{x}\,\delta c^{+A}_{y}}\bigg )^{-1}\bigg |_{\mathrm{fields}=0}\ \rightarrow \ - \delta _{AA} \frac{{\bar{D}}_x^2 D_y^2}{4\partial ^2} \ln G_c\, \delta ^8_{xy} \end{aligned}$$
(92)

in the expression (88). Next, according to the algorithm described in Sect. 4, \(\theta ^4 (v^B)^2\) should be inserted to the point x. After this, it is necessary to make the replacementFootnote 9

$$\begin{aligned} \ln G_c\cdot \delta _{AA}\ \rightarrow \ \ln G_c\cdot r C_2\, \frac{\partial ^2}{\partial q_\mu \partial q^\mu } \end{aligned}$$
(93)

and apply the operator (38) to the resulting expression. Note that the differentiation with respect to \(\ln \Lambda \) 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)

$$\begin{aligned} \Delta _{\mathrm{ghost}}\bigg (\frac{\beta }{\alpha _0^2}\bigg )= & {} \frac{i\pi C_2}{4{{\mathcal {V}}}_4} \frac{d}{d\ln \Lambda } \int d^8x\,d^8y\, \left( \theta ^4\right) _x \left( v^B\right) ^2_x\nonumber \\&\times \int \frac{d^4q}{(2\pi )^4}\, \delta ^8_{xy}(q) \ln G_c\, \frac{\partial ^2}{\partial q_\mu \partial q^\mu } \bigg (\frac{1}{q^2}\bigg ) \nonumber \\&\times {\bar{D}}_x^2 D_x^2 \delta ^8_{xy}. \end{aligned}$$
(94)

(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,

$$\begin{aligned}&\frac{d}{d\ln \Lambda } \int \frac{d^4q}{(2\pi )^4}\, \delta ^8_{xy}(q) \ln G_c\, \frac{\partial ^2}{\partial q_\mu \partial q^\mu } \bigg (\frac{1}{q^2}\bigg ) \nonumber \\&\quad \rightarrow \ - \frac{i}{4\pi ^2} \delta ^4(\theta _x-\theta _y) \gamma _c(\alpha _0,\lambda _0,Y_0). \end{aligned}$$
(95)

With the help of this equation the considered contribution to the function (34) can be brought to the form

$$\begin{aligned} \Delta _{\mathrm{ghost}}\bigg (\frac{\beta }{\alpha _0^2}\bigg )= & {} \frac{C_2}{16\pi {{\mathcal {V}}}_4} \gamma _c(\alpha _0,\lambda _0,Y_0) \nonumber \\&\times \int d^8x\,d^8y\,\left( \theta ^4\right) _x \left( v^B\right) ^2_x\, \delta ^4(\theta _x-\theta _y) \nonumber \\&\times {\bar{D}}_x^2 D_x^2 \delta ^8_{xy}. \end{aligned}$$
(96)

Then we use the identity (84) and calculate the integral over \(d^8y\) with the help of resulting \(\delta ^8_{xy}\),

$$\begin{aligned} \Delta _{\mathrm{ghost}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) = \frac{C_2}{4\pi {{\mathcal {V}}}_4} \gamma _c(\alpha _0,\lambda _0,Y_0) \int d^8x\,\bigg (\theta ^4\bigg )_x \bigg (v^B\bigg )^2_x.\nonumber \\ \end{aligned}$$
(97)

According to Eq. (86), the remaining integral present in this equation is equal to \(4{{\mathcal {V}}}_4\), so that the sum of singular contributions produced by cuts of the Faddeev–Popov ghost propagators takes the form

$$\begin{aligned} \Delta _{\mathrm{ghost}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) = \frac{C_2}{\pi } \gamma _c(\alpha _0,\lambda _0,Y_0). \end{aligned}$$
(98)

This result agrees with the one found in Ref. [64] by a different method and coincides with the term containing the anomalous dimension \(\gamma _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 antighost \({\bar{c}}\)) in the adjoint representation of the gauge group.

5.4 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

$$\begin{aligned} \Gamma ^{(2)}_V= & {} -\frac{1}{4} \int d^8x\, V^A G_V \partial ^2\Pi _{1/2} V^A \nonumber \\&- \frac{1}{32 \xi _0} \int d^8x\, D^2 V^A K {\bar{D}}^2 V^A, \end{aligned}$$
(99)

where the regulator function K depends on \(\partial ^2/\Lambda ^2\), and the argument \(q^2\) of function \(G_V\) should be replaced by \(-\partial ^2\). Differentiating this expression with respect to the quantum gauge superfield, we obtain the corresponding two-point Green function

$$\begin{aligned}&\frac{\delta ^2\Gamma }{\delta V_x^A \delta V_y^B}\bigg |_{\mathrm{fields}=0} \nonumber \\&\quad = \delta _{AB} \bigg [ - \frac{1}{2} G_V \partial ^2\Pi _{1/2} - \frac{1}{32 \xi _0} K \bigg (D^2{\bar{D}}^2 + {\bar{D}}^2 D^2\bigg )\bigg ] \delta ^8_{xy}.\nonumber \\ \end{aligned}$$
(100)

By definition, the inverse function satisfies the equation

$$\begin{aligned} \int d^8z\,\bigg (\frac{\delta ^2\Gamma }{\delta V_x^A \delta V_z^C}\bigg )^{-1} \frac{\delta ^2\Gamma }{\delta V_z^C \delta V_y^B}\bigg |_{\mathrm{fields}=0} = \delta _{AB} \delta ^8_{xy}. \end{aligned}$$
(101)

Using the identity (see, e.g., [80])

$$\begin{aligned} 0 = \partial ^2 + \partial ^2\Pi _{1/2} + \frac{1}{16}\bigg (D^2 {\bar{D}}^2 + {\bar{D}}^2 D^2\bigg ), \end{aligned}$$
(102)

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

$$\begin{aligned}&\bigg (\frac{\delta ^2\Gamma }{\delta V_x^A \delta V_y^B}\bigg )^{-1}\bigg |_{\mathrm{fields}=0} \nonumber \\&\quad = - \delta _{AB}\, \bigg [\,\frac{2}{G_V \partial ^4} \partial ^2\Pi _{1/2} + \frac{\xi _0}{8 \partial ^4 K} \bigg (D^2{\bar{D}}^2 + {\bar{D}}^2 D^2\bigg )\bigg ] \delta ^8_{xy}.\nonumber \\ \end{aligned}$$
(103)

Thus, we see that the effective propagator is proportional to \(Q^{-4}\), where \(Q_\mu \) is the corresponding Euclidean momentum.

5.5 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

$$\begin{aligned}&\frac{i}{4}\int d^8x \nonumber \\&\qquad \times \bigg [ \partial ^2\Pi _{1/2} R(\partial ^2/\Lambda ^2) + \frac{1}{16\xi _0}\bigg ({\bar{D}}^2 D^2 + D^2 {\bar{D}}^2\bigg ) K(\partial ^2/\Lambda ^2) \bigg ]_x \nonumber \\&\qquad \times \bigg (\frac{\delta ^2\Gamma }{\delta V^A_{x}\delta V^A_y}\bigg )^{-1} \bigg |_{y=x;\ \mathrm{fields}=0}, \end{aligned}$$
(104)

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(\partial ^2/\Lambda ^2)\) in the tree approximation. That is why it is convenient to present it in the form

$$\begin{aligned} G_V = R + \Delta G_V, \end{aligned}$$
(105)

where \(\Delta 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

$$\begin{aligned} \left( \Pi _{1/2}\right) ^2 = - \Pi _{1/2};\qquad \Pi _{1/2} D^2 = 0;\qquad \Pi _{1/2} {\bar{D}}^2 = 0,\nonumber \\ \end{aligned}$$
(106)

and (76). To simplify the resulting expression, we introduce the notation

$$\begin{aligned} P \equiv \frac{1}{\partial ^4 R} \partial ^2\Pi _{1/2} + \frac{\xi _0}{16 \partial ^4 K} \bigg (D^2{\bar{D}}^2 + {\bar{D}}^2 D^2\bigg ). \end{aligned}$$
(107)

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

$$\begin{aligned}&\bigg (\frac{\delta ^2\Gamma }{\delta V_x^A \delta V_y^A}\bigg )^{-1}\bigg |_{\mathrm{fields}=0} \nonumber \\&\quad = - r\, \bigg [\,\frac{2}{G_V \partial ^4} \partial ^2\Pi _{1/2} + \frac{\xi _0}{8 \partial ^4 K} \bigg (D^2{\bar{D}}^2 + {\bar{D}}^2 D^2\bigg )\bigg ] \delta ^8_{xy} \nonumber \\&\quad = - 2 r\,P\bigg (1 - \partial ^2\Pi _{1/2} \Delta G_V P + \partial ^2\Pi _{1/2} \Delta G_V P \nonumber \\&\qquad \times \partial ^2\Pi _{1/2} \Delta G_V P - \cdots \bigg )\delta ^8_{xy}. \end{aligned}$$
(108)

The factors \(\partial ^2\Pi _{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 Sect. 5.1, for constructing a contribution to the function (34), one should first divide a term with the n-th power of \(\Delta 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

$$\begin{aligned}&- 2 r\, P\bigg ( - \partial ^2\Pi _{1/2} \Delta G_V P + \frac{1}{2}\partial ^2\Pi _{1/2} \Delta G_V P \nonumber \\&\qquad \times \partial ^2\Pi _{1/2} \Delta G_V P - \cdots \bigg )\delta ^8_{xy}\nonumber \\&\quad = - 2 r\, P \sum \limits _{n=1}^\infty \frac{1}{n} \bigg (-\partial ^2\Pi _{1/2} \Delta G_V P\bigg )^n \delta ^8_{xy} \nonumber \\&\quad = - 2 r\, \frac{1}{\partial ^4 R} \partial ^2\Pi _{1/2}\sum \limits _{n=1}^\infty \frac{1}{n} \bigg ( -\frac{\Delta G_V}{R}\bigg )^n \delta ^8_{xy}. \end{aligned}$$
(109)

Calculating the sum of this series, we obtain that, at the first step, it is necessary to make in Eq. (104) the formal substitution

$$\begin{aligned} \bigg (\frac{\delta ^2\Gamma }{\delta V_x^A \delta V_y^A}\bigg )^{-1}\bigg |_{\mathrm{fields}=0}\ \rightarrow \ \frac{2 r}{\partial ^4 R}\, \partial ^2\Pi _{1/2} \ln \frac{G_V}{R} \delta ^8_{xy}. \end{aligned}$$
(110)

After this, with the help of Eq. (106) we obtain

$$\begin{aligned}&(104)\ \rightarrow \ -\frac{i r}{2}\int d^8x\, \partial ^2\Pi _{1/2} \frac{1}{\partial ^2}\, \ln \frac{G_V}{R} \delta ^8_{xy}\bigg |_{y=x} \nonumber \\&\quad = \frac{i r}{2}\int d^8x\, d^8y\, \int \frac{d^4q}{(2\pi )^4} \delta ^8_{xy}(q)\, \partial ^2\Pi _{1/2} \frac{1}{q^2}\, \ln \frac{G_V}{R} \delta ^8_{xy}.\nonumber \\ \end{aligned}$$
(111)

It should be noted that all possible bad terms proportional to \(1/q^{2k}\) with \(k\ge 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 \(\theta ^4 (v^B)^2\) to the point x, apply the operator \(C_2 \partial ^2/\partial q_\mu \partial q^\mu \) to \(1/q^2\) coming from the gauge propagator which is cut in the point x, multiply the result by \(-2\pi /(r{{\mathcal {V}}}_4)\), and differentiate it with respect to \(\ln \Lambda \). The sum of the gauge singularities constructed according to this procedure is given by

$$\begin{aligned}&\Delta _{\mathrm{gauge}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) \nonumber \\&\quad = -\frac{i\pi C_2}{{{\mathcal {V}}}_4} \frac{d}{d\ln \Lambda } \int d^8x\,d^8y\, \left( \theta ^4\right) _x \left( v^B\right) ^2_x\nonumber \\&\qquad \times \int \frac{d^4q}{(2\pi )^4}\, \delta ^8_{xy}(q)\, \ln \frac{G_V}{R}\, \frac{\partial ^2}{\partial q_\mu \partial q^\mu }\left( \frac{1}{q^2}\right) \nonumber \\&\qquad \times \left( \partial ^2\Pi _{1/2}\right) _x \delta ^8_{xy}. \end{aligned}$$
(112)

As earlier, we calculate the momentum integral in the Euclidean space after the Wick rotation taking into account that \(R(0)=1\),

$$\begin{aligned}&\frac{d}{d\ln \Lambda } \int \frac{d^4q}{(2\pi )^4}\, \delta ^8_{xy}(q)\, \ln \frac{G_V}{R}\, \frac{\partial ^2}{\partial q_\mu \partial q^\mu } \bigg (\frac{1}{q^2}\bigg )\nonumber \\&\quad \ \rightarrow \ -\frac{i}{4\pi ^2}\, \delta ^4(\theta _x-\theta _y)\, \frac{d}{d\ln \Lambda } \ln G_V\bigg |_{Q=0} \nonumber \\&\quad = -\frac{i}{2\pi ^2}\, \delta ^4(\theta _x-\theta _y)\, \gamma _V(\alpha _0,\lambda _0,Y_0). \end{aligned}$$
(113)

This implies that

$$\begin{aligned}&\Delta _{\mathrm{gauge}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) \nonumber \\&\quad = - \frac{C_2}{2\pi {{\mathcal {V}}}_4} \gamma _V(\alpha _0,\lambda _0,Y_0) \nonumber \\&\qquad \times \int d^8x\,d^8y\, \left( \theta ^4\right) _x \left( v^B\right) ^2_x\, \delta ^4(\theta _x-\theta _y) \left( \partial ^2\Pi _{1/2}\right) _x \delta ^8_{xy}.\nonumber \\ \end{aligned}$$
(114)

For calculating the remaining integral we use the identity

$$\begin{aligned} \delta ^4(\theta _x-\theta _y) \left( \partial ^2\Pi _{1/2}\right) _x \delta ^8_{xy}= & {} -\frac{1}{8}\, \delta ^4(\theta _x-\theta _y) \left( D^a {\bar{D}}^2 D_a\right) _x \delta ^8_{xy} \nonumber \\= & {} - \frac{1}{2}\,\delta ^8_{xy}. \end{aligned}$$
(115)

Therefore,

$$\begin{aligned}&\int d^8x\,d^8y\, \left( \theta ^4\right) _x \left( v^B\right) ^2_x\, \delta ^4(\theta _x-\theta _y) \left( \partial ^2\Pi _{1/2}\right) _x \delta ^8_{xy} \nonumber \\&\quad = - \frac{1}{2} \int d^8x\, d^8y\, \left( \theta ^4\right) _x \left( v^B\right) ^2_x\, \delta ^8_{xy} = - 2 {{\mathcal {V}}}_4. \end{aligned}$$
(116)

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,

$$\begin{aligned} \Delta _{\mathrm{gauge}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) = \frac{C_2}{\pi } \gamma _V(\alpha _0,\lambda _0,Y_0). \end{aligned}$$
(117)

We see that it coincides with the corresponding term in Eq. (28) in exact agreement with the guess made in Ref. [2].

5.6 The \(\beta \)-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,

$$\begin{aligned}&\frac{\beta (\alpha _0,\lambda _0,Y_0)}{\alpha _0^2} - \frac{\beta _{\mathrm{1-loop}}(\alpha _0)}{\alpha _0^2} \nonumber \\&\quad = \Delta _{\mathrm{matter}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) + \Delta _{\mathrm{ghost}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) + \Delta _{\mathrm{gauge}}\bigg (\frac{\beta }{\alpha _0^2}\bigg )\qquad \nonumber \\&\quad = - \frac{1}{2\pi r} C(R)_i{}^j \left( \gamma _\phi \right) _j{}^i(\alpha _0,\lambda _0,Y_0) \nonumber \\&\qquad + \frac{1}{\pi } C_2 \gamma _c(\alpha _0,\lambda _0,Y_0) + \frac{1}{\pi } C_2 \gamma _V(\alpha _0,\lambda _0,Y_0). \end{aligned}$$
(118)

Taking into account that the one-loop contribution to the \(\beta \)-function is given by Eq. (35), we see that the resulting expression for the \(\beta \)-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 [57] due to a different structure of loop integrals [67].

Note that, in fact, in the previous sections we have also proved the Eq. (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 \Lambda \) at fixed values of renormalized couplings and obtain the equation

$$\begin{aligned}&\beta (\alpha _0,\lambda _0,Y_0) \nonumber \\&\quad = 2\alpha _0\bigg (\gamma _c(\alpha _0,\lambda _0,Y_0) + \gamma _V(\alpha _0,\lambda _0,Y_0)\bigg ). \end{aligned}$$
(119)

Excluding the sum \(\gamma _c+\gamma _V\) from Eqs. (28) and (119) we obtain the NSVZ relation in the original form of the relation between the \(\beta \)-function and the anomalous dimension of the matter superfields,

$$\begin{aligned}&\frac{\beta (\alpha _0,\lambda _0,Y_0)}{\alpha _0^2} \nonumber \\&\quad = - \frac{3 C_2 - T(R) + C(R)_i{}^j \left( \gamma _\phi \right) _j{}^i(\alpha _0,\lambda _0,Y_0)/r}{2\pi (1- C_2\alpha _0/2\pi )}.\nonumber \\ \end{aligned}$$
(120)

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. Note that this regularization is not uniquely defined, because for introducing it we use three arbitrary higher-derivative regulator functions R(x), F(x), and K(x), see Eqs. (12) and (17), and two numerical parameters \(a=M/\Lambda \) and \(a_\varphi = M_\varphi /\Lambda \), which are the ratios of the Pauli–Villars masses to the dimensionful parameter in the higher derivative terms. The NSVZ equations (28) and (120) for RGFs defined in terms of the bare couplings hold for any choice of all these functions (provided they are equal to 1 at \(x=0\) and rapidly tend to infinity at \(x\rightarrow \infty \)) and parameters. However, all RGFs entering these equations in general can depend on the functions R(x) and F(x). Also they can depend on a and \(a_\varphi \). The explicit form of this dependence for the two-loop anomalous dimension and for the three-loop \(\beta \)-functionFootnote 10 can be found in Ref. [94]. The function K(x) appears in the gauge fixing term and is also present in the action for the Nielsen–Kallosh ghosts. The calculations in the lowest orders (see, e.g., [74, 93]) demonstrate that the \(\beta \)-function and the anomalous dimension of the matter superfields do not depend on it. Possibly, this is a consequence of general theorems about the gauge dependence of the effective action (see, e.g., [95] and references therein), although the detailed superfield analyses in the supersymmetric case has not yet been done. Note that the gauge independence of the matter superfield anomalous dimension presumably follows from the fact that in supersymmetric theories it is related to the mass anomalous dimension due to the non-renormalization of the superpotential [1]. However, the anomalous dimensions \(\gamma _c\) and \(\gamma _V\) seem to depend on the function K.Footnote 11 Nevertheless, both sides of the NSVZ equations have the same dependence on all regulator functions and parameters, so that the results discussed above are valid independently of their particular form.

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.

6 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 Sect. 2. Note that these RGFs are scheme-independent [56], 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 \(\overline{\text {DR}}\)-scheme is not the NSVZ scheme [27,28,29]. 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

$$\begin{aligned} {\widetilde{\beta }}(\alpha ,\lambda ,Y)\equiv & {} \left. \frac{d \alpha }{d\ln \mu }\right| _{\alpha _0,\lambda _0,Y_0 = \mathrm{const}}; \nonumber \\ {\widetilde{\gamma }}_V(\alpha ,\lambda ,Y)\equiv & {} \left. \frac{d\ln Z_V}{d\ln \mu }\right| _{\alpha _0,\lambda _0,Y_0 = \mathrm{const}};\nonumber \\ {\widetilde{\gamma }}_c(\alpha ,\lambda ,Y)\equiv & {} \left. \frac{d\ln Z_c}{d\ln \mu }\right| _{\alpha _0,\lambda _0,Y_0 = \mathrm{const}}; \nonumber \\ ({\widetilde{\gamma }}_\phi )_i{}^j(\alpha ,\lambda ,Y)\equiv & {} \left. \frac{d(\ln Z_\phi )_i{}^j}{d\ln \mu }\right| _{\alpha _0,\lambda _0,Y_0 = \mathrm{const}}. \end{aligned}$$
(121)

Unlike Eqs. (29)–(32), the derivatives in Eq. (121) are taken with respect to \(\ln \mu \) (instead of \(\ln \Lambda \)) at fixed values of the bare couplings (instead of the renormalized ones).

The key observation made in [56] 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 \Lambda /\mu \)) take the form

$$\begin{aligned} \alpha (\alpha _0,\lambda _0,Y_0,\ln \Lambda /\mu \rightarrow x_0)= & {} \alpha _0; \nonumber \\ Z_V(\alpha _0,\lambda _0,Y_0,\ln \Lambda /\mu \rightarrow x_0)= & {} 1; \nonumber \\ Z_c(\alpha _0,\lambda _0,Y_0,\ln \Lambda /\mu \rightarrow x_0)= & {} 1; \nonumber \\ (Z_\phi )_i{}^j(\alpha _0,\lambda _0,Y_0,\ln \Lambda /\mu \rightarrow x_0)= & {} \delta _i^j; \nonumber \\ Y(\alpha _0,\lambda _0,Y_0,\ln \Lambda /\mu \rightarrow x_0)= & {} Y_0. \end{aligned}$$
(122)

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

$$\begin{aligned} \frac{d}{d\ln \Lambda } (Z_\alpha ^{-1/2} Z_c Z_V) = 0 \end{aligned}$$
(123)

has the only solution \(Z_\alpha ^{-1/2} Z_c Z_V=1\). Similarly, it is easy to see that the condition \(Z_\xi Z_V^2=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 \Lambda /\mu \) and all finite constants which define a subtraction scheme vanish. The schemes corresponding to \(x_0\ne 0\) are related to the HD+MSL scheme by the redefinition of the renormalization point \(\mu = \exp (-x_0)\mu _{\mathrm{HD+MSL}}\) exactly as in the Abelian case considered in [36].

Under the boundary conditions (122) RGFs (121) can be related to RGFs (29)–(32) by the equations

$$\begin{aligned}&{\widetilde{\beta }}(\alpha ,\lambda ,Y) = \beta (\alpha _0\rightarrow \alpha , \lambda _0\rightarrow \lambda ,Y_0\rightarrow Y); \nonumber \\&{\widetilde{\gamma }}_V(\alpha ,\lambda ,Y) = \gamma _V(\alpha _0\rightarrow \alpha , \lambda _0\rightarrow \lambda ,Y_0\rightarrow Y); \nonumber \\&{\widetilde{\gamma }}_c(\alpha ,\lambda ,Y) = \gamma _c(\alpha _0\rightarrow \alpha , \lambda _0\rightarrow \lambda ,Y_0\rightarrow Y); \nonumber \\&({\widetilde{\gamma }}_\phi )_i{}^j(\alpha ,\lambda ,Y) = (\gamma _\phi )_i{}^j(\alpha _0\rightarrow \alpha , \lambda _0\rightarrow \lambda ,Y_0\rightarrow Y),\nonumber \\ \end{aligned}$$
(124)

where the arrows point that it is necessary to make a formal change of the argument, say, to write \(\alpha \) instead of \(\alpha _0\), etc.

All these equations can be proved repeating the argumentation of [56] (in which similar equations were derived in the Abelian case). For example, taking into account that \(\alpha = \alpha _0 Z_\alpha \), we obtain

$$\begin{aligned} {\widetilde{\beta }}(\alpha ,\lambda ,Y)= & {} \frac{d\alpha }{d\ln \mu } \bigg |_{\alpha _0,\lambda _0,Y_0 = \mathrm{const}} \nonumber \\= & {} \frac{d}{d\ln \mu } \bigg (\alpha _0 Z_\alpha (\alpha ,\lambda ,Y,\ln \Lambda /\mu )\bigg )\bigg |_{\alpha _0,\lambda _0,Y_0 = \mathrm{const}} \nonumber \\= & {} \alpha _0 Z_\alpha \bigg (\frac{\partial \ln Z_\alpha }{\partial \ln \mu } + \frac{\partial \ln Z_\alpha }{\partial \alpha }\, \frac{d\alpha }{d\ln \mu } \nonumber \\&+ \frac{\partial \ln Z_\alpha }{\partial \lambda ^{ijk}}\, \frac{d\lambda ^{ijk}}{d\ln \mu } + \frac{\partial \ln Z_\alpha }{\partial \lambda ^*_{ijk}} \nonumber \\&\qquad \times \frac{d\lambda ^*_{ijk}}{d\ln \mu } + \frac{\partial \ln Z_\alpha }{\partial Y}\, \frac{dY}{d\ln \mu }\bigg ). \end{aligned}$$
(125)

Here the derivative \(d/d\ln \mu \) acts on both \(\ln \Lambda /\mu \) inside \(\alpha \), \(\lambda \), Y and the explicitly written \(\ln \Lambda /\mu \). The partial derivative \(\partial /\partial \ln \mu \), by contrast, acts only on the explicitly written \(\ln \Lambda /\mu \) and does not act on \(\ln \Lambda /\mu \) inside \(\alpha \), \(\lambda \), and Y. Therefore,

$$\begin{aligned} \frac{\partial \ln Z_\alpha }{\partial \ln \mu }= & {} -\frac{\partial \ln Z_\alpha }{\partial \ln \Lambda } \equiv -\frac{d\ln Z_\alpha }{d\ln \Lambda }\bigg |_{\alpha ,\lambda ,Y=\mathrm{const}} \nonumber \\= & {} \frac{d\ln (Z_\alpha )^{-1}}{d\ln \Lambda }\bigg |_{\alpha ,\lambda ,Y=\mathrm{const}}. \end{aligned}$$
(126)

Next, we consider Eq. (125) in the point \(\ln \Lambda /\mu = x_0\). Due to the boundary conditions (122) \(\ln Z_\alpha (\alpha ,\lambda ,Y,\ln \Lambda /\mu \rightarrow x_0) = 0\) for all values of \(\alpha \), \(\lambda \), and Y. This implies that

$$\begin{aligned}&\frac{\partial \ln Z_\alpha }{\partial \alpha }\bigg |_{\ln \Lambda /\mu =x_0} = 0; \nonumber \\&\frac{\partial \ln Z_\alpha }{\partial \lambda ^{ijk}}\bigg |_{\ln \Lambda /\mu =x_0} = 0; \nonumber \\&\frac{\partial \ln Z_\alpha }{\partial \lambda ^*_{ijk}}\bigg |_{\ln \Lambda /\mu =x_0} = 0; \nonumber \\&\frac{\partial \ln Z_\alpha }{\partial Y}\bigg |_{\ln \Lambda /\mu =x_0} = 0. \end{aligned}$$
(127)

Taking into account that the left hand side of Eq. (125) does not depend on the value of \(\ln \Lambda /\mu \) at fixed values of the renormalized couplings due to the renormalization group equation, from Eqs. (125), (126), and (127) we obtain

$$\begin{aligned} {\widetilde{\beta }}(\alpha ,\lambda ,Y)= & {} \alpha \left. \frac{d\ln (Z_\alpha )^{-1}}{d\ln \Lambda }\right| _{\alpha ,\lambda ,Y=\mathrm{const}} \nonumber \\= & {} \left. \frac{d\alpha _0}{d\ln \Lambda }\right| _{\alpha ,\lambda ,Y=\mathrm{const}} = \beta (\alpha _0,\lambda _0,Y_0)\bigg |_{\ln \Lambda /\mu = x_0}.\nonumber \\ \end{aligned}$$
(128)

However, according to Eq. (122), in the point \(\ln \Lambda /\mu = x_0\) the values of the bare and renormalized couplings coincide,

$$\begin{aligned}&\alpha _0\bigg |_{\ln \Lambda /\mu = x_0} = \alpha ;\qquad \lambda _0^{ijk}\bigg |_{\ln \Lambda /\mu = x_0} = \lambda ^{ijk}; \nonumber \\&Y_0\bigg |_{\ln \Lambda /\mu = x_0} = Y. \end{aligned}$$
(129)

(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

$$\begin{aligned} {\widetilde{\beta }}(\alpha ,\lambda ,Y) = \beta (\alpha _0\rightarrow \alpha ,\lambda _0\rightarrow \lambda ,Y_0\rightarrow Y). \end{aligned}$$
(130)

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 \mu \) we can present the anomalous dimension of the quantum gauge superfield in the form

$$\begin{aligned} {\widetilde{\gamma }}_V(\alpha ,\lambda ,Y)= & {} \frac{\partial \ln Z_V}{\partial \ln \mu } + \frac{\partial \ln Z_V}{\partial \alpha }\, \frac{d\alpha }{d\ln \mu } \nonumber \\&+ \frac{\partial \ln Z_V}{\partial \lambda ^{ijk}} \frac{d\lambda ^{ijk}}{d\ln \mu } + \frac{\partial \ln Z_V}{\partial \lambda ^*_{ijk}} \nonumber \\&\times \frac{d\lambda ^*_{ijk}}{d\ln \mu } + \frac{\partial \ln Z_V}{\partial Y}\, \frac{dY}{d\ln \mu }. \end{aligned}$$
(131)

Again we consider this equation in the point \(\ln \Lambda /\mu = x_0\) and take into account that due to Eq. (122) \(\ln Z_V(\alpha ,\lambda ,Y,\ln \Lambda /\mu \rightarrow x_0) = 0\). Then, repeating the above argumentation, we obtain

$$\begin{aligned} {\widetilde{\gamma }}_V(\alpha ,\lambda ,Y)= & {} \left. - \frac{d\ln Z_V}{d\ln \Lambda }\right| _{\alpha ,\lambda ,Y=\mathrm{const}} \nonumber \\= & {} \gamma _V(\alpha _0\rightarrow \alpha ,\lambda _0\rightarrow \lambda ,Y_0\rightarrow Y). \end{aligned}$$
(132)

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 \(\alpha _0\rightarrow \alpha \), \(\lambda _0 \rightarrow \lambda \), \(Y_0\rightarrow 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

$$\begin{aligned}&\frac{{\widetilde{\beta }}(\alpha ,\lambda ,Y)}{\alpha ^2} \nonumber \\&\quad = - \frac{1}{2\pi }\bigg (3 C_2 - T(R) - 2C_2 {\widetilde{\gamma }}_c(\alpha ,\lambda ,Y) \nonumber \\&\qquad - 2C_2 {\widetilde{\gamma }}_V(\alpha ,\lambda ,Y) \nonumber \\&\qquad + \frac{1}{r} C(R)_i{}^j \left( {\widetilde{\gamma }}_\phi \right) _j{}^i(\alpha ,\lambda ,Y)\bigg );\nonumber \\&{\widetilde{\beta }}(\alpha ,\lambda ,Y) \nonumber \\&\quad = - \frac{\alpha ^2\bigg (3 C_2 - T(R) + C(R)_i{}^j \left( {\widetilde{\gamma }}_\phi \right) _j{}^i(\alpha ,\lambda ,Y)/r\bigg )}{2\pi (1- C_2\alpha /2\pi )}.\nonumber \\ \end{aligned}$$
(133)

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 \(\mu \). Certainly, in these schemes the NSVZ relation is also valid in all loops.Footnote 12 These statements agree with the explicit calculations (of some scheme dependent terms in RGFs) made in Refs. [50, 51].

Let us also note that the HD+MSL scheme can supplement various versions of the higher covariant derivative regularization, which differ in a particular form of the regulator functions R, F, and K and values of the parameters a and \(a_\varphi \). This implies that we obtain a set of in general different NSVZ renormalization prescription, which can certainly be related by finite renormalizations.

7 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 Sect. 4.1 the two-loop \(\beta \)-function for a general \({{\mathcal {N}}}=1\) supersymmetric gauge theory with a simple gauge group has been calculated in Ref. [74]. Also using this method the parts of the three-loop \(\beta \)-function containing the Yukawa couplings and ghost loops have been found in Refs. [63] and [52], respectively. Note that before this the part of the three-loop \(\beta \)-function containing the Yukawa couplings has also been calculated with the help of the standard technique in Refs. [50, 51]. For \({{\mathcal {N}}}=1\) SQED with \(N_f\) flavors the complete three-loop \(\beta \)-function in a general \(\xi \)-gauge has been calculated by the algorithm of Sect. 4.1 in Ref. [93]. Here all these calculations are used for checking the exact results described in the previous sections. Note that we will not verify that the Eqs. (28), (120), and (133) really hold, because this has already been done in Refs. [50,51,52, 74, 93]. 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.

7.1 The two-loop approximation

First, we reanalyse the two-loop calculation made in Ref. [74]. The corresponding contribution to the \(\beta \)-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 \({\varvec{V}}\)-legs in all possible ways. The method of Ref. [63] 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.

Fig. 3
figure 3

Vacuum supergraphs generating the two-loop contribution to the \(\beta \)-function

Let us start with the supergraph B1. The expression for it has been calculated in Ref. [63]. If we include the factor \(-2\pi /(r{{\mathcal {V}}}_4)\cdot d/d\ln \Lambda \), the result can be written as

$$\begin{aligned} \begin{aligned} \text {B1}\,&\rightarrow \ -\frac{2\pi }{r{{\mathcal {V}}}_4} \frac{d}{d\ln \Lambda } \cdot \frac{2}{3} {{\mathcal {V}}}_4 \\&\quad \times \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{\lambda _0^{ijk} \lambda ^*_{0ijk}}{Q^2 F_Q K^2 F_K (Q+K)^2 F_{Q+K}}, \end{aligned} \end{aligned}$$
(134)

where \(F_K\equiv F(K^2/\Lambda ^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 \(\delta \)-symbol \(\delta _m^n\)) by \(4\pi ^2 C(R)_m{}^n \delta ^4(P)\). This gives the contribution to the function (34) of the form

$$\begin{aligned} \Delta _{\mathrm{Yukawa}}\bigg (\frac{\beta }{\alpha _0^2}\bigg )= & {} -\frac{4\pi }{r} C(R)_i{}^m \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \nonumber \\&\times 4\pi ^2 \delta ^4(Q)\frac{\lambda _0^{ijk} \lambda ^*_{0mjk}}{F_Q K^2 F_K (Q+K)^2 F_{Q+K}}\qquad \nonumber \\= & {} -\frac{1}{\pi r} C(R)_i{}^j \frac{d}{d\ln \Lambda } \nonumber \\&\times \int \frac{d^4K}{(2\pi )^4} \frac{\lambda _0^{imn} \lambda ^*_{0jmn}}{K^4 F_K^2}. \end{aligned}$$
(135)

This result exactly agrees with the one obtained in Ref. [50] by the direct calculation of two-point superdiagrams with two external \({\varvec{V}}\)-legs.

The supergraphs B2 and B3 containing a matter loop produce two different contributions. If a matter loop corresponds to the superfields \(\phi _i\) and to the Pauli–Villars superfields \(\Phi _i\), then (for an arbitrary value of \(\xi _0\)) the result for the (properly modified) vacuum supergraphs multiplied by the operator (38) is given by

$$\begin{aligned}&\text {B2} + \text {B3}\ \rightarrow \ \frac{4\pi }{r}\, \text {tr}\, C(R)\, \frac{d}{d\ln \Lambda } \nonumber \\&\quad \times \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{e_0^2}{K^2 R_K} \nonumber \\&\quad \times \bigg \{\frac{1}{2Q^2(Q+K)^2} + \frac{1}{((K+Q)^2-Q^2)}\nonumber \\&\quad \times \bigg [\frac{F_{K+Q}^2}{2((K+Q)^2 F_{K+Q}^2 + M^2)} - \frac{F_Q^2}{2(Q^2 F_Q^2 + M^2)} \nonumber \\&\quad - \frac{M^2 F_{K+Q}'}{\Lambda ^2 F_{K+Q} ((K+Q)^2 F_{K+Q}^2 + M^2)}\nonumber \\&\quad + \frac{M^2 F_Q'}{\Lambda ^2 F_Q (Q^2 F_Q^2 + M^2)}\bigg ] \bigg \}, \end{aligned}$$
(136)

where \(R_K\equiv R(K^2/\Lambda ^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 \(\xi _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\pi ^2 C_2 \delta ^4(K)\) or by replacing \(\delta _i^j/Q^2\) or \(\delta _i^j/(Q+K)^2\) (coming from the propagators of the superfields \(\phi _i\)) by \(4\pi ^2 C(R)_i{}^j \delta ^4(Q)\) or \(4\pi ^2 C(R)_i{}^j \delta ^4(Q+K)\), respectively. Note that no replacement should be made for the non-singular Pauli–Villars propagators. The above procedure gives the result

$$\begin{aligned}&\Delta _{\mathrm{matter}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) \nonumber \\&\quad = \frac{64\pi ^4}{r}\,\text {tr}\left( C(R)^2\right) \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{\alpha _0}{K^2 R_K}\nonumber \\&\qquad \times \bigg \{\frac{1}{2(Q+K)^2}\,\delta ^4(Q) + \frac{1}{2Q^2}\, \delta ^4(Q+K) \bigg \}\nonumber \\&\qquad + \frac{64\pi ^4}{r}\,C_2\, \text {tr}\,C(R)\, \frac{d}{d\ln \Lambda } \nonumber \\&\qquad \times \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4}\, \delta ^4(K) \frac{\alpha _0}{R_K} \nonumber \\&\qquad \times \bigg \{\frac{1}{2Q^2(Q+K)^2} + \frac{1}{((K+Q)^2-Q^2)} \nonumber \\&\quad \qquad \times \bigg [\frac{F_{K+Q}^2}{2((K+Q)^2 F_{K+Q}^2 + M^2)} - \frac{F_Q^2}{2(Q^2 F_Q^2 + M^2)} \nonumber \\&\qquad \qquad - \frac{M^2 F_{K+Q}'}{\Lambda ^2 F_{K+Q} ((K+Q)^2 F_{K+Q}^2 + M^2)} \nonumber \\&\qquad \qquad + \frac{M^2 F_Q'}{\Lambda ^2 F_Q (Q^2 F_Q^2 + M^2)}\bigg ] \bigg \}. \end{aligned}$$
(137)

After calculating the integrals of the \(\delta \)-functions and some transformations we obtain the expression

$$\begin{aligned} \Delta _{\mathrm{matter}}\bigg (\frac{\beta }{\alpha _0^2}\bigg )= & {} \frac{4}{r}\,\text {tr}\left( C(R)^2\right) \frac{d}{d\ln \Lambda } \int \frac{d^4K}{(2\pi )^4} \frac{\alpha _0}{K^4 R_K} \qquad \nonumber \\&+ \frac{1}{2r} \,C_2\, \text {tr}\,C(R)\,\frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{\partial ^2}{\partial Q_\mu ^2} \nonumber \\&\times \bigg [\frac{\alpha _0 }{Q^2} \ln \bigg (1+\frac{M^2}{Q^2 F_Q^2}\bigg )\bigg ], \end{aligned}$$
(138)

which exactly coincides with the one found in Ref. [74] 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 \(\varphi _a\). In this case the calculation of the vacuum superdiagrams with an insertion of \(\theta ^4 (v^B)^2\) multiplied by the operator (38) for an arbitrary value of \(\xi _0\) gives

$$\begin{aligned}&\text {B2} + \text {B3}\ \rightarrow \ 4\pi C_2 \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{e_0^2}{K^2 R_K} \nonumber \\&\quad \times \bigg \{\frac{1}{(Q^2+M_\varphi ^2)((Q+K)^2+M_\varphi ^2)} - \frac{1}{((K+Q)^2-Q^2)} \nonumber \\&\qquad \times \bigg [\frac{R_{K+Q}^2}{2((K+Q)^2 R_{K+Q}^2 + M_\varphi ^2)} - \frac{R_Q^2}{2(Q^2 R_Q^2 + M^2)}\nonumber \\&\qquad - \frac{M_\varphi ^2 R_{K+Q}'}{\Lambda ^2 R_{K+Q} ((K+Q)^2 R_{K+Q}^2 + M_\varphi ^2)} \nonumber \\&\qquad + \frac{M_\varphi ^2 R_Q'}{\Lambda ^2 R_Q (Q^2 R_Q^2 + M_\varphi ^2)} + \frac{R_{K+Q}'}{\Lambda ^2 R_{K+Q}} - \frac{R_Q'}{\Lambda ^2 R_Q} \bigg ] \bigg \}. \end{aligned}$$
(139)

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\ \rightarrow \ 4\pi ^2 C_2 \delta ^4(K)\), after which the contribution of the Pauli–Villars superfields \(\varphi _a\) to the function (34) takes the form

$$\begin{aligned} \Delta _\varphi \bigg (\frac{\beta }{\alpha _0^2}\bigg )= & {} C_2^2 \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \nonumber \\&\times \bigg \{\frac{6\alpha _0}{Q^4} - \frac{\partial ^2}{\partial Q_\mu ^2}\bigg [\frac{\alpha _0}{Q^2}\ln \bigg (1+ \frac{M_\varphi ^2}{Q^2}\bigg )\nonumber \\&\quad - \frac{\alpha _0}{2Q^2} \ln \bigg (1+\frac{M_\varphi ^2}{Q^2 R_Q^2}\bigg ) - \frac{\alpha _0}{Q^2} \ln R_Q \bigg ]\bigg \}.\nonumber \\ \end{aligned}$$
(140)

This expression again coincides with the one found in Ref. [74] with the help of a different technique.Footnote 13

The contribution of the vacuum supergraphs B4, B5, B6, and B7 is written in the form

$$\begin{aligned}&\text {B4} + \text {B5} + \text {B6} + \text {B7}\ \rightarrow \ 4\pi C_2\, \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{e_0^2}{R_K R_Q} \nonumber \\&\quad \times \bigg \{-\frac{R_K}{2 Q^2 K^2 (K+Q)^2} - \frac{1}{2Q^2 K^2} \bigg (\frac{R_Q-R_K}{Q^2-K^2}\bigg ) \nonumber \\&\qquad - \frac{1}{R_{K+Q} K^2} \bigg (1-\frac{Q^2}{2(K+Q)^2}\bigg ) \bigg (\frac{R_Q-R_K}{Q^2-K^2}\bigg ) \nonumber \\&\qquad \times \bigg (\frac{R_{K+Q}-R_Q}{(K+Q)^2-Q^2}\bigg ) - \frac{1}{R_{K+Q} (K+Q)^2}\nonumber \\&\qquad \times \bigg (\frac{R_Q-R_K}{Q^2-K^2}\bigg )^2 +\frac{2}{K^2 \left( (K+Q)^2-Q^2\right) ^2} \nonumber \\&\qquad \times \bigg [ R_{K+Q} - R_Q - R_Q' \bigg (\frac{(K+Q)^2}{\Lambda ^2} - \frac{Q^2}{\Lambda ^2}\bigg )\bigg ] \nonumber \\&\qquad - \frac{Q_\mu K^\mu }{Q^2 K^2}\nonumber \\&\qquad \times \bigg [\frac{R_{K+Q}}{\left( (K+Q)^2-K^2\right) \left( (K+Q)^2-Q^2\right) } \nonumber \\&\qquad \quad + \frac{R_{K}}{\left( K^2-(K+Q)^2\right) \left( K^2-Q^2\right) } \nonumber \\&\qquad \quad + \frac{R_{Q}}{\left( Q^2-(K+Q)^2\right) \left( Q^2-K^2\right) }\bigg ] \bigg \}. \end{aligned}$$
(141)
Fig. 4
figure 4

The three-loop vacuum supergraphs containing vertices with the Yukawa couplings

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 Sect. 5. Also all terms dependent on the gauge parameter \(\xi _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 \(\delta \)-function multiplied by \(4\pi ^2 C_2\) and sum up all expression thus obtained. The result is

$$\begin{aligned}&\Delta _{\mathrm{gauge+ghost}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) \nonumber \\&\quad = 64\pi ^4 C_2^2\, \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{\alpha _0}{R_K R_Q}\nonumber \\&\qquad \times \bigg \{-\frac{R_K}{2 K^2 (K+Q)^2}\, \delta ^4(Q) \nonumber \\&\qquad \quad -\frac{R_K}{2 Q^2 (K+Q)^2}\, \delta ^4(K) -\frac{R_K}{2 Q^2 K^2}\, \delta ^4(K+Q) \nonumber \\&\qquad \quad - \bigg (\frac{1}{2 K^2}\,\delta ^4(Q) + \frac{1}{2Q^2} \delta ^4(K) \bigg )\bigg (\frac{R_Q-R_K}{Q^2-K^2}\bigg ) - \frac{1}{R_{K+Q}} \nonumber \\&\quad \qquad \times \bigg [\,\delta ^4(K) \bigg (1-\frac{Q^2}{2(K+Q)^2}\bigg ) - \frac{Q^2}{2 K^2}\delta ^4(K+Q) \bigg ] \nonumber \\&\quad \qquad \times \bigg (\frac{R_Q-R_K}{Q^2-K^2}\bigg ) \bigg (\frac{R_{K+Q}-R_Q}{(K+Q)^2-Q^2}\bigg )\nonumber \\&\qquad \quad - \frac{1}{R_{K+Q}}\, \delta ^4(K+Q) \bigg (\frac{R_Q-R_K}{Q^2-K^2}\bigg )^2 +\frac{2}{\left( (K+Q)^2-Q^2\right) ^2}\,\delta ^4(K) \nonumber \\&\quad \qquad \times \bigg [ - R_Q' \bigg (\frac{(K+Q)^2}{\Lambda ^2} - \frac{Q^2}{\Lambda ^2}\bigg ) +R_{K+Q} - R_Q \bigg ] \nonumber \\&\quad \qquad - \bigg (\frac{Q_\mu K^\mu }{K^2}\,\delta ^4(Q) + \frac{Q_\mu K^\mu }{Q^2}\,\delta ^4(K) \bigg ) \nonumber \\&\quad \qquad \times \bigg [\frac{R_{K+Q}}{\left( (K+Q)^2-K^2\right) \left( (K+Q)^2-Q^2\right) } \nonumber \\&\quad \qquad + \frac{R_{K}}{\left( K^2-(K+Q)^2\right) \left( K^2-Q^2\right) } \nonumber \\&\quad \qquad + \frac{R_{Q}}{\left( Q^2-(K+Q)^2\right) \left( Q^2-K^2\right) }\bigg ] \bigg \}. \end{aligned}$$
(142)

After calculating the integrals of the \(\delta \)-functions this expression can be rewritten as

$$\begin{aligned}&\Delta _{\mathrm{gauge+ghost}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) \nonumber \\&\quad = 4 C_2^2\, \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \nonumber \\&\qquad \times \bigg \{-\frac{3\alpha _0}{2 Q^4} - \frac{\alpha _0}{\Lambda ^4 R_Q^2} (R_Q')^2 + \frac{\alpha _0}{\Lambda ^4 R_Q} R_Q''\bigg \}\qquad \nonumber \\&\quad = C_2^2 \frac{d}{d\ln \Lambda }\int \frac{d^4Q}{(2\pi )^4}\bigg \{-\frac{6\alpha _0}{Q^4} + \frac{\partial ^2}{\partial Q_\mu ^2} \bigg [\frac{\alpha _0}{Q^2}\ln R_Q\bigg ]\bigg \}.\nonumber \\ \end{aligned}$$
(143)

Again we see that the result coincides with the one obtained by the method of Ref. [63] in Ref. [74]. This confirms the correctness of the argumentation used in this paper for deriving the NSVZ equation.

The overall two-loop expression for the \(\beta \)-function reads as

$$\begin{aligned}&\frac{\beta (\alpha _0,\lambda _0,Y_0)}{\alpha _0^2} - \frac{\beta _{\mathrm{1-loop}}(\alpha _0)}{\alpha _0^2} \nonumber \\&\quad = \int \frac{d^4Q}{(2\pi )^4} \frac{d}{d\ln \Lambda } \nonumber \\&\qquad \times \bigg \{ - \alpha _0 C_2^2\, \frac{\partial ^2}{\partial Q_\mu ^2} \nonumber \\&\qquad \quad \times \bigg [\frac{1}{2Q^2}\ln \bigg (1+\frac{M_\varphi ^2}{Q^2 R_Q^2}\bigg ) - \frac{1}{Q^2} \ln \bigg (1+\frac{M_\varphi ^2}{Q^2}\bigg )\bigg ] \nonumber \\&\quad \qquad + \frac{\alpha _0}{2r}\, C_2\,\text {tr}\, C(R)\, \frac{\partial ^2}{\partial Q_\mu ^2} \nonumber \\&\quad \qquad \times \bigg [\frac{1}{Q^2}\ln \bigg (1+\frac{M^2}{Q^2 F_Q^2}\bigg )\bigg ] + \frac{4\alpha _0}{r}\, \text {tr}\left( C(R)^2\right) \nonumber \\&\quad \qquad \times \frac{1}{Q^4 R_Q} - \frac{1}{\pi r} \lambda _0^{imn} \lambda ^*_{0jmn} C(R)_i{}^j \frac{1}{Q^2 F_Q^2}\bigg \} \nonumber \\&\quad + O(\alpha _0^2,\alpha _0\lambda _0^2,\lambda _0^4) \end{aligned}$$
(144)

and coincides with the one obtained in Ref. [74]. Certainly, this expression satisfies the NSVZ equations (28) and (120) and agrees with the previous calculations (first made in Ref. [97]).

7.2 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 \(\beta \)-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 \({\varvec{V}}\) has been made in Refs. [50, 51] in the Feynman gauge \(\xi _0=1\). Subsequently, the result has been reobtained by the method described in Sect. 4.1 in Ref. [63]. 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

$$\begin{aligned}&B8\ \rightarrow \ \frac{4\pi }{r} \frac{d}{d\ln \Lambda } \nonumber \\&\quad \times \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \lambda _0^{ijk} \lambda ^*_{0ijl} \lambda _0^{mnl} \lambda ^*_{0mnk} \nonumber \\&\quad \times \frac{1}{Q^2 F_Q (K+Q)^2 F_{K+Q} L^2 F_L (K+L)^2 F_{K+L} K^2 F_K^2},\nonumber \\ \end{aligned}$$
(145)

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 \(\delta \)-symbols coming from the corresponding propagator by momentum \(\delta \)-functions multiplied by \(4\pi ^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 \(\delta \)-functions we obtain the required contribution

$$\begin{aligned} \Delta _{\mathrm{B8}}\bigg (\frac{\beta }{\alpha _0^2}\bigg )= & {} \frac{1}{\pi r} C(R)_p{}^l \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \nonumber \\&\times \lambda _0^{ijk} \lambda ^*_{0ijl} \lambda _0^{mnp} \lambda ^*_{0mnk} \frac{1}{Q^4 F_Q^2 L^4 F_L^2}\nonumber \\&+ \frac{4}{\pi r} C(R)_i{}^p \frac{d}{d\ln \Lambda } \int \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \nonumber \\&\times \lambda _0^{ijk} \lambda ^*_{0pjl} \lambda _0^{mnl} \lambda ^*_{0mnk} \frac{1}{K^4 F_K^3 L^2 F_L (K+L)^2 F_{K+L}},\nonumber \\ \end{aligned}$$
(146)

which exactly coincides with the corresponding result of Refs. [50, 63] found by different methods.

According to [63] the properly modified vacuum supergraph B9 is given by the expression

$$\begin{aligned}&B9\ \rightarrow \ -\frac{8\pi }{r} \frac{d}{d\ln \Lambda } \nonumber \\&\quad \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4}\, e_0^2 \lambda _0^{ijk} \lambda ^*_{0imn} (T^A)_j{}^m (T^A)_k{}^n\nonumber \\&\quad \times \frac{N(Q,K,L)}{K^2 R_K L^2 F_L Q^2 F_Q (Q+K)^2 F_{Q+K} (Q-L)^2 F_{Q-L}}\nonumber \\&\quad \times \frac{1}{(Q+K-L)^2 F_{Q+K-L}}, \end{aligned}$$
(147)

where \(K_\mu \) corresponds to the propagator of the quantum gauge superfield and, following Ref. [51], we use the notation

$$\begin{aligned} N(Q,K,L)\equiv & {} L^2 F_{Q+K} F_{Q+K-L} - Q^2 \bigg ((Q+K)^2 - L^2\bigg ) \nonumber \\&\times F_{Q+K-L} \frac{F_{Q+K}-F_Q}{(Q+K)^2-Q^2}\nonumber \\&- (Q-L)^2 \bigg ((Q+K-L)^2-L^2\bigg ) F_{Q+K} \nonumber \\&\times \frac{F_{Q+K-L}-F_{Q-L}}{(Q+K-L)^2-(Q-L)^2} + Q^2 (Q-L)^2 \nonumber \\&\times \bigg (L^2 - (Q+K)^2 - (Q+K-L)^2\bigg ) \nonumber \\&\times \bigg (\frac{F_{Q+K}-F_Q}{(Q+K)^2-Q^2}\bigg ) \nonumber \\&\times \bigg (\frac{F_{Q+K-L}-F_{Q-L}}{(Q+K-L)^2-(Q-L)^2}\bigg ). \end{aligned}$$
(148)

Again, Eq. (147) does not contain bad terms. Next, we proceed according to the algorithm of Sect. 4.2. In this case the relevant replacements are \(1/K^2 \rightarrow 4\pi ^2 C_2 \delta ^4(K)\) and \(\delta _i^j/P^2 \rightarrow 4\pi ^2 C(R)_i{}^j \delta ^4(P)\), where \(P_\mu \) stands for momenta of the matter superfields, namely, \(Q_\mu \), \((Q+K)_\mu \), \((Q-L)_\mu \), and \((Q+K-L)_\mu \). After some transformations involving the identities

$$\begin{aligned}&\lambda _0^{ijk} \lambda ^*_{0imn} (T^A)_j{}^m (T^A)_k{}^n = -\frac{1}{2} \lambda _0^{imn} \lambda ^*_{0jmn} C(R)_i{}^j; \end{aligned}$$
(149)
$$\begin{aligned}&\lambda _0^{ijk} \lambda ^*_{0imn} C(R)_j{}^l (T^A)_l{}^m (T^A)_k{}^n = -\frac{1}{2} \lambda _0^{imn} \lambda ^*_{0jmn} \left( C(R)^2\right) _i{}^j;\nonumber \\ \end{aligned}$$
(150)
$$\begin{aligned}&\lambda _0^{ijk} \lambda ^*_{0lmn} C(R)_i{}^l (T^A)_j{}^m (T^A)_k{}^n \nonumber \\&\quad = \frac{1}{2} \lambda _0^{imn} \lambda ^*_{0jmn} \left( C(R)^2\right) _i{}^j - \lambda _0^{ijm} \lambda ^*_{0klm} C(R)_i{}^k C(R)_j{}^l,\nonumber \\ \end{aligned}$$
(151)

which follow from Eq. (10), we get the result

$$\begin{aligned}&\Delta _{\mathrm{B9}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) = \frac{1}{\pi r}\, \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \nonumber \\&\quad \times e_0^2 \lambda _0^{imn}\lambda ^*_{0jmn} C_2 C(R)_i{}^j \frac{N(Q,0,L)}{Q^4 F_Q^2 L^2 F_L (Q-L)^4 F_{Q-L}^2}\nonumber \\&\quad - \frac{1}{\pi r}\, \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4}\, e_0^2 \nonumber \\&\quad \times \bigg (\lambda _0^{imn}\lambda ^*_{0jmn} \left( C(R)^2\right) _i{}^j - 2 \lambda _0^{ijm} \lambda ^*_{0klm} C(R)_i{}^k C(R)_j{}^l\bigg ) \nonumber \\&\quad \times \frac{N(Q,K,0)}{K^2 R_K~Q^4 F_Q^2 (Q+K)^4 F_{Q+K}^2} \nonumber \\&\quad + \frac{4}{\pi r}\, \frac{d}{d\ln \Lambda } \int \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \nonumber \\&\quad \times e_0^2 \lambda _0^{imn}\lambda ^*_{0jmn} \left( C(R)^2\right) _i{}^j \frac{1}{K^4 R_K F_K L^4} \nonumber \\&\quad \times \frac{N(0,K,L)}{F_L^2 (K-L)^2 F_{K-L}}, \end{aligned}$$
(152)

which agrees with the one found in [51, 63] by different methods.

The supergraph B10 is given by the expression

$$\begin{aligned}&B10\ \rightarrow \ -\frac{8\pi }{r} \frac{d}{d\ln \Lambda } \nonumber \\&\quad \times \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4}\, e_0^2 \lambda _0^{ijk} \lambda ^*_{0ijl} (T^A)_k{}^m (T^A)_m{}^l\nonumber \\&\quad \times \frac{L(Q,Q+K)}{K^2 R_K Q^2 F_Q^2 (Q+L)^2 F_{Q+L} (Q+K)^2 F_{Q+K} L^2 F_L},\nonumber \\ \end{aligned}$$
(153)

where \(K_\mu \) denotes the momentum of the quantum gauge superfield propagator and, again following Ref. [51],

$$\begin{aligned} L(Q,P)\equiv & {} F_Q F_P + \frac{F_P-F_Q}{P^2-Q^2}\, \bigg (F_Q Q^2 + F_P P^2\bigg ) \nonumber \\&+ 2 Q^2 P^2 \bigg (\frac{F_P-F_Q}{P^2-Q^2}\bigg )^2. \end{aligned}$$
(154)

In this case the algorithm of Sect. 4.2 produces the expression

$$\begin{aligned}&\Delta _{\mathrm{B10}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) = -\frac{2}{\pi r} \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \nonumber \\&\qquad \times e_0^2\, \lambda _0^{imn} \lambda ^*_{0jmn} C_2 C(R)_i{}^j \frac{L(Q,Q)}{Q^4 F_Q^3 L^2 F_L (Q+L)^2 F_{Q+L}} \nonumber \\&\qquad -\frac{2}{\pi r} \frac{d}{d\ln \Lambda } \int \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \nonumber \\&\qquad \times e_0^2\, \lambda _0^{imn} \lambda ^*_{0jmn} \left( C(R)^2\right) _i{}^j \nonumber \\&\qquad \times \bigg (\frac{L(0,K)}{K^4 R_K F_K L^4 F_L^2} + \frac{L(K,0)}{K^4 R_K F_K^2 L^2 F_L} \nonumber \\&\qquad \times \frac{1}{(K+L)^2 F_{K+L}}\bigg ) - \frac{4}{\pi r} \frac{d}{d\ln \Lambda } \nonumber \\&\qquad \times \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4}\, e_0^2\, \lambda _0^{ijm} \lambda ^*_{0klm} C(R)_i{}^k C(R)_j{}^l \nonumber \\&\qquad \frac{L(Q,Q+K)}{K^2 R_K Q^4 F_Q^3~(Q+K)^2 F_{Q+K}}, \end{aligned}$$
(155)

which also coincides with the results of Refs. [51, 63].

The expression for the last supergraph B11 reads as

$$\begin{aligned}&B11\ \rightarrow \ \frac{8\pi }{r} \frac{d}{d\ln \Lambda } \nonumber \\&\quad \times \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4}\, e_0^2 \lambda _0^{ijk} \lambda ^*_{0ijl} (T^A)_k{}^m (T^A)_m{}^l\nonumber \\&\quad \times \frac{K(Q,K)}{K^2 R_K Q^2 F_Q^2 L^2 F_L (Q+L)^2 F_{Q+L}}, \end{aligned}$$
(156)

where the momentum of the gauge superfield propagator is denoted by \(K_\mu \) and

$$\begin{aligned} K(Q,K)\equiv & {} \frac{F_{Q+K}-F_Q-2Q^2 F_Q'/\Lambda ^2}{(Q+K)^2 - Q^2} \nonumber \\&+ \frac{2Q^2 (F_{Q+K}-F_Q)}{((Q+K)^2-Q^2)^2}. \end{aligned}$$
(157)

Replacing the squared inverse momenta and the corresponding \(\delta \)-symbols according to the prescription given in Sect. 4.2 we obtain the part of the function (34),

$$\begin{aligned}&\Delta _{\mathrm{B11}}\bigg (\frac{\beta }{\alpha _0^2}\bigg ) = \frac{2}{\pi r} \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \nonumber \\&\qquad \times e_0^2\, \lambda _0^{imn} \lambda ^*_{0jmn} C_2 C(R)_i{}^j \frac{K(Q,0)}{Q^2 F_Q^2 L^2 F_L (Q+L)^2 F_{Q+L}}\nonumber \\&\qquad + \frac{2}{\pi r} \frac{d}{d\ln \Lambda } \int \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \nonumber \\&\qquad \times e_0^2\, \lambda _0^{imn} \lambda ^*_{0jmn} \left( C(R)^2\right) _i{}^j \frac{K(0,K)}{K^2 R_K L^4 F_L^2}\nonumber \\&\qquad + \frac{4}{\pi r} \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \nonumber \\&\qquad \times e_0^2\, \lambda _0^{ijm} \lambda ^*_{0klm} C(R)_i{}^k C(R)_j{}^l \frac{K(Q,K)}{K^2 R_K Q^4 F_Q^3}. \end{aligned}$$
(158)

Again, it agrees with the calculations made previously.

7.3 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. [52]. These supergraphs are presented in Fig. 5. Note that in this approximation the first nonlinear term in the function \({{\mathcal {F}}}(V)\) (see Eq. (8)) is essential. It generates the vertex

$$\begin{aligned} - \frac{3}{4} e_0^2\, y_0\, G^{ABCD} \int d^8x\, ({\bar{c}}^A + {\bar{c}}^{+A}) V^C V^D (c^B - c^{+B})\nonumber \\ \end{aligned}$$
(159)

present in the supergraph B21, see Ref. [87] for details. This vertex is denoted by a cross.

Fig. 5
figure 5

The three-loop vacuum supergraphs containing a ghost loop (or ghost loops) and the corresponding superdiagrams contributing to the ghost and matter anomalous dimensions. The gray circles denote insertions of the one-loop polarization operator of the quantum gauge superfield, see Ref. [90] for details

Expressions for the supergraphs in Fig. 5 contain \(1/K^4\) singularities, where \(K_\mu \) 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 \(\gamma _c(\alpha _0,\lambda _0,Y_0)\) and \((\gamma _\phi )_i{}^j(\alpha _0,\lambda _0,Y_0)\). With the help of the method described in Sect. 4.1 this has been done in Ref. [52]. Also it is possible to use these results for checking the modification of the algorithm discussed in Sect. 4.2. This is made as follows:

  1. 1.

    First, we consider a vacuum supergraph containing a ghost loop (or loops) with an insertion of \(\theta ^4 (v^B)^2\) and calculate it using the D-algebra.

  2. 2.

    Next, we replace one of inverse squared momenta (say, \(1/Q^2\)) of a ghost propagator by \(4\pi ^2 C_2 \delta ^4(Q)\). If matter propagators are present in the considered supergraph, then \(\delta _i^j/Q^2\) coming from a matter propagator is replaced by \(4\pi ^2 C(R)_i{}^j \delta ^4(Q)\). All expressions obtained after the replacements of one propagator should be summed.

  3. 3.

    After calculating the integrals of the \(\delta \)-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

$$\begin{aligned}&\Delta _{\mathrm{B}i} \bigg (\frac{\beta }{\alpha _0^2}\bigg )\ \rightarrow \ \frac{C_2}{\pi } \Delta _{\mathrm{B}i} \gamma _c,\qquad \text {for}\qquad i=12,\ldots ,21; \nonumber \\ \end{aligned}$$
(160)
$$\begin{aligned}&\Delta _{\mathrm{B}22} \bigg (\frac{\beta }{\alpha _0^2}\bigg ) + \Delta _{\mathrm{B}23} \bigg (\frac{\beta }{\alpha _0^2}\bigg ) \nonumber \\&\quad \rightarrow \ \frac{C_2}{\pi } \Delta _{\mathrm{B}22 + \mathrm{B}23} \gamma _c - \frac{1}{2\pi r} C(R)_i{}^j (\Delta _{\mathrm{B}22 + \mathrm{B}23}\gamma _\phi )_j{}^i,\nonumber \\ \end{aligned}$$
(161)

in agreement with Eqs. (87) and (98). Note that this verification is different from the one made in Ref. [52], because the algorithm used for constructing the contribution to the function (34) is different.

7.4 \({{\mathcal {N}}}=1\) SQED in the three-loop approximation

As one more example for checking the method proposed in Sect. 4.2 we can consider \({{\mathcal {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 \rightarrow \delta _{\alpha \beta } \cdot 1_2\), where \(\alpha ,\beta =1,\ldots , N_f\) and \(1_2\) denotes an identity matrix of the size \(2\times 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 \({{\mathcal {N}}}=1\) SQED is a particular case of \({{\mathcal {N}}}=1\) supersymmetric gauge theories considered earlier, this example is not trivial, because for this theory the complete expression for the three-loop \(\beta \)-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.

Using the results of Ref. [93] the sum of two- and three-loop vacuum supergraphs modified by an insertion of \(\theta ^4 (v^B)^2\) and multiplied by \(-2\pi /{{\mathcal {V}}}_4\cdot d/d\ln \Lambda \) in the general \(\xi _0\)-gauge is written as

$$\begin{aligned}&4\pi N_f \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{e_0^2}{K^2 R_K} \nonumber \\&\quad \times \bigg (\frac{1}{Q^2 (Q+K)^2} - \frac{1}{(Q^2+M^2)((Q+K)^2+M^2)}\bigg )\nonumber \\&\quad - 4\pi N_f^2 \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \frac{e_0^4}{K^2 R_K^2} \nonumber \\&\quad \times \bigg (\frac{1}{Q^2 (Q+K)^2} - \frac{1}{(Q^2+M^2)((Q+K)^2+M^2)}\bigg )\nonumber \\&\quad \times \bigg (\frac{1}{L^2 (L+K)^2} - \frac{1}{(L^2+M^2)((L+K)^2+M^2)}\bigg ) \nonumber \\&\quad + 8\pi N_f \frac{d}{d\ln \Lambda } \int \frac{d^4Q}{(2\pi )^4} \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4}\frac{e_0^4}{K^2 R_K L^2 R_L} \nonumber \\&\quad \times \bigg \{\frac{1}{Q^2 (Q+K)^2 (Q+L)^2} - \frac{K^2}{Q^2 (Q+K)^2 (Q+L)^2 (Q+K+L)^2} \nonumber \\&\qquad + \frac{K^2+M^2}{(Q^2+M^2)((Q+K)^2+M^2)((Q+L)^2+M^2)((Q+K+L)^2+M^2)} \nonumber \\&\qquad - \frac{1}{(Q^2+M^2)((Q+K)^2+M^2)((Q+L)^2+M^2)}\nonumber \\&\qquad + \frac{2M^2}{(Q^2+M^2)^2((Q+K)^2+M^2)((Q+L)^2+M^2)}\bigg \}, \end{aligned}$$
(162)

and does not contain inverse momenta to the fourth power. According to the method discussed in Sect. 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\pi ^2 \delta ^4(P)\), where \(P_\mu \) 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_\mu \), while in the last term the momenta of the quantum gauge superfield propagators are \(K_\mu \) and \(L_\mu \). Then the above described procedure gives

$$\begin{aligned}&\frac{\beta (\alpha _0)}{\alpha _0^2} - \frac{\beta _{\mathrm{1-loop}}(\alpha _0)}{\alpha _0^2} \nonumber \\&\quad = \frac{2 N_f}{\pi } \frac{d}{d\ln \Lambda } \int \frac{d^4K}{(2\pi )^4} \frac{e_0^2}{K^4 R_K} \nonumber \\&\qquad - \frac{4N_f^2}{\pi } \frac{d}{d\ln \Lambda } \int \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \frac{e_0^4}{K^4 R_K^2}\nonumber \\&\qquad \times \bigg (\frac{1}{L^2 (L+K)^2} - \frac{1}{(L^2+M^2)((L+K)^2+M^2)}\bigg ) \nonumber \\&\qquad - \frac{N_f}{\pi } \frac{d}{d\ln \Lambda } \nonumber \\&\qquad \times \int \frac{d^4K}{(2\pi )^4} \frac{d^4L}{(2\pi )^4} \frac{e_0^4}{R_K R_L}\qquad \nonumber \\&\qquad \times \bigg (\frac{4}{K^2 L^4 (K+L)^2} - \frac{2}{K^4 L^4}\bigg ) + O(e_0^6). \end{aligned}$$
(163)

Again, this expression correctly reproduces the result obtained earlier by different methods, see Refs. [53, 93].

Therefore, for all considered supergraphs the results obtained with the help of the technique discussed in Sect. 4.2 (which is a certain modification of the one proposed in [63]) 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 \(\beta \)-function.

8 Conclusion

In this paper we have finished the all-order perturbative derivation of the exact NSVZ \(\beta \)-function (1) for non-Abelian \({{\mathcal {N}}}=1\) supersymmetric gauge theories which was started in Refs. [2, 63, 64]. 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. [64], with this regularization all loop integrals giving the \(\beta \)-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 \(\beta \)-function starting from the two-loop approximation. (The one-loop quantum corrections should be considered separately. This has been done in Ref. [78].) 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 \(\beta \)-function. The NSVZ equation in the form (5) relates this contribution to the \(\beta \)-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 \(\beta \)-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 for arbitrary higher-derivative regulator functions R, F, and K and arbitrary values of the regularization parameters \(a=M/\Lambda \) and \(a_\varphi =M_\varphi /\Lambda \) independent of couplings. 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. [64] 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. The original NSVZ \(\beta \)-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 an arbitrary version of 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 \Lambda /\mu \) (or, equivalently, set all finite constants to 0). Note that the HD+MSL schemes corresponding to different versions of the higher covariant derivative regularization constitute a continuous set of the NSVZ schemes, which are in general different.

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., [50,51,52, 93]. Also the results of this paper confirm the correctness of the expression for the three-loop \(\beta \)-function of a general \({{\mathcal {N}}}=1\) supersymmetric gauge theory with matter superfields and a simple gauge group derived in Ref. [94] from the NSVZ equation.