Two-loop renormalization of the matter superfields and finiteness of ${\cal N}=1$ supersymmetric gauge theories regularized by higher derivatives

The two-loop anomalous dimension of the chiral matter superfields is calculated for a general ${\cal N}=1$ supersymmetric gauge theory regularized by higher covariant derivatives. We obtain both the anomalous dimension defined in terms of the bare couplings, and the one defined in terms of the renormalized couplings for an arbitrary renormalization prescription. For the one-loop finite theories we find a simple relation between the higher derivative regulators under which the anomalous dimension defined in terms of the bare couplings vanishes in the considered approximation. In this case the one-loop finite theory is also two-loop finite in the HD+MSL scheme. Using the assumption that with the higher covariant derivative regularization the NSVZ equation is satisfied for RGFs defined in terms of the bare couplings, we construct the expression for the three-loop $\beta$-function. Again, the result is written both for the $\beta$-function defined in terms of the bare couplings and for the one defined in terms of the renormalized couplings for an arbitrary renormalization prescription.


Introduction
Quantum properties of supersymmetric theories have a lot of very interesting features. The maximally extended N = 4 supersymmetric Yang-Mills theory (SYM) is finite in all loops [1,2,3,4], while N = 2 supersymmetric gauge theories are finite beyond the one-loop approximation [1,4,5]. This implies that N = 2 theories finite in the one-loop approximation are also finite in all orders [6]. The one-loop finiteness can be achieved by a special choice of a gauge group and a representation for the matter superfields. In N = 1 supersymmetric theories the divergent quantum corrections can exist in all orders. However, supersymmetry leads to some interesting relations between various renormalization constants. For example, due to the finiteness of the superpotential [7] the renormalizations of masses and Yukawa couplings are related to the renormalization of the chiral matter superfields. Similarly, the non-renormalization of the triple gauge-ghost vertices [8] allows choosing a renormalization prescription for which where Z α , Z c , and Z V are the renormalization constants for the gauge coupling constant, the Faddeev-Popov ghosts, and the quantum gauge superfield, respectively. However, the most interesting quantum feature of N = 1 supersymmetric gauge theories is the existence of a relation between the β-function and the anomalous dimension of the matter superfields called "the exact NSVZ β-function" [9,10,11,12]. Usually it is written as 1 β(α, λ) = − α 2 (3C 2 − T (R) + C(R) i j (γ φ ) j i (α, λ)/r) where r is the dimension of the gauge group G with the structure constants f ABC , and the group Casimirs are defined as In our notation T A are the generators of the representation to which the matter superfields belong. They should be distinguished from the generators of the fundamental representation denoted by t A , which are assumed to be normalized by the condition tr(t A t B ) = δ AB /2. If N = 2 supersymmetric gauge theories are considered as a particular case of N = 1 theories, then the NSVZ relation leads to the finiteness beyond the one-loop approximation, provided the quantization is manifestly N = 2 supersymmetric [13,14]. A natural way to provide this is to use harmonic superspace [15,16] and an invariant regularization [17]. In this case we will also obtain the all-loop finiteness of N = 4 SYM theory as a consequence of the NSVZ equation.
It is important to recall that the NSVZ equation is valid only for certain renormalization prescriptions which are usually called "the NSVZ schemes". In particular, the most popular DR scheme (which is obtained in the case of using dimensional reduction [18] supplemented by modified minimal subtraction [19]) is not NSVZ [20,21,22,23,24]. The MOM scheme is also not NSVZ [25]. Nevertheless, the DR calculations implicitly confirm the NSVZ equation, because they demonstrate the validity of scheme-independent consequences following from the NSVZ equation [25,26]. These scheme independent relations appear because some terms in the renormalization group functions (RGFs) remain invariant under finite renormalizations, see, e.g., [27]. The fact that these relations are satisfied indicates the existence of NSVZ schemes. In fact, there are an infinite number of the NSVZ schemes which constitute a continuous set. In the Abelian case this set has been described in Ref. [28] and, in particular, includes the on-shell [29] and HD+MSL [30] schemes. The latter scheme is obtained if a theory is regularized by higher covariant derivatives [31,32] (see Refs. [33,34] for supersymmetric versions of this regularization) and the divergences are removed by minimal subtractions of logarithms, when only powers of ln Λ/µ are included into renormalization constants [35,36]. Equivalently, this scheme can be introduced by imposing certain boundary conditions on the renormalization constants [30]. Presumably, the HD+MSL scheme is also NSVZ in the non-Abelian case [8]. This is confirmed by some explicit calculations made in such an approximation where the scheme dependence is essential [37,38,39].
The reason why the HD+MSL scheme turns out to be NSVZ is that the NSVZ equation is satisfied by RGFs defined in terms of the bare couplings for theories regularized by higher covariant derivatives. Such RGFs are independent of a renormalization prescription for a fixed regularization, so that Eq. (2) holds for them in an arbitrary substraction scheme. However, the calculation of Ref. [40] indicates that with dimensional reduction these RGFs do not satisfy the NSVZ equation starting from the three-loop approximation for the β-function. For Abelian theories regularized by higher derivatives the validity of the NSVZ equation for RGFs defined in terms of the bare couplings has been proved in all loops in Refs. [41,42]. This proof is based on the fact that the loop integrals which give the β-function defined in terms of the bare coupling constant in supersymmetric theories are integrals of double total derivatives with respect to the loop momenta. The factorization into total and double total derivatives has first been noted in calculating the lowest-order quantum corrections for N = 1 supersymmetric electrodynamics (SQED) in Refs. [43] and [44], respectively. The subsequent generalizations of the proof made in Ref. [41] allowed deriving the all-loop NSVZ-like equations for the Adler D-function [45] in N = 1 supersymmetric chromodynamics [46,47] and for the renormalization of the photino mass in softly broken N = 1 SQED [48]. (In softly broken supersymmetric theories the NSVZ-like equation for the renormalization of the gaugino mass has first been found in [49,50,51].) In both cases the HD+MSL scheme is NSVZ [52,53], because in this scheme RGFs defined in terms of the renormalized couplings coincide with RGFs defined in terms of the bare couplings up to the renaming of arguments, see Eq. (30) below.
In the non-Abelian case the all-loop proof of the NSVZ equation by a similar method has not yet been completed, although its main steps are at present quite clear. First, the NSVZ equation should be rewritten in the equivalent form [8] β(α, λ) with the help of Eq. (1), where γ c and γ V are the anomalous dimensions of the Faddeev-Popov ghosts and the quantum gauge superfield, respectively. According to [54] the β-function is given by integrals of double total derivatives with respect to the loop momenta in all loops if the higher derivatives are used for regularization. (Again, in the case of using dimensional reduction it is not so [55].) Certainly, this is confirmed by a large number of explicit calculations [37,38,56,57,58]. Due to this structure of the loop integrals the β-function beyond the one-loop approximation is given by a sum of δ-singularities, which appear when double total derivatives act on an inverse squared momentum. The all-loop sums of singularities which occur when the double total derivatives act on inverse squared momenta of the matter superfields and of the Faddeev-Popov ghosts give the corresponding anomalous dimensions in Eq. (4) [59]. For remaining singularities produced by momenta of the quantum gauge superfield the corresponding paper is in preparation. Thus, there are strong indications that the NSVZ equation is satisfied by RGFs defined in terms of the bare couplings for theories regularized by higher covariant derivatives independently of a way of renormalization. The above discussion demonstrates that the higher covariant derivative regularization helps to reveal the underlying structure of quantum corrections which is responsible for the appearance of the NSVZ equation in perturbation theory. That is why it is especially interesting to use it for investigating N = 1 finite theories (see Ref. [60] for a recent review of the theoretical aspects and the phenomenological applications). The direct calculation of Ref. [61] made in the DR scheme demonstrated that if an N = 1 supersymmetric gauge theory is finite in the one-loop approximation, then it is also finite in the two-loop approximation. The same result was obtained from arguments based on anomalies [62,63]. According to [64], for theories finite in the L-th loop the β-function vanishes in the (L+1)-th loop. The same statement immediately follows from the NSVZ equation. We know that the NSVZ equation naturally appears with the higher covariant derivative regularization and is not valid in the DR scheme. Therefore, we are tempted to suggest that the higher derivative regularization could reveal some features of one-loop finite N = 1 supersymmetric theories leading to their possible all-loop finiteness. For this purpose one should investigate the anomalous dimension of the matter superfields. The calculation made in Ref. [65] demonstrated that in the DR scheme the three-loop anomalous dimension does not vanish. However, it has explicitly been verified [66] that it is possible to tune a subtraction scheme in such a way that a one-loop finite theory will also be finite in the three-loop approximation. According to the general argumentation of Refs. [67,68,69,70], a scheme in which a one-loop finite theory is finite in all loops should exist, but at present there is no simple prescription for constructing it. Possibly, the use of the higher covariant derivative regularization could help to solve this problem.
One more interesting subject is the possible existence of the exact expression for the anomalous dimension of the matter superfields for theories obeying the P = 1 3 Q constraint proposed by Jack and Jones in Ref. [71]. In our notation it can be written as According to Ref. [66], for theories which satisfy the condition (5) the anomalous dimension of the matter superfields can possibly be written in the Jack, Jones, North (JJN) form while the β-function does not depend on the Yukawa couplings and is given by the geometric series (Note that here we again do not specify a definition of RGFs, a regularization, and a renormalization prescription.) Although the three-loop calculation of Ref. [66] presumably excludes this possibility, this particular case seems to be very interesting and deserving a further investigation. As a justification, it is possible to suggest that Eq. (6) can be valid in higher loops if some more constraints are imposed on the theory together with Eq. (5).
In this paper we consider a general renormalizable N = 1 supersymmetric theory regularized by higher covariant derivatives, which is described in Sect. 2. In Sect. 3 for this theory we calculate the two-loop anomalous dimension of the matter superfields defined in terms of the bare couplings. It is demonstrated that there is a simple relation between the higher derivative regulators and the Pauli-Villars masses for which it vanishes for the one-loop finite theories. Also there is a regularization for which Eq. (6) is valid for the anomalous dimension defined in terms of the bare couplings under the P = 1 3 Q constraint. Certainly, the same statements are valid for the anomalous dimension defined in terms of the renormalized couplings in the HD+MSL scheme. For a general renormalization prescription the expression for the anomalous dimension defined in terms of the renormalized couplings is found in Sect. 4. Using the statement that the NSVZ equation is presumably valid for RGFs defined in terms of the bare couplings with the higher derivative regularization, the expression for the three-loop β-function is written in Sect. 5. Again this is done for the β-functions defined both in terms of the bare couplings and in terms of the renormalized couplings. The particular cases of the one-loop finite theories and theories satisfying the constraint (5) are investigated. Also in Sect. 5 we demonstrate that for the one-loop finite theories the NSVZ equation in the considered approximation is valid in an arbitrary subtraction scheme.

The theory under consideration
We consider the N = 1 SYM theory with a simple gauge group G interacting with the chiral matter superfields φ i in a representation R, which can in general be reducible. At the classical level this theory in the massless limit is described by the superfield action [72,73,74] where the supersymmetric gauge field strength is given by the chiral superfield W a = D 2 e −2V D a e 2V /8. Note that V = e 0 V A t A in the first term of the action (8), while V = e 0 V A T A in the second one.
To quantize the theory, one should take into account that the quantum gauge superfield is renormalized in a nonlinear way [75,76,77]. Also it is convenient to use the background field method [78,79,80] in the superfield formulation [1,72], because it produces a manifestly gauge invariant effective action. All this can be achieved by making the replacement where V and V are the background and quantum gauge superfields, respectively. Note that in our notation the latter superfield satisfies the constraint V + = e −2V V e 2V . The parameters of the nonlinear renormalization are included into the function The first nonlinear term in this function has been calculated in Refs. [81,82]. Subsequently it was demonstrated that the presence of the nonlinear renormalization is very essential for the renormalization group equations to be satisfied [83].
After the replacement (9), the action regularized by higher covariant derivatives in the massless limit can be written as where the left and right covariant spinor derivatives are given by the equations respectively. The regulator functions F (x) and R(x) are infinite at infinity and approach 1 at x = 0. Note that in Eq. (11) the gauge superfield strength is defined as In this paper we will use the gauge fixing term containing the parameter ξ 0 . Due to the presence of the background covariant derivatives it is invariant under the background gauge symmetry. According to [58] the one-loop renormalization of the gauge parameter is described by the equation where a 1 is a finite constant which originates from the arbitrariness in choosing a renormalization prescription. In this paper we will use the Feynman gauge ξ = 1 in which The quantization procedure (see, e.g., [73]) also requires introducing the Faddeev-Popov and Nielsen-Kallosh ghosts. Their actions (S FP and S NK , respectively) can be found, e.g., in Ref. [87]. Due to the presence of the higher derivative regulators R(x) and F (x) in the actions (11) and (14) all divergences disappear beyond the one-loop approximation. This is a general feature of the higher covariant derivative regularization, see, e.g., [84]. For regularizing the remaining one-loop divergences one has to introduce the Pauli-Villars determinants [85]. Following Refs. [58,86], we define the generating functional as where the Pauli-Villars determinants are given by the functional integrals Here ϕ a is a set of three commuting chiral superfields in the adjoint representation with the action and Φ i is a multiplet of commuting chiral superfields in a representation R P V that admits a gauge invariant mass term with the action We assume that the invariant tensor M ij satisfies the condition and the Pauli-Villars masses are proportional to the constant Λ in the higher derivative terms, where a ϕ and a are independent of couplings. Then the one-loop divergences and subdivergences will be regularized if the constant c in the generating functional (18) is equal to T (R)/T (R P V ). In this paper we will calculate the anomalous dimension of the matter superfields. To construct it, first, one should consider the part of the effective action which corresponds to self-energy diagrams for the chiral matter superfields. It can written in the form where (G φ ) i j = δ j i + (∆G φ ) i j , and (∆G φ ) i j is of order α 0 or λ 2 0 . Let α ≡ e 2 /4π and λ ijk denote the renormalized gauge and Yukawa coupling constants, respectively. The bare and renormalized chiral matter superfields are related by the equation where the renormalization constant (Z φ ) i j is determined by requiring the finiteness of the product (Z φ ) i k (G φ ) k j expressed in terms of the renormalized couplings in the limit Λ → ∞. Note that this requirement does not allow fixing this constant uniquely. The remaining arbitrariness is removed by choosing a renormalization prescription.
Divergences of a theory are encoded in RGFs. According to [30], it is necessary to distinguish between RGFs defined in terms of the bare couplings and the ones (standardly) defined in terms of the renormalized couplings. In terms of the bare couplings, the definitions of the β-function and the anomalous dimension of the chiral matter superfields read as For a fixed regularization these RGFs are independent of a renormalization procedure. That is why they (presumably) satisfy the NSVZ relation in an arbitrary subtraction scheme supplementing the higher covariant derivative regularization. In terms of the renormalized couplings RGFs are defined by the equations where the derivatives are taken with respect to the logarithm of the renormalization point µ at fixed values of the bare couplings. The β-function (28) and the anomalous dimension (29) are scheme-dependent starting from the three-and two-loop approximation, respectively. Up to the renaming of arguments both definitions of RGFs coincide in the HD+MSL scheme, in which the renormalization constants contain only powers of ln Λ/µ and all finite constants fixing a subtraction scheme are set to 0, (30) RGFs defined in terms of the bare couplings are related to the corresponding Green functions. For example, using the finiteness of the expression (Z φ G φ ) i j the anomalous dimension (27), which will be calculated in this paper, can be presented in the equivalent form where the limit q → 0 is necessary here to get rid of the terms vanishing in the limit Λ → ∞.
To use Eq. (31), one needs to know how the renormalized and bare couplings are related at lower orders. In particular, for calculating the two-loop anomalous dimension of the matter superfields, one needs the one-loop relation between the bare and renormalized couplings. Moreover, we will also need the two-loop relation between the bare and renormalized gauge coupling constants for investigating the three-loop β-function defined in terms of the renormalized couplings in Sect. 5.2. It can be found by integrating Eq. (26) in which we substitute the two-loop β-function. With the regularization considered in this paper it has been calculated in Ref. [87]. The result contains some arbitrary integration constants b i , which reflect the arbitrariness in choosing a subtraction scheme, Due to the non-renormalization theorem of Ref. [7], the superpotential does not receive divergent radiative corrections. This implies that the renormalization of the Yukawa couplings is related to the renormalization of the matter superfields by the equation In this paper we will consider only such renormalization schemes for which it is valid (although, in general, it is possible to construct subtraction schemes breaking this equation). The one-loop expression for the renormalization constant (Z φ ) i j can be written as where g 11 and g 12 are arbitrary constants similar to b i in Eq. (32). Eqs. (33) and (34) determine the one-loop renormalization of the Yukawa couplings, By definition, in the HD+MSL renormalization scheme all constants b i and g i are set to 0.

Two-loop anomalous dimension
Superdiagrams contributing to the two-loop anomalous dimension of the matter superfields can be divided into two parts. The first part contains the superdiagrams without Yukawa vertices. They are presented in Fig. 1. The gray circles in the superdiagrams (17) and (18) denote insertions of the one-loop polarization operator of the quantum gauge superfield, which is equal to the sum of the supergraphs presented in Fig. 2. They have been calculated in Ref. [86]. The superdiagrams of the second part presented in Fig. 3 contain the Yukawa vertices. They have already been calculated in Refs. [37,38] for the higher derivative regulators R(x) = 1 + x m and F (x) = 1 + x n . In this paper the integrals written in Refs. [37,38] are calculated for an arbitrary form of the functions R(x) and F (x). The details of this calculation are presented in Appendix A.
The superdiagrams presented in Fig. 1 have been calculated in Ref. [52] for N = 1 SQCD. It is essential that the action of N = 1 SQCD does not contain the Yukawa term cubic in the matter superfields, so that there is no need for the regulator F (x). That is why the calculation of (1) Figure 1: Superdiagrams without Yukawa vertices contributing to the two-loop anomalous dimension of the matter superfields. In the diagrams (17) and (18) the gray circles denote insertions of the one-loop polarization operator of the quantum gauge superfield.
Ref. [52] has been done for F (x) = 1. However, for theories containing the Yukawa interaction the higher covariant derivative regularization should include F (x) = 1. In this paper we consider the general case. The presence of the function F generates new vertices, which have to be taken into account. For example, the triple gauge-matter vertex is written as where the coefficients f α are determined by the equation Remarkably, the sum of the supergraphs (1) - (16) in Fig. 1 turns out to be independent of F . This regulator is present only in the superdiagrams (17) and (18) inside the expression for the polarization operator. Explicitly writing the sum of the superdiagrams depicted in Fig. 1, the two-loop anomalous dimension of the matter superfields can be presented as (1) ( +λ-dependent terms + higher orders. The explicit expression for the λ-dependent terms can be found in Ref. [38], see also Appendix A. The functions f (k/Λ) and h(k/Λ) are related to the one-loop polarization operator of the quantum gauge superfield in the Feynman gauge, In our notation the polarization operator Π(α 0 , λ 0 , k 2 /Λ 2 ) is defined by the equation where Γ The explicit expressions for the functions f (k/Λ) and h(k/Λ) are rather large. They can be found in Ref. [86]. However, in this paper we need only their asymptotic behavior at small k, where o(1) denotes terms that vanish in the limit k → 0. In Eq. (37) the differentiation with respect to ln Λ is to be performed at fixed values of the renormalized couplings α and λ before the momentum integration. This makes the integrals well-defined and finite in the infrared region.
The integrals giving the two-loop anomalous dimension are calculated in Appendix A. We calculate both the integrals explicitly written in Eq. (37) and the integrals present in the λdependent terms for an arbitrary form of the functions R(x) and F (x). The result for the two-loop anomalous dimension defined in terms of the bare couplings obtained in Appendix A is written as where the constants A and B are given by the integrals For the regulators R(x) = 1 + x m and F (x) = 1 + x n (which were used in Refs. [38,88]) these integrals can be taken, For the particular case of N = 1 SQED with N f flavors Eq. (43) gives This expression agrees with the results of Refs. [26,88,89] obtained for R(x) = 1 + x m (and, therefore, A = 0) and generalizes them to the case of an arbitrary regulator function. Note that according to Ref. [89] Eq. (46) is valid for an arbitrary ξ-gauge. The anomalous dimension (43) is considerably simplified for theories finite in the one-loop approximation, which are obtained if The first equation (which is an analog of the Banks-Zaks prescription [90]) leads to the vanishing one-loop β-function, while the second one follows from the vanishing of the one-loop anomalous dimension of the matter superfields. In this case we obtain This implies that the one-loop finite theories remain finite in the two-loop approximation if the higher derivative regularization is chosen in such a way that A = B and a = a ϕ . The first condition is automatically satisfied if R(x) = F 2 (x). That is why a version of the higher covariant derivative regularization leading to the two-loop finiteness is obtained by imposing the conditions Possibly, this could help to reveal if the finiteness of the considered class of supersymmetric theories takes place in all loops and what is its underlying reason. Another interesting particular case corresponds to the theories which satisfy the P = 1 3 Q condition. In terms of the bare couplings it is written as In this case the anomalous dimension (43) takes the form From this result we see that in the considered approximation the anomalous dimension satisfies the exact JJN equation (6) if the parameters of the higher covariant derivative regularization satisfy the conditions In particular, the first equality is valid for R(x) = F 2 (x), while the second one can be obtained by a special choice of the Pauli-Villars masses. In this case

Two-loop anomalous dimension defined in terms of the renormalized couplings
To obtain the anomalous dimension (standardly) defined in terms of the renormalized couplings, it is necessary to rewrite the left hand side of the equation (27) in terms of the renormalized couplings with the help of Eqs. (32) and (35) and integrate the result with respect to ln Λ. Then we obtain the function ln Z i j (α, λ, ln Λ/µ), which should be rewritten in terms of the bare couplings using Eqs. (32) and (35). According to Eq. (29), the result should be differentiated with respect to ln µ. The final expression for the anomalous dimension is obtained after expressing the bare couplings in terms of the renormalized ones again using Eqs. (32) and (35). As a correctness check, one can verify cancellation of all ln Λ/µ, which follows from the general theory of the renormalization group (see, e.g., [91]). The result of the above described calculation is written as It depends on the finite constants b i and g i , which determine the renormalization prescription. Certainly, this agrees with the statement that the anomalous dimension defined in terms of the renormalized couplings is scheme dependent starting from the two-loop order. However, the two-loop terms proportional to C(R) 2 i j and to λ * iac λ jab λ * bde λ cde are scheme independent in agreement with Refs. [25] and [37], respectively. Now, let us consider particular cases of the general expression (54). First, we write the expression for the anomalous dimension in the HD+MSL scheme, in which all finite constants (namely, g 11 , g 12 , b 11 , and b 12 in the considered approximation) should be set to 0, so that In agreement with the general statement [30], this expression coincides with Eq. (43) up to the renaming of arguments α → α 0 , λ → λ 0 .
Eq. (54) should also agree with the result in the DR-scheme, which can be found, e.g., in Ref. [20]. This means that there must be such finite constants that the anomalous dimension takes the form This is really true for the finite constants satisfying the equations b 11 − g 11 = ln a ϕ + 1 2 (1 + A); b 12 − g 11 = ln a + 1 2 (1 + A); The existence of these values can be considered as a non-trivial check of the calculation correctness. Note that the constant g 11 remains unfixed due to the arbitrariness of choosing the renormalization point µ. Its value can be found by comparing the renormalized one-loop twopoint Green functions for the matter superfields calculated with the higher covariant derivative regularization and in the DR scheme, see Appendix A and Ref. [92] for details. This gives For N = 1 SQED with N f flavors from Eq. (54) we obtain the expression which agrees with Ref. [26]. The N = 1 supersymmetric theories finite in the one-loop approximation are of a special interest. Written in terms of the renormalized couplings the finiteness conditions are In this case the expression (54) takes the form We see that in general the one-loop finiteness does not lead to the two-loop finiteness. At first sight, this result contradicts the anomaly based consideration of Refs. [62,63]. However, it seems that in Refs. [62,63] the use of the DR scheme is implicitly assumed [93]. For this scheme the finite constants are given by Eq. (57), and the anomalous dimension (61) really vanishes. Also it is possible to find other subtraction schemes in which it is so. In particular, a more interesting example is the HD+MSL scheme obtained for the higher covariant derivative regularization (49). The reason for this is that the HD+MSL scheme seems to be NSVZ in all loops. Moreover, the restrictions on the choice of the regularization can presumably reveal a deep structure of a theory needed for the finiteness, say, a certain symmetry underlying it. One more interesting special case is the theories which satisfy the P = 1 3 Q condition (5), under which the anomalous dimension is given by the expression Again we see that there are two ways to obtain the JJN equation (6) in the considered approximation. Namely, it is possible to choose finite constants for which Eq. (57) is valid. The corresponding subtraction scheme is equivalent to the DR scheme. Another way is to use the HD+MSL scheme with the regularization satisfying Eq. (52).

The β-function defined in terms of the bare couplings
As we have already mentioned, there are strong indications that the NSVZ equation is valid in all loops for RGFs defined in terms of the bare couplings in the case of using the higher covariant derivative regularization. If we believe that this is really so, then it is possible to construct the expression for the three-loop β-function, where γ φ,1-loop i j and γ φ,2-loop i j are the one-and two-loop parts of the anomalous dimension, respectively. Substituting these contributions we obtain Certainly, this expression is independent of a renormalization prescription (for a fixed regularization) as any RGF defined in terms of the bare couplings. In the case of N = 1 SQED with N f flavors it gives We see that the three-loop β-function vanishes for the higher covariant derivative regularization with A = B and a ϕ = a, and, in particular, for the regularization (49). For theories satisfying the P = 1 3 Q constraint (50) the three-loop β-function defined in terms of the bare couplings takes the form From this equation we see that in the considered approximation the JJN equation (7) is satisfied by the β-function defined in terms of the bare couplings for the higher derivative regularization with the parameters obeying the conditions (52).

The β-function defined in terms of the renormalized couplings
To find the β-function defined in terms of the renormalized couplings, first, it is necessary to rewrite the right hand side of Eq. (64) in terms of the renormalized couplings and, then, integrate it with respect to ln Λ. The relations between the bare and renormalized gauge and Yukawa couplings are given by Eqs. (32) and (35), respectively. After integrating Eq. (64) with respect to ln Λ at fixed values of the renormalized couplings we find the three-loop expression for the renormalized gauge coupling constant as a function of the bare couplings. Next, the result is differentiated with respect to ln µ (at fixed values of the bare couplings). Finally, the bare couplings are expressed in terms of the renormalized ones. Note that, according to the general theory of the renormalization group, all ln Λ/µ should disappear. This fact can be used as a correctness check.
The procedure described above gives the three-loop β-function defined in terms of the renormalized couplings for a general renormalization prescription supplementing the higher covariant derivative regularization, We see that the terms corresponding to the three-loop contribution depend on the constants b i and g i . Certainly, this follows from the fact that the β-function defined in terms of the renormalized couplings is scheme dependent starting from the three-loop approximation. The terms containing the Yukawa couplings exactly coincide with the ones obtained by the direct calculations of Refs. [37,38]. By definition, in the HD+MSL scheme all finite constants b i and g i are set to 0, so that the expression for the β-function takes the form and up to the renaming of arguments α → α 0 , λ → λ 0 coincides with the β-function defined in terms of the bare couplings given by Eq. (64).
For certain values of the finite constants b i and g i the expression (68) should also reproduce the DR-result first found in Refs. [20,21], which in our notation takes the form In particular, this implies that all scheme independent terms in Eq. (68) should coincide with the corresponding terms in Eq. (70). We see that it is really so. The other terms coincide if From these equations and Eq. (57) one can express all finite constants in terms of g 11 given by Eq. (58), b 11 = g 11 + ln a ϕ + 1 2 (1 + A) = ln a ϕ ; b 12 = g 11 + ln a + 1 2 (1 + A) = ln a; b 21 = g 11 + ln a ϕ + 3 4 The values (72) of the finite constants b i and g i correspond to the DR renormalization scheme.
Comparing Eqs. (68) and (54) we see that for a general renormalization prescription the NSVZ relation does not take place. However, if the finite constants b i and g i satisfy the equations then the NSVZ β-function in the considered approximation is valid for RGFs defined in terms of the renormalized couplings. Evidently, Eq. (73) does not unambiguously fix the values of all finite constants. This implies that there is a class of NSVZ schemes similar to the one in the Abelian case which was described in Ref. [28]. From Eq. (68) we see that the β-function of N = 1 SQED with N f flavors is given by the expression which agrees with Ref. [26].
For the one-loop finite theories (which satisfy the conditions (60)) the three-loop β-function defined in terms of the renormalized couplings takes the form Comparing this expression with Eq. (61) we see that in this case the NSVZ equation is valid in the considered approximation for a general renormalization prescription in agreement with the general statement that for a theory finite in a certain loop the β-function always vanishes in the next loop [64]. Note that in Ref. [64] the DR scheme was essentially used. However, in general, the NSVZ equation is not valid in the DR scheme (see, e.g., [20]). Now we see how this seeming contradiction is solved in the considered approximation. Also we see that for theories satisfying Eq. (60) the three-loop β-function vanishes if ln a ϕ a In particular, these equations are valid in the DR scheme and in the HD+MSL scheme supplementing the regularization with R(x) = F 2 (x) and a ϕ = a. For theories which satisfy the P = 1 3 Q condition (5) the β-function defined in terms of the renormalized couplings is Therefore, Eq. (7) in the considered approximation is valid, for example, in the HD+MSL scheme with the regularization satisfying the conditions (52). Another possibility is the NSVZ scheme constructed from the DR scheme by a special redefinition of the gauge coupling constant in Refs. [20,21,22]. This scheme corresponds to the finite constants

Conclusion
In this paper the two-loop anomalous dimension of the matter superfields has been calculated for a general N = 1 supersymmetric gauge theory regularized by higher covariant derivatives. The result has been found both for the anomalous dimension defined in terms of the bare couplings and for the one defined in terms of the renormalized couplings for an arbitrary renormalization prescription. For theories finite in the one-loop approximation the two-loop anomalous dimension defined in terms of the bare couplings does not in general vanish. However, for a version of the higher derivative regularization with the parameters R(x) = F 2 (x) and a = a ϕ this is so. Possibly, this could help to understand deeper reasons responsible for the finiteness of a theory. Moreover, we analysed a possibility of satisfying the exact equation (6) for the anomalous dimension which was proposed by Jack, Jones, and North for theories obeying the P = 1 3 Q condition. It appears that in the considered approximation for this equation to be valid the regularization parameters should satisfy the constraints (52).
Using the statement that for RGFs defined in terms of the bare couplings the NSVZ equation is valid in the case of using the higher covariant derivative regularization we also construct the expression for the three-loop β-function. Again, we present the results for the β-function defined in terms of the bare couplings and for the one defined in terms of the renormalized couplings. For the one-loop finite theories the three-loop β-function defined in terms of the bare couplings vanishes under the same conditions (49) as the two-loop anomalous dimension. For theories satisfying Eq. (5) the exact expression for the β-function in the form of the geometric series (7) is valid with the regularization (52).
For RGFs defined in terms of the renormalized couplings the above results are valid in the HD+MSL scheme. Also we have demonstrated that for the one-loop finite theories the NSVZ equation in the considered approximation is valid for an arbitrary renormalization prescription.
RGFs obtained in this paper with the help of the higher covariant derivative regularization are in agreement with the ones in the DR scheme in a sense that there are certain values of finite constants fixing the renormalization prescription for which the results of this paper reproduce the DR results.
The (rather large) explicit expressions for the functions f (k/Λ) and h(k/Λ) inside the integrals I 11 and I 12 can be found in Ref. [86].
The integrals I 6 , I 7 , and I 10 have been calculated in Refs. [37], [38], and [29], respectively. The results are written as Calculating the integrals I 8 and I 9 in the four-dimensional spherical coordinates one can relate them to the constants A and B defined by Eq. (44), Earlier the integrals I 3 , I 4 , I 5 and the integrals I 11 , I 12 have been found only for the higher derivative regulators R(x) = 1 + x m and F (x) = 1 + x n , see Refs. [38] and [83], respectively. However, it is possible to generalize the results to the case of arbitrary functions R(x) and F (x). For this purpose we will use the equation where the last equalities in Eqs. (98) and (99) follow from Eqs. (41) and (42), respectively. Now, let us calculate the integral I 3 for a general higher derivative regulator R(x). Introducing the Feynman parameter x and making the substitution k µ → k µ − xl µ the integral in Eq. (96) can be rewritten as where we took into account that Next, it is convenient to change the integration variable to z = k 2 /Λ 2 and introduce the notation ǫ ≡ x(1 − x) l 2 /Λ 2 . Taking into account that ǫ is small, we obtain where the constant A is defined by Eq. (44). Integrating this expression with respect to x and omitting terms vanishing in the limit Λ → ∞, the considered integral can be written as The result for the integral I 3 is calculated by substituting this expression into Eq. (96), According to Ref. [38], because the last two integrals are well defined in the limit k → 0. Therefore, This expression is completely analogous to the integral I 3 . Repeating the calculation described above, we obtain where the constant B is defined by Eq. (44). Thus, the integrals (81) -(92) are given by the following expressions: Substituting them into Eq. (80) the anomalous dimension can be written as Rewriting the result in terms of the bare couplings using Eqs. (32) and (35) we obtain the anomalous dimension (43). Note that all ln Λ/µ and all finite constants cancel each other in the resulting expression. Certainly, this can be considered as a correctness check.