The three-loop anomalous dimension and the four-loop β-function for N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 1 SQED regularized by higher derivatives

For N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 1 SQED with Nf flavors regularized by higher derivatives in the general ξ-gauge we calculate the three-loop anomalous dimension of the matter superfields defined in terms of the bare coupling constant and demonstrate its gauge independence. After this the four-loop β-function defined in terms of the bare coupling constant is obtained with the help of the NSVZ equation, which is valid for these renormalization group functions in all loops. Next, we calculate the three-loop anomalous dimension and the four-loop β-function defined in terms of the renormalized coupling constant for an arbitrary subtraction scheme supplementing the higher derivative regularization. Then we construct a renormalization prescription for which the results coincide with the ones in the DR¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{DR}} $$\end{document}-scheme and describe all NSVZ schemes in the considered approximation. Also we demonstrate the existence of a subtraction scheme in which the anomalous dimension does not depend on Nf, while the β-function contains only terms of the first order in Nf. This scheme is obtained with the help of a finite renormalization compatible with a structure of quantum corrections and is NSVZ. The existence of such an NSVZ scheme is also proved in all loops.


Introduction
Calculations of higher order quantum corrections in N = 1 supersymmetric theories are important for both theory and phenomenological applications [1]. Most of these calculations (see, e.g., [2][3][4][5][6]) were made in the DR scheme, when a theory is regularized by dimensional reduction [7] and divergences are removed by modified minimal subtraction [8]. However, the dimensional reduction is not mathematically consistent [9] and, as a result, can JHEP04(2022)108 break supersymmetry in very higher loops [10][11][12]. A much better regularization in the supersymmetric case is the higher covariant derivative method [13][14][15]. In the superfield formulation [16,17] (see also [18,19]) it is mathematically consistent, does not break supersymmetry in all orders, and, most importantly, allows revealing some special features of quantum corrections, see, e.g., [20,21]. For instance, it was demonstrated that in the case of using the higher covariant derivative regularization the loop integrals giving the β-function are integrals of double total derivatives with respect to the loop momenta in all orders [22,23]. 1 This fact turns out to be the key to explaining the appearance of the NSVZ equation [35][36][37][38], which relates the β-function of N = 1 supersymmetric gauge theories to the anomalous dimension of the matter superfields. Its perturbative proof in the Abelian case has been done in [22] and [39] by two different methods. A much more complicated proof in the non-Abelian case was constructed in [40], see also [41], on the base of the non-renormalization theorem for the triple gauge ghost vertices [42]. 2 It turned out that the NSVZ equation is valid in all loops for the renormalization group functions (RGFs) defined in terms of the bare couplings in the case of using the higher covariant derivative regularization for an arbitrary renormalization prescription. For standard RGFs, which are defined in terms of the renormalized couplings, the NSVZ equation is valid in the HD+MSL scheme [46] when the theory is regularized by higher covariant derivatives and the renormalization is made with the help of minimal subtractions of logarithms [47,48]. In the Abelian case another all-loop NSVZ renormalization prescription is the on-shell scheme [49].
The higher covariant derivative regularization also allows deriving some NSVZ-like relations and constructing the corresponding all-loop NSVZ renormalization prescriptions. For example, the NSVZ-like relation for the Adler D-function in N = 1 SQCD [50,51] (see also [52]) can also be derived using the method proposed in [22]. In the softly broken SQED the NSVZ-like relation for the renormalization of the gaugino mass [53][54][55] can also be obtained by this method [56,57].
Nevertheless, although the higher covariant derivative regularization has a lot of attractive features, its application is technically difficult due to the complicated structure of vertices coming from the higher derivative terms and complicated loop integrals. At present it was used only for calculating the two-loop anomalous dimension and the three-loop βfunction, see, e.g., [58,59]. Even in this case the calculation of loop integrals was rather nontrivial, see, e.g., [24,60]. Therefore, the question arises as to whether it is possible to use the higher covariant derivative regularization for making explicit calculations in higher orders.
In this paper we consider the N = 1 supersymmetric electrodynamics (SQED) with N f flavors (which is the simplest N = 1 supersymmetric gauge theory) and calculate for it the three-loop anomalous dimension of the matter superfields and the four-loop β-function.

JHEP04(2022)108
In section 2 we briefly describe how this theory is formulated and regularized by higher derivatives. After this in section 3 we obtain the anomalous dimension defined in terms of the bare coupling constant. The corresponding supergraphs were generated and reduced to the momentum integrals with the help of a special C++ computer programm written by I.S. The resulting integrals are calculated with the help of the Chebyshev polynomial technique [61] in appendix A. After obtaining the result for the three-loop anomalous dimension defined in terms of the bare coupling constant we construct the expression for the four-loop β-function defined in terms of the bare coupling constant. For this purpose we use the general statement that RGFs defined in terms of the bare coupling constant satisfy the NSVZ equation in the case of using the higher covariant derivative regularization, see, e.g., [22]. Note that this statement is valid for an arbitrary renormalization prescription because RGFs defined in terms of the bare couplings are scheme independent for a fixed regularization, but, certainly, depend on a regularization. RGFs standardly defined in terms of the renormalized coupling constant are obtained in section 4. They depend on a regularization and a renormalization prescription and satisfy the NSVZ equation only in some special subtractions schemes (which in particular include the HD+MSL and on-shell schemes). Analysing the scheme dependence we see that all terms proportional to N f in the three-loop anomalous dimension and proportional to (N f ) 2 in the four-loop β-function can be set to 0 by a special choice of a renormalization prescription. In the resulting subtraction scheme (which we will call "the minimal scheme") the expressions for RGFs have the simplest form and satisfy the NSVZ equation. In section 5 we demonstrate that the minimal scheme exists in all orders. Note that the proof is also valid in the nonsupersymmetric case, say, for QED with N f flavors or for scalar QED with N f flavors. For N = 1 SQED we also prove that the minimal scheme is NSVZ in all orders.

N = 1 SQED with N f flavors
In N = 1 superspace the action of massless N = 1 SQED with N f flavors is written as where φ α and φ α are N f pairs of chiral matter superfields with opposite U(1) charges, and V is a real gauge superfield. In our notation the bare and renormalized coupling constants are denoted by e 0 and e, respectively. We regularize the theory (2.1) by the higher covariant derivatives following the ideas proposed in [13][14][15] and generalized to the supersymmetric case. Here we mainly use the same conventions as in [58], but, for completeness, briefly describe the details of the regularization and quantization. The regularization is introduced by adding the higher derivative terms to the action (2.1), after which the regularized action

JHEP04(2022)108
will include the regulator function R(x). This function grows at infinity and is equal to 1 at x = 0. Another similar regulator function appears in the gauge fixing term where ξ 0 is the bare gauge parameter. The minimal (Feynman) gauge corresponds to ξ 0 = 1 and R(x) = K(x). However, we will make calculations for an arbitrary ξ 0 and K(x) = R(x). The generating functional of the regularized theory also contains the Pauli-Villars determinant needed for removing residual one-loop divergences, Here the action for the massive Pauli-Villars superfields is given be the expression and it is important that the ratio of the Pauli-Villars mass M to the regularization parameter Λ should not depend on the coupling constant, Due to the renormalizability of the considered theory the ultraviolet divergences can be absorbed into the renormalization of the coupling constant, of the gauge parameter, and of the chiral matter superfields φ α and φ α . Note that all chiral superfields have the same renormalization constant Z, such that φ α, It is convenient to encode the ultraviolet divergences in RGFs. Following [46], we distinguish between RGFs defined in terms of the bare coupling constant, (2.8) and the ones standardly defined in terms of the renormalized coupling constant by the equations where µ is a renormalization point. RGFs (2.8) are independent of a renormalization prescription for a fixed regularization, but depend on a regularization. According to [22,39] in the case of using the higher derivative regularization described above they satisfy the NSVZ equation [62,63] β(α 0 )

JHEP04(2022)108
in all loops for an arbitrary renormalization prescription. (In fact, both parts of this  equation do not depend on this prescription in accordance with what was mentioned above.) In contrast, RGFs defined by eq. (2.9) depend on both a regularization and a subtraction scheme. For them the equation similar to (2.10) is valid only for some certain renormalization prescriptions called "the NSVZ schemes". Some of them are given by the HD+MSL renormalization prescription [46,64,65], when the theory is regularized by higher derivatives and divergences are removed by minimal subtractions of logarithms [47,48].

The three-loop anomalous dimension and the four-loop β-function defined in terms of the bare coupling constant
Calculating superdiagrams with two external lines of the matter superfields one can obtain the function G related to the corresponding part of the effective action by the equation If the function G is known, then the anomalous dimension defined in terms of the bare coupling constant can be obtained with the help of the equation where the condition p → 0 removes terms proportional to powers of p/Λ. The superdiagrams contributing to the three-loop anomalous dimension of the matter superfields were generated by a special C++ program written by I.S. Also this program calculated the D-algebra and reduced the expression for a superdiagram to a momentum integral using the standard technique for calculating supergraphs, see, e.g., [66][67][68]. The result obtained with the help of this program is written as where R K ≡ R(K 2 /Λ 2 ). Note that the calculation was made for an arbitrary value of the gauge parameter ξ 0 , but the gauge dependence disappeared as should be according to the general theorems, see [69] and references therein. Certainly, this fact can be considered as a useful correctness check.
Calculating the expression (3.3) one should express the bare coupling constant in terms of the renormalized one and only after this make the differentiation with respect to ln Λ, which should certainly be done before integrations. In the lowest orders the bare and renormalized coupling constants (α 0 = e 2 0 /4π and α = e 2 /4π, respectively) are related by the equation (see, e.g., [58]) where b 1,0 and b 2,0 are finite constants which (partially) determine a subtraction scheme in the one-and two-loop approximations, respectively. Note that similar finite constants appear in the renormalization constant for the chiral matter superfields. For instance, in the lowest approximation The integrals present in the expression (3.3) (rewritten in terms of the renormalized coupling constant with the help of eq. (3.4)) are calculated in appendix A using the Chebyshev polynomials method [61]. The Chebyshev polynomials are defined as C n (cos θ) ≡ sin ((n + 1)θ) sin θ (3.6) and (for t < 1) satisfy the important equation Consequently, the function (K − L) −2 = (K 2 − 2KL cos θ + L 2 ) −1 , where θ is the angle between the Euclidian four-vectors K µ and L µ , can be presented in the form

JHEP04(2022)108
Then the angular parts of the loops integrals can be calculated with the help of the useful identities where dΩ is the element of a solid angle on a sphere S 3 in the momentum space. (The subscript Q indicates that this sphere lives in the space with the Cartesian coordinates The calculation of the integrals in eq. (3.3) is described in appendix A. The resulting expression for the three-loop anomalous dimension of the matter superfields defined in terms of the bare coupling constant is written as where A 1 , A 2 , C, D 1 , and D 2 are numerical parameters depending on the regulator function R(x). They are defined by the following equations: Note that the expression (3.11) depends only on the regularization parameters, but is independent of the renormalization parameters which fix a subtraction scheme in the considered approximation. Certainly, this occurs because we consider the anomalous dimension defined in terms of the bare coupling constant, which is scheme independent for a fixed regularization [46].
In the particular case R(x) = 1 + x n some of the regularization parameters can be calculated analytically, where ψ(z) is the logarithmic derivative of the gamma function Γ(z), (3.14) -7 -

JHEP04(2022)108
For n = 1, 2 the expression for the constant C can be simplified, Unfortunately we did manage to obtain an analytic result for the constant D 2 even for However, this constant can easily be calculated numerically, for instance, To find the β-function defined in terms of the bare coupling constant, we substitute the expression (3.11) into the NSVZ equation (2.10). According to [22,39], for RGFs defined in terms of the bare coupling constant it is valid in all loops for an arbitrary renormalization prescription supplementing the higher derivative regularization. Therefore, the four-loop β-function takes the form Certainly, as any renormalization group function defined in terms of the bare couplings, this expression depends only on the regularization parameters, but is independent of the parameters which determine a subtraction scheme (for a fixed regularization).

The algorithm for obtaining RGFs defined in terms of the renormalized coupling constant
Starting from RGFs defined in terms of the bare coupling constant one can obtain RGFs defined in terms of the renormalized coupling constant. For this purpose it is necessary to do the following: 1. Using the expressions for α(α 0 ) and ln Z in the previous orders, the functions β(α 0 ) and γ(α 0 ) should be rewritten in terms of the renormalized coupling constant α. The results are substituted into the left hand side of the equations (2.8).
5. Resulting RGFs should be written in terms of the renormalized coupling constant. After this all ln Λ/µ should cancel each other providing a test of the calculation correctness.

The three-loop anomalous dimension defined in terms of the renormalized coupling constant
For the anomalous dimension of the matter superfields this procedure is as follows: integrating the equation we obtain the expression for ln Z written in terms of the renormalized coupling constant, Note that in this expression the finite constants fixing a subtraction scheme are introduced in such a way that the terms containing them will be polynomials in N f with the same maximal power as the quantum correction in the corresponding approximation. This implies that the renormalization prescription is compatible with the structure of quantum corrections. Rewriting ln Z in terms of the bare coupling constant and differentiating the result with respect to ln µ we construct the anomalous dimension defined in terms of the renormalized coupling constant, For the values of finite constants

JHEP04(2022)108
the expression (4.3) produces the result in the DR scheme, which was found in ref. [5]. Note that the terms which do not contain N f are scheme independent and should be the same in all subtraction schemes [64]. The coincidence of these terms in eqs. (4.3) and (4.5) can be considered as a test of the calculation correction.
Another interesting observation is that the finite constants which determine the renormalization of the matter superfields can be chosen so as to set all scheme dependent terms (which are proportional to (N f ) k with k ≥ 1) to 0. Really, if then all scheme dependent terms vanish and the anomalous dimension of the matter superfields takes the simplest possible form Below this renormalization prescription will be called "the minimal scheme". Actually, by a special tuning of a renormalization prescription we obtained the result coinciding with the one in the so-called "conformal symmetry limit" if we adopt the terminology of ref. [70]. However, instead of omitting certain supergraphs here we set the corresponding contributions to 0 by a proper choice of finite constants fixing a subtraction scheme. Note that one can make such a finite renormalization that all contributions to the anomalous dimension beyond the one-loop approximation vanish. This choice corresponds to the so-called "t'Hooft scheme" [71,72], see also [73] and references therein. For example, the two-and three-loop contributions to the anomalous dimension can be set to 0 if we choose The presence of the terms proportional to (N f ) −1 in these expressions tells us that the finite renormalization which gives the t'Hooft scheme is not compatible with the structure of quantum corrections, because such a dependence on N f cannot appear in any supergraph contributing to the considered renormalization group function. Therefore, it is eq. (4.7) that is the simplest expression for the anomalous dimensions among those that can be obtained by finite renormalizations compatible with the structure of quantum corrections.

The four-loop β-function defined in terms of the renormalized coupling constant
To obtain the four-loop β-function defined in terms of the renormalized coupling constant, we first integrate the equation Into the right hand side of this equation we substitute the β-function (3.18), which was derived from the NSVZ relation. As a result we obtain the dependence α 0 (α) which again contains some finite constants which specify a subtraction scheme in the considered approximation, Solving this equation for 1/α and differentiating the result with respect to ln µ we construct the required β-function, The DR result [4] β DR (α) is obtained from this expression if the finite constants which determine a subtraction scheme take the values Certainly, all scheme independent terms in eq. (4.11) and (4.12) coincide. However, this is rather evident, because such terms coincided in the expression for the three-loop anomalous dimension. Really, the scheme independent terms in the L-loop β-function and in the (L − 1)-loop anomalous dimension are related by the NSVZ equation. In some schemes the NSVZ equation takes place according to the results of [22,46]. Therefore, it should be valid for the terms which do not depend on a renormalization presciption.
Also we see that all scheme dependent terms can be set to 0 if the finite constants are chosen so as to satisfy the equations (4.14) In this minimal scheme the β-function (certainly defined in terms of the renormalized coupling constant) takes the simplest form Evidently, the minimal scheme is NSVZ. In the considered approximation this can be seen from eqs. (4.7) and (4.15). In all loops this is also true, because due to the existence of the all-loop NSVZ schemes (constructed, e.g., in [46,49]) the scheme independent consequences of the NSVZ relation should be satisfied [64].
For completeness, we also present the finite renormalization producing the t'Hooft scheme (in which all terms in the β-function starting from the three-loop approximation are set to 0) in the considered approximation, Certainly, due to the presence of the terms which do not contain N f this finite renormalization is not compatible with the structure of quantum correction, because any supergraph contributing to the β-function of the considered theory is evidently proportional to N f . Therefore, the simplest expression for the β-function in the four-loop approximation is given by eq. (4.15).

The NSVZ relation for RGFs defined in terms of the renormalized coupling constant
The scheme dependent RGFs defined in terms of the renormalized coupling constant satisfy the NSVZ equation

JHEP04(2022)108
only for certain renormalization prescriptions called "the NSVZ schemes". According to [74] (see also [43,75,76] for the generalization to the non-Abelian case) these schemes constitute a continuous set and are related by finite renormalizations α = α (α), Z = z(α)Z which satisfy the equation where B is a finite constant.
Comparing the four-loop β-function (4.11) and the three-loop anomalous dimension (4.3) we see that the NSVZ equation (4.17) is valid if and only if the finite constants b i and g i satisfy the relations It is easy to see that these equations follows from eq. (4.18). Really, in the HD+MSL scheme all finite constants vanish and RGFs defined in terms of the bare coupling constant coincide with the ones defined in terms of the renormalized coupling constant after the formal replacement α 0 → α. (In the considered approximation this can be verified explicitly by comparing the corresponding expressions.) Certainly, the HD+MSL scheme is NSVZ. This follows from the general result of [46] which is based on the statement that RGFs defined in terms of the bare couplings satisfy the NSVZ equation in all loops in the case of using the higher derivative regularization. The scheme determined by eqs. (4.2) and (4.10) can be obtained from the HD+MSL scheme as a result of the finite renormalization with

General statements
In the previous sections we saw that all terms in the anomalous dimension proportional to (N f ) k with k ≥ 1 and all terms in the β-function propotional to (N f ) k with k ≥ 2 can be set to 0 by a special choice of a subtraction scheme. Therefore, there is a certain finite renormalization α → α (α); Z → Z = z(α)Z which removes these terms. It is very important that the form of this finite renormalization is compatible with the structure of quantum corrections. In (S)QED with N f flavors this implies that

JHEP04(2022)108
(The second condition means that z(α) should not contain terms with negative powers of N f .) Here we will demonstrate that the above statement is valid not only in the considered approximations (three loops for the anomalous dimension and four loops for the β-function), but also in all orders. As a starting point we split RGFs into the parts proportional to various powers of N f , The equations which describe how RGFs transform under finite renormalizations were constructed in [77,78] and are written as According to [64] the functions β 1 (α) and γ 0 (α) are scheme independent and remain invariant under these finite renormalizations. Moreover, in the supersymmetric case they satisfy the relation which follows from the NSVZ equation and is its scheme independent consequence. Below will construct a renormalization scheme in which all scheme dependent terms in RGFs (defined in terms of the renormalized coupling constant) are set to 0. For the β-function this implies that β (α ) = N f β 1 (α ). However, it is also very important that the difference 1/α − 1/α should be proportional to N f , because only in this case the considered finite renormalization is compatible with the structure of quantum corrections. To see this, we note that . Therefore, eq. (5.8) (for the considered theory) can equivalently be rewritten in the form where c k are constants which do not depend on N f . The integration constant in (5.9) was denoted by N f b 1,0 /π in order that eq. (4.20) will be valid in the lowest order. From eq. (5.9) it is already clear that the difference 1/α − 1/α is proportional to N f , so that JHEP04(2022)108 the finite renormalization determined by eq. (5.8) is really compatible with the structure of quantum corrections. In the next section we will illustrate this general statement by an explicit calculation. Similarly, if we set to 0 all scheme dependent terms in the anomalous dimension of the matter superfields, then and from eq. (5.5) we obtain the required finite renormalization of the matter superfields, It is evident that both the numerator and the denominator of the integrand in the right hand side are proportional to α 2 N f . Taking into account that in the one-loop approximation β(α) = α 2 N f /π we see that the right hand side of eq. (5.11) can be presented as a series in α in which the coefficients do not contain negative powers of N f , so that Thus, we conclude that using the finite renormalizations compatible with the structure of quantum corrections it is possible to choose such a subtraction scheme in which in all orders

The minimal scheme in the lowest orders
Let us illustrate the general argumentation presented in the previous section by an explicit calculation. Namely, we would like to construct such a finite renormalization that transfer the HD+MSL scheme into the minimal scheme and compare the result with eqs. (4.6) and (4.14). In the HD+MSL scheme all finite constants g i and b i are equal to 0, so that RGFs in this scheme are given by the expressions The scheme in which the anomalous dimension and the β-function are given by eqs. (4.3) and (4.11) is obtained from the HD+MSL scheme after the finite renormalization given by eqs. (4.20) and (4.21). Here our purpose is to find the values of the finite constants g i and b i which correspond to the minimal scheme using the method described in section 5.1. By

JHEP04(2022)108
other words, the new three-loop anomalous dimension and the new four-loop β-function should be given by the expressions where α is determined by eq. (4.20). Substituting the expressions (5.15) and (5.17) into eq. (5.8) and integrating two resulting series we obtain the equation The coupling constant α in the left hand side should be expressed in terms of α with the help of eq. (4.20). After this, equating the coefficients at different powers of alpha we find the values of the coefficients b i , which coincide with the ones defined by eq. (4.14). Similarly, substituting the expressions (5.14), (5.15), and (5.16) into eq. (5.11) after some transformations we construct the expression for ln z(α) in the considered approximation, Comparing this expression with eq. (4.21) we obtain the values of the finite constants corresponding to the minimal scheme for the anomalous dimension of the matter superfields,

JHEP04(2022)108
Again, it is easy to see that this system of equations is equivalent to eq. (4.6) derived directly from this requirement. Therefore, the general argumentation presented in this section really works and can really reproduce the results in the lowest approximations.

Conclusion
In this paper we have calculated the three-loop anomalous dimension of the matter superfields and the four-loop β-function for N = 1 SQED with N f flavors regularized by higher derivatives. As a starting point we construct the integrals which give the threeloop anomalous dimension defined in terms of the bare coupling constant with the help of a special computer program written by I.S. As a correctness test we have verified the cancellation of the gauge dependence. In the supersymmetric case it takes place because (due to the non-renormalization of the superpotential [79]) the anomalous dimension of the chiral matter superfields is related to the anomalous dimension of their mass. The obtained integrals have been calculated using the Chebyshev polynomial method [61]. The result for the scheme independent part of the three-loop contribution (which does not contain N f ) coincided with the one obtained with the use of dimensional reduction [5]. The parts proportional to N f and (N f ) 2 turn out to depend on regularization parameters. Some of them can be calculated analytically for the simplest regulator functions, but one of them (in the (N f ) 2 contribution) can be obtained only numerically.
The four-loop β-function defined in terms of the bare coupling constant can immediately be found using the general statement [22,39] that the NSVZ equation is valid for RGFs defined in terms of the bare coupling constant in all orders in the case of using the higher derivative regularization. (Due to the scheme independence of these RGFs this is true for an arbitrary renormalization prescription supplementing this regularization.) The three-loop anomalous dimension of the matter superfields and the four-loop βfunction defined in terms of the renormalized coupling constant have been obtained from the corresponding RGFs defined in terms of the bare coupling constant in an arbitrary subtraction scheme. For a certain renormalization prescription the results coincided with the ones in DR scheme. Also it turned out that by a special choice of a renormalization prescription compatible with the structure of quantum corrections one can set to 0 all terms proportional to (N f ) k with k ≥ 1 in the anomalous dimension and all terms proportional to (N f ) k with k ≥ 2 in the β-function. We have demonstrated that this can be done in all orders. Note that the proof is also valid for other Abelian gauge theories, say, for usual QED with N f flavors. In this "minimal" scheme RGFs defined in terms of the renormalized couplings take the simplest form. For N = 1 SQED this minimal scheme is NSVZ in all orders, and only scheme independent terms survive in it.
We believe that it would be also interesting to generalize the results of this paper to the non-Abelian case. Possibly, they will also be useful for investigating fixed points in supersymmetric theories, see, e.g., [80][81][82] and references therein.

A.1 Integrals which do not contain N f
First, let us calculate a part of the three-loop contribution to the anomalous dimension of the matter superfields which does not contain N f . It is important that all terms coming from the renormalization of the coupling constant in the previous orders are proportional to N f or (N f ) 2 . Therefore, the considered contribution can be written in the form where we took into account that the difference between e 6 0 and e 6 is of the order e 8 and can be ignored.
Making in the two parts of the last term the changes of integration variables K µ → K µ − L µ ; L µ → −L µ and Q µ → Q µ − L µ we can present this term in the form (A.2) To simplify this expression, we note that because both integrals (without the derivative with respect to ln Λ) are finite in both infrared and ultraviolet regions. This implies that they are finite constants which do not depend on Λ and, therefore, vanish after the differentiation with respect to ln Λ.

JHEP04(2022)108
In the first term the angle integration is trivial, so that it can equivalently be rewritten as It is convenient to break the integration domain into 6 parts defined by the conditions Q > K > L, Q > L > K, K > Q > L, K > L > Q, L > K > Q, L > Q > K.
Evidently, due to the symmetry of the integrand integrations over all these domains give the same results, so that the expression (A.6) can be presented as The second term in the expression (A.5) can be calculated using the technique of the Chebyshev polynomials. Namely, using eqs. (3.8) (for (Q + L) −2 ) and (3.10) it can be presented as In the first term Q > L and, therefore, the integration domain consists of three parts: K > Q > L, Q > K > L, and Q > L > K. Similarly, in the second term L > Q and the integration domain consists of the parts K > L > Q, L > K > Q, and L > Q > K. Breaking the integral into the sum of integrals over these domains, changing the order of integrations, and renaming the integration variables we can be rewrite the expression (A.8) in the form The third term in the expression (A.5) can also be calculated with the help of Chebyshev polynomilas. Using eqs. (3.8), (3.9), and (3.10) it can be presented in the form Summing up the integrals (A.7), (A.9), and (A.10) we obtain that the considered contribution to the three-loop anomalous dimension is given by the expression (A.11)

JHEP04(2022)108
After the change of the integration variable Q = xL it can be rewritten as (A. 12) In this expression we again make the change of the integration variable L = yK. The result takes the form All regulator functions R present in this expression depend on K and Λ only in the combination K/Λ. Therefore, the derivative with respect to ln Λ can be expressed in terms of the derivative with respect to ln K, .14) and the resulting integral over K turns out to be the integral of a total derivative, Taking into account that the regulator function R(x) is equal to 1 at x = 0 and rapidly increases at infinity, we obtain the resulting expression for a part of the three-loop anomalous dimension which does not depend on N f , According to [64], this contribution is scheme independent and, therefore, should coincide with the one obtained in the DR scheme. This coincidence really takes place and confirmes the correctness of the calculation.

A.2 Integrals linear in N f containing an insertion of the one-loop polarization operator
Next, let us consider a part of the contribution proportional to (N f ) 1 which comes from the supergraphs containing an insertion of the one-loop polarization operator. Note that in this case it is necessary to take into account the renormalization of the coupling constant in the one-loop approximation. The corresponding contribution is obtained when in the two-loop contribution the bare coupling constant is expressed in terms of the renormalized one with the help of eq. (3.4). Then the resulting expression for considered part of the three-loop anomalous dimension can be written as

JHEP04(2022)108
It is convenient to introduce the function From the explicit form of the function p(K) it s easy to see that Let us rewrite the expression (A.17) in the form and calculate the integrals entering it. First, using the technique of Chebyshev polynomials we present the integral in the form Making the change of the integration variable L ≡ xK and taking into account that the function (R K R xK ) −1 depends on K and Λ only in the combination K/Λ we can calculate this integral for an arbitrary regulator function R, Next, we consider the integral and calculate its angular part using the Chebyshev polynomials. Then, after some transformations, it can be rewritten in the form (A.25)

JHEP04(2022)108
After the change of the integration variable L ≡ xK this integral can be reduced to the integral of the total derivative with respect to K, Calculating it and taking eq. (A.19) into account we obtain The remaining integral in the expression (A.20), after performing integration by angles with the help of the Chebyshev polynomials takes the form The part of this expression containing ln x can be calculated for an arbitrary regulator function R, However, unfortunately, we did not manage to do this for the remaining contribution. After the change of the integration variable y = K 2 /Λ 2 it can be reduced to the expression proportional to the constant C defined by eq. (3.12), For the regulator function R(x) = 1 + x n the expression for this constant can be expressed in terms of the logarithmic derivative of the gamma-function. For this purpose we first calculate the integral over y,

JHEP04(2022)108
and then make the substitution t = x 2n , after which the considered expression takes the form (A.34)

A.3 Integrals linear in N f containing an insertion of the two-loop polarization operator
The second part of the contribution to the three-loop anomalous dimension proportional to (N f ) 1 comes from superdiagrams containing an insertion of the two-loop polarization operator of the quantum gauge superfield. For calculating it we should also take into account the corresponding term in the renormalization of the coupling constant. Then the expression for the considered contribution can be written as It is convenient to introduce the function Then the expression (A.35) can equivalently be presented in the form To calculate the integrals which enter it, we use the identity

JHEP04(2022)108
valid for a nonsingular function f (K/Λ) which rapidly decreases at infinity. With the help of eq. (A.38) we immediately obtain The remaining integral 40) is reduced to the constant A 1 (defined by eq. (3.12)) after the change of the integration variable x = K 2 /Λ 2 , For the regulator function R(x) = 1 + x n this constant vanishes. Therefore, it remains only to calculate the constant q(0). For this purpose we first present q(K) as a sum q(K) = q 1 (K) + q 2 (K) , (A. 42) where the functions q 1 (K) and q 2 (K) are defined by the equations − 2L µ (Q + L) µ L 2 Q 2 (Q + L) 2 (Q + K + L) 2 ; (A.44) Taking the limit K → 0 in the expression for the function q 1 (K) we obtain that the value q 1 (0) can be written as a vanishing integral of a total derivative in the momentum space, The term with the Pauli-Villars masses vanishes due to the absence of singular contributions, while the first two terms cancel each other. (The first one is reduced to the integral over an infinitely large sphere S 3 ∞ , while the second one is reduced to the surface integral over an infinitely small sphere S 3 ε surrounding the point Q µ = 0.)

JHEP04(2022)108
Thus, to obtain the value q(0), it remains to calculate q 2 (0). After the changes of the integration variables Q µ → Q µ −L µ , L µ → −L µ this expression can be presented in the form The angular integrations in this expression can be made with the help of the Chebyshev polynomials. After some transformations the result can be presented in the form Some integrals in this expression can be calculated. (In some of them it is necessary to take into account that lim K→0 R(K 2 /Λ 2 ) = 1.) Then we obtain (A. 48) In the limit K → 0 the last integral can be taken after the substitution L = xK, Also we note that the constant A 1 can equivalently be rewritten in the form (A.50) where we made the change of the integration variable x = L 2 /Λ 2 and took ε = K 2 /Λ 2 . Using eqs. (A.49) and (A.50) we obtain so that Therefore, the considered contribution to the anomalous dimension of the matter superfields takes the form I 1,2 = − ln Λ µ + b 2,0 I 1,2,1 + I 1,2,2 + I 1,2,3 = − ln Λ µ each other. This fact can be considered as test of the calculation correctness. This is really so, and the final result for the three-loop anomalous dimension defined in terms of the bare coupling constant is given by the expression (3.11), which is also presented here for completeness, As it should be [46], this expression does not depend on the finite constants b 1,0 and b 2,0 which (partially) determine the subtraction scheme in the lowest orders.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.