The NSVZ relations for ${\cal N}=1$ supersymmetric theories with multiple gauge couplings

We investigate the NSVZ relations for ${\cal N}=1$ supersymmetric gauge theories with multiple gauge couplings. As examples, we consider MSSM and the flipped $SU(5)$ model, for which they easily reproduce the results for the two-loop $\beta$-functions. For ${\cal N}=1$ SQCD interacting with the Abelian gauge superfield we demonstrate that the NSVZ-like equation for the Adler $D$-function follows from the NSVZ relations. Also we derive all-loop equations describing how the NSVZ equations for theories with multiple gauge couplings change under finite renormalizations. They allow describing a continuous set of NSVZ schemes in which the exact NSVZ $\beta$-functions are valid for all gauge coupling constants. Very likely, this class includes the HD+MSL scheme, which is obtained if a theory is regularized by Higher covariant Derivatives and divergences are removed by Minimal Subtractions of Logarithms. That is why we also discuss how one can construct the higher derivative regularization for theories with multiple gauge couplings. Presumably, this regularization allows to derive the NSVZ equations for such theories in all loops. In this paper we make the first step of this derivation, namely, the NSVZ equations for theories with multiple gauge couplings are rewritten in a new form which relates the $\beta$-functions to the anomalous dimensions of the quantum gauge superfields, of the Faddeev--Popov ghosts, and of the matter superfields. The equivalence of this new form to the original NSVZ relations follows from the extension of the non-renormalization theorem for the triple gauge-ghost vertices, which is also derived in this paper.


Introduction
Although supersymmetry has not yet been discovered in direct experiments, there are some indications that physics beyond the Standard model is supersymmetric [1]. These indications are obtained from the analysis of various quantum corrections, which are compared with experimental data and theoretical considerations. Certainly, the most convincing evidence in favor of supersymmetry is the unification of gauge coupling constants in supersymmetric extensions of the Standard model [2][3][4]. The evolution of couplings is encoded in the β-functions. To see that the gauge coupling unification becomes much better in the supersymmetric case, it is sufficient to consider the lowest (one-loop) contributions to the β-functions. However, the higher order corrections are also very important. For N = 1 supersymmetric gauge theories there is an exact expression for the β-function, which relates it to the anomalous dimension of the matter superfields [5][6][7][8] and is usually called "the exact NSVZ β-function". As a rule, this equation was analysed for theories with a simple gauge group. This implies that in this case there is the only gauge coupling constant α. If Yukawa couplings are denoted by λ, then the NSVZ equation for a theory with a single gauge coupling constant is written in the form where r = dim G is a dimension of the gauge group, Here T A denotes the generators of the representation in which the chiral matter superfields lie. They should be distinguished from the generators of the fundamental representation denoted by t A , which are assumed to be normalized by the condition From the NSVZ equation we see that the all-loop β-function for the pure N = 1 supersymmetric Yang-Mills (SYM) theory is given by the geometric series [6]. Moreover, it is possible to prove [9,10] that the finiteness of N = 4 SYM [11][12][13][14][15] and the finiteness of N = 2 supersymmetric gauge theories beyond the one-loop approximation [14][15][16] also follow from Eq. (1). 1 Note that due to the N = 2 non-renormalization theorem it is possible to construct N = 2 supersymmetric theories finite in all loops [21] choosing a gauge group and a representation for the matter superfields in such a way that the one-loop β-function vanishes. 2 The NSVZ relation can also be applied for investigating N = 1 finite theories [25][26][27][28][29] (see [30] for a recent review). In particular, it allows to prove vanishing of the (L + 1)-loop β-function for theories finite in L-loops [31,32] in a simple way [33]. The exact NSVZ β-function is also useful for studying fixed points in N = 1 supersymmetric theories, see, e.g., [34,35]. Moreover, in softly broken supersymmetric theories the renormalization of the gaugino masses is described by the NSVZ-like relations [36][37][38]. It is important that the NSVZ equation is valid only for certain (NSVZ) renormalization prescriptions, which constitute a continuous set [39][40][41]. In particular, it does not hold in the DR-scheme [42][43][44][45][46], when the theory is regularized by dimensional reduction [47] and divergences are removed by modified minimal subtraction [48]. Also it is not valid in the MOM-scheme [49,50]. Certainly, it is possible to construct an NSVZ scheme by making specially tuned finite renormalizations in each order of the perturbation theory [42][43][44][45][46], see also [51,52]. However, there is a simple renormalization prescription giving some NSVZ schemes in all orders. To obtain it, one should regularize a theory by higher covariant derivatives [53,54] in a supersymmetric way [55,56]. It is important that this regularization also includes the insertion of the Pauli-Villars determinants for removing one-loop divergences which survive after adding a term with higher derivatives to the action [57]. In the supersymmetric case such a construction can be found in Refs. [58,59]. With the Higher covariant Derivative regularization an NSVZ scheme can be obtained if divergences are removed by Minimal Subtractions of Logarithms or, in other words, in the HD+MSL scheme. In this scheme only powers of ln Λ/µ are included into renormalization constants, where Λ is the dimensionful cut-off parameter of the regularized theory and µ is the renormalization point [60,61]. Equivalently, the HD+MSL scheme can be introduced by imposing certain boundary conditions on renormalization constants [62]. Note that minimal subtractions of logarithms can supplement various versions of the higher derivative regularization (which differ in the form of the higher derivative terms and in the Pauli-Villars masses), so that the HD+MSL prescription in general produces a certain set of NSVZ schemes.
The all-loop proof that the HD+MSL prescription gives some NSVZ schemes [63] is based on the all-loop derivation of the NSVZ equation in Refs. [63][64][65], see also Ref. [66]. 3 The NSVZ equation is obtained as follows: 1. According to [64] the triple gauge-ghost vertices (with two ghost external lines and one external line of the quantum gauge superfield) are finite in all orders. (This statement has been verified in the one-and two-loop approximations by explicit calculations made in [64,72]. ) 2. Using the finiteness of the triple gauge-ghost vertices the NSVZ β-function can rewritten in an equivalent form [64], which relates the β-function to the anomalous dimensions of the quantum gauge superfield, of the Faddeev-Popov ghosts, and of the matter superfields, 3. In [65] it was proved that with the higher covariant derivative regularization the integrals giving the β-function are integrals of double total derivatives in the momentum space. 4 4. The integrals of double total derivatives are nontrivial due to singularities of the integrands, which appear due to the identity The sum of these singularities has been calculated in [63]. (The sums of ghost and matter singularities have also been calculated by a different method in Ref. [66].) As the result, the NSVZ equation (4) is obtained for the renormalization group functions (RGFs) defined in terms of the bare couplings for an arbitrary renormalization prescription supplementing the higher covariant derivative regularization. 5. For RGFs defined in terms of the renormalized couplings some NSVZ schemes are obtained using the fact that in the HD+MSL scheme RGFs defined in terms of the bare couplings coincide with the ones defined in terms of the renormalized couplings up to a formal change of the argument [62].
Although the NSVZ equation was usually studied for theories with a single gauge coupling constant, its analogs can be written for theories with multiple gauge couplings (in which gauge groups are given by direct products of simple and/or U (1) factors) [83]. For example, the allloop equations describing the renormalization of the Minimal Supersymmetric Standard Model (MSSM) couplings written in [84] are actually the integrated form of the NSVZ relations in this case. It was verified [85] that they really agree with the explicit two-loop calculations. The NSVZ and NSVZ-like relations were applied in MSSM and its extensions for investigating various quantum corrections [86,87], although for obtaining the higher-loop results in the DR-scheme a proper tuning of a subtraction scheme was certainly needed. Also the NSVZ equations were used for studying the reduction of couplings in theories with multiple gauge couplings, see [30,88,89] and references therein. In Refs. [90,91] the equation similar to the Abelian NSVZ relation [92,93] was written for the Adler D-function [94] in N = 1 supersymmetric quantum chromodynamics (SQCD) interacting with the Abelian gauge superfield. (The theory considered in [90,91] has the gauge group SU (N ) × U (1) and, therefore, two gauge coupling constants.) Thus, the NSVZ equations are important for phenomenologically interesting theories with multiple gauge couplings and deserve a more detailed consideration in this case. In particular, it is necessary to derive them using the methods of the perturbation theory and specify renormalization prescriptions under which they are valid. In this paper we start this investigation considering as examples MSSM (see, e.g., [1]) and the flipped SU (5) Grand Unification Theory (GUT) [95][96][97][98].
The paper is organized as follows. In Sect. 2 we recall a form of the NSVZ relations for theories with multiple gauge couplings. In the next Sect. 3 these equations are written for N = 1 SQCD interacting with the N = 1 supersymmetric Abelian gauge superfield, MSSM, and the flipped SU (5) model. The scheme dependence of the NSVZ relations for theories with multiple gauge couplings is investigated in Sect. 4. We present equations describing how they change under finite renormalizations of couplings and matter superfields, and construct finite renormalizations which transfer an NSVZ scheme into another NSVZ scheme. In particular, this implies that the NSVZ schemes constitute a continuous sets. In theories with a single gauge coupling some schemes of this class in all loops are obtained with the help of the HD+MSL prescription. That is why in Sect. 5 we explain how one can regularize theories with multiple gauge couplings by higher covariant derivatives. The general construction is illustrated by the examples of MSSM and the flipped SU (5) model. This regularization can be used a starting point for the all-loop derivation of the NSVZ equations for the theories under consideration. In theories with a single gauge coupling such a derivation produces a new form of the NSVZ equation (4). Its generalization to the case of multiple gauge couplings is constructed in Sect. 6 using the extension of the non-renormalization theorem for the triple gauge-ghost vertices, which is also derived in this section.

The NSVZ equations for theories with multiple gauge couplings
In this paper we will consider a general renormalizable N = 1 supersymmetric theory with the gauge group where for each K = 1, . . . , n the subgroup G K is either a simple compact group or U (1). The chiral matter superfields can be split into sets (numerated by the index a) in such a way that each of these sets transforms under certain irreducible representations R aK of the simple subgroups G K or has certain charges q aK with respect to G K = U (1). This implies that for a fixed a the matter superfields are also numerated by indices corresponding to various subgroups G K . It is convenient to introduce the index i numerating all matter superfields φ i as the set It is important that for different values of a the sets of i 1 , . . . i n values are in general different. In this notation the generators of the gauge group G corresponding to the subgroup G K can be written as where T A K a i K j K are either the generators of G K in the representation R aK for simple subgroups G K or the charges of the superfields φ i ≡ φ a; i 1 i 2 ...in with respect to G K = U (1). Note that we always assume that the generators are normalized by the conditions where t A K are generators of the fundamental representation of a simple subgroup G K . N = 1 supersymmetric theories with the gauge group (6) in the massless limit are described by the action which is written in terms of N = 1 superfields. Such a formulation is especially convenient, because in this case N = 1 supersymmetry is a manifest symmetry of the theory. Using a proper regularization and quantization procedure [99][100][101], N = 1 supersymmetry can also be made manifest at all steps of calculating quantum corrections. In Eq. (11) V denotes the gauge superfield. In the matter part of the action it is given by the expression where the bare gauge coupling constants are denoted by e 0K . Below we will also use the notation The gauge superfield strengths corresponding to various subgroups G K in the direct product (6) are defined as All expressions (W a ) A K (W a ) A K are invariant under gauge transformations of the whole group G, so that there are n invariants of this structure. Therefore, the theory (11) has n gauge coupling constants α K = e 2 K /4π, each of them corresponding to a certain subgroup of the group G. In terms of the renormalized couplings the corresponding n β-functions are defined by the equations where the subscript 0 marks the bare gauge and Yukawa couplings, and µ stands for the renormalization point.
Taking into account that chiral matter superfields φ a belong to irreducible representations of all simple subgroups G K , the renormalization constants for them have no indices, where the subscript R denotes the renormalized superfields. In terms of the renormalized coupling constant the corresponding anomalous dimensions are defined by the equation Earlier we denoted the generators of the representation R aK by T A K a . This implies that It is important that the right hand side of this equation is proportional to the δ-symbol because, by construction, the representation R aK is irreducible.
→ q aK are numbers.) Also it is useful to define the analog of T (R) with the help of the equation Then the matter contribution to the one-loop β-function corresponding to the coupling constant α K will be proportional to where we introduced the notation Really, calculating the function β K in the one-loop approximation it is necessary to sum over all indices of chiral superfields corresponding to all subgroups except for G K , while the trace of two generators of the subgroup G K gives the factor (20). It is reasonable to suggest (see, e.g., [83]) that the exact β-function for a coupling constant corresponding to a factor G K in the gauge group G is given by the expression (Note that this equation is also valid even if the subgroup G K coincides with U (1). Certainly, in this case it is necessary to set C 2 (G K ) = 0.) Really, if G is a simple group, then the representation to which the matter superfields belong can be presented as a direct sum of R a , so that Therefore, the NSVZ relation (1) in this case can be written in the form which is similar to Eq. (23). Thus, Eq. (23) is a natural generalization of Eq. (1) to the case of N = 1 supersymmetric theories with multiple gauge couplings.

Examples
3.1 N = 1 SQCD + SQED As a simplest example we consider N = 1 SQCD with N f flavors interacting with the electromagnetic field in the supersymmetric way. This theory is based on the group G × U (1) and contains N f quark flavors in a certain irreducible representation R of the group G. (For "true" SQCD G = SU (3), N f = 6, and R is the fundamental representation.) Each flavor is composed from two chiral matter superfields φ a and φ a in the representations R andR, respectively, with the opposite U (1) charges q a and −q a . At the classical level the action in the massless limit is given by the expression where V and V are the gauge superfields corresponding to the subgroups G and U (1), respectively, and the supersymmetric gauge superfield strengths are defined as Evidently, there are two coupling constants in the theory. The analog of the strong coupling constant is α s ≡ g 2 /4π, while the analog of the electromagnetic coupling constant is α ≡ e 2 /4π. We will denote the corresponding β-functions by β s (α s , α) and β(α s , α), respectively. To construct NSVZ expressions for both these β functions, we need the values T aK defined by Eq. (22). It is easy to see that Substituting these expressions into Eq. (23) and using the notation C 2 = C 2 (G) we obtain the exact expressions for the β-functions corresponding to the strong and electromagnetic coupling constants, β(α s , α) Certainly, they are valid only under certain renormalization prescriptions, which will be discussed below and presumably include the HD+MSL scheme.
Earlier the theory (26) was considered in Refs. [90,91] (see also [102]), where the corrections to the electromagnetic coupling constant generated by loops of quarks, gluons, and their superpartners were analyzed. These corrections are encoded in the Adler D-function [94], for which an all-order expression has been obtained. The D-function calculated in [90,91] is very similar to the β-function for the electromagnetic coupling constant. Up to a normalization factor, the Adler function can be obtained from the function β(α s , α)/α 2 by taking the limit α → 0, Note that in the limit α → 0 all chiral matter superfields have the same anomalous dimension, Really, the anomalous dimensions of various flavors are different because the electromagnetic charges q a are different for different flavors. In the limit α → 0 it is necessary to omit all superdiagrams containing propagators of the electromagnetic gauge superfield V , so that the dependence on the electromagnetic charges disappears. Therefore, using Eq. (32) the exact Adler function (31) can be rewritten in the form This expression reproduces the result of Refs. [90,91] in the particular case when G = SU (N ) and R is the fundamental representation (so that dim R = N ). 5 Also Eq. (33) agrees with the result of Ref. [102] if we take into account that for an irreducible representation γ i j = δ j i · γ. A detailed analysis under what renormalization prescriptions Eq. (33) is valid has been done in Ref. [52]. In particular, it holds in the HD+MSL scheme [102].

NSVZ relations for MSSM
As another illustration, we consider NSVZ equations for MSSM. This model is a softly broken N = 1 supersymmetric theory with the gauge group Certainly, this implies that there are three gauge coupling constants ; Due to the factor 5/3 in the definition of α 1 the gauge coupling unification condition is written in the simple form α 1 = α 2 = α 3 . Actually, this factor encodes the value of the Weinberg angle at the scale of Grand Unification, sin 2 θ W = 3/8. The chiral matter superfields of MSSM include three generations of quarks and leptons and two Higgs douplets The model under consideration does not contain right neutrinos. This is not essential for us, because they acquire large masses of the GUT scale order and do not affect running of couplings at low energies. Therefore, each generation includes the chiral superfields U a , D a , E, The right upper quarks, right down quarks, and right charged leptons are components of U a , D a , and E, respectively. The superfields U a and D a lie in the fundamental representation of SU (3), while E-s are SU (3) singlets. All superfields containing right quarks and leptons are SU (2) singlets. The (charge conjugated) left quarks and leptons are components of Q a and L, respectively. These superfields transform under the fundamental representation of SU (2). With respect to SU (3) the superfields Q a and L belong to the antifundamental and trivial representations, respectively. All chiral matter superfields nontrivially transform under U (1). The corresponding charges (i.e., hypercharges) are presented in Table 1 together with the SU (3) and SU (2) representation. The subscripts 1, 2, and 3 in this table correspond to the generations. The index a numerates the chiral matter superfields in MSSM that transform under irreducible representations of SU (3) and SU (2).
Values of T aK for all a and K = SU (3), SU (2), U (1) calculated according to the definition (22) are presented in Table 2. Note that for quarks and leptons they correspond to the superfields of a single generation. Using these values it is easy to calculate the (well-known) coefficients which are present in the one-loop MSSM β-functions, However, we are interested in exact all-order expressions for the β-functions, which can be constructed with the help of Eq. (23). Using the data of Table 2 from this equation we conclude that they can be written as where the index I numerates generations. Note that the factor 3/5 in the last equation appears due to the factor 5/3 in the definition of α 1 , see Eq. (35). Certainly, Eqs. (41) - (43) are in agreement with the analogous equations written in [84] in a different form. For completeness, in Appendix A.1 we verify that these equations reproduce the correct two-loop expressions for the MSSM β-functions (which was first made in [85] starting from the result of Ref. [84]). For this purpose in Appendix A.1 we listed the one-loop contributions to the anomalous dimensions entering these equations and substitute them into the NSVZ relations. The result of this calculation exactly coincides with the known expressions for the two-loop βfunctions, see, e.g., [104]. Certainly, the coincidence occurs due to the scheme independence of the two-loop contributions to the gauge β-functions, which will be discussed below in Sect. 4.1.

NSVZ relations for the flipped SU(5) model
As another interesting example we consider the so-called flipped SU (5) model [95][96][97][98]. This model is based on the gauge group SU (5)×U (1) and, therefore, contains two coupling constants. We will denote them by α 5 and α 1 . Certainly, the normalization of α 1 should be fixed by a certain condition, which will be described below.
The matter content of the flipped SU (5) theory includes three generations of quarks and leptons which are placed into the representations In our notation a first symbol denotes an SU (5) representation, while a U (1) charge is written in brackets. Unlike the usual SU (5) model, in the flipped SU (5) model the charged leptons are singlets with respect to SU (5), the representation 10 contains D and N (instead of U and E, respectively), and the representation 5 contains U (instead of D), The SU (5) × U (1) symmetry can be broken down to SU (3) × SU (2) × U (1) Y by the vacuum expectation values of two chiral Higgs superfields H and H in the representations 10(−1) and 10(1), respectively. The group U (1) Y is obtained as a superposition of the SU (5) transformations with (where 1 3 and 1 2 are 3×3 and 2×2 identity matrices, respectively) and the U (1) transformations with where q is a U (1) charge of the corresponding superfield normalized as in Eq. (44). Also it is necessary to include Higgs superfields responsible for breaking SU (2) × U (1) Y down to U (1) em . This is done by introducing two chiral superfields h and h in the representations 5(2) and5(−2), respectively. Moreover, (for N G = 3 generations) the model includes N G + 1 = 4 chiral superfields φ R which are singlets with respect to SU (5) and have the U (1) charges equal to 0. The superpotential of the resulting theory is presented in Appendix A.2.
Note that the normalization of the U (1) charge adopted in Eq. (44) is not convenient for physical applications. To construct a more convenient prescription, we recall that the considered model can naturally be embedded into the SO(10) theory, which has a single coupling constant. Therefore, a better requirement for fixing a normalization of α 1 is obtained by the requirement that the gauge couplings are unified as α 1 = α 5 . To find the relation between these two different normalization conditions, we note that the generators of the SU (5) subgroup embedded into SO(10) and normalized by Eq.
In this equation t A r and t A i (with A = 1, . . . , 12) denote real and purely imaginary generators of the SU (5) fundamental representation, respectively. This implies that e 5 = e 10 / √ 2. The generator of the U (1) subgroup embedded into SO(10), again, normalized by Eq. (3) is written as where 1 5 is the 5 × 5 identity matrix. The image of this generator in the spinor representation (see, e.g., [1]) multiplied by e 10 acting on the SU (5) representations 10, 5, and 1 gives respectively. This implies that if normalize α 1 by the condition α 1 = α 5 for the SO(10) symmetric theory, then the U (1) coupling constant should be defined as (1)  10  45  25  10  20  0   Table 3: Values of T aK (where K = SU (5), U (1)) for the chiral matter superfields of the flipped SU (5) model.
The coefficients T aK for the flipped SU (5) model are presented in Table 3. Note that for the representation 10 the generators are written as where (t A ) i j are the generators of the fundamental representation, so that Certainly, the number 3/2 here and similar numbers for other representations are well-known (see, e.g., [105]), but the explicit calculation is more visual.
Substituting the values T aK presented in Table 3 into Eq. (23) and taking into account the factor 40 in Eq. (50) we obtain the exact β-functions for the flipped SU (5) model in the form (Certainly, these equations are written for the region with unbroken SU (5) × U (1) symmetry, i.e., above all thresholds.) In the one-loop approximation Eqs. (53) and (54) exactly agree with the results of Ref. [98]. For calculating the two-loop contributions to the β-function we need to know the one-loop expressions for the anomalous dimensions present in these equations. They are collected in Appendix A.2. Substituting them into the exact NSVZ relations after some transformations we obtain the two-loop β-functions, which are also presented in Appendix A.2. Their parts which do not contain Yukawa couplings agree with the corresponding expressions presented in [98]. (Unfortunately, we did not find the expressions for the parts containing the Yukawa couplings for this model.) Therefore, we obtain a nontrivial evidence that the NSVZ equations are also valid for theories with multiple gauge couplings.

NSVZ relations and the ambiguity of choosing a renormalization prescription 4.1 Scheme dependence of NSVZ relations for theories with multiple gauge couplings
It is well known [42][43][44][45][46] that the NSVZ relation does not hold for an arbitrary renormalization prescription, and even the most popular DR scheme is not NSVZ. The scheme dependence of the NSVZ equation was analysed in [106] for theories with a single (gauge) coupling constant. For theories with simple gauge groups containing Yukawa couplings a similar analysis has been done in [49]. In this paper we extend these results to the case of theories with multiple gauge couplings.
It is known [107] that various renormalization schemes are related by finite renormalizations, which can be written in the form where α ′ (α, λ), λ ′ (α, λ), and z a (α, λ) are finite functions. Note that here we will not assume that the finite renormalization of Yukawa couplings is related to the finite renormalization of matter superfields. Under the finite renormalization (55) the β-functions change as [108,109] where the summation over the index L is performed from 1 to n (see Eq. (6)), and the Yukawa β-function is defined as For theories with a single gauge coupling the two-loop β-function is scheme independent [109]. For multicharge theories containing Yukawa couplings the Yukawa β-functions are schemedependent starting from the two-loop approximation. However, the gauge β-functions even for theories with multiple gauge couplings are scheme-independent in the two-loop approximation [104]. To prove this statement carefully, we should take into account that finite renormalizations should have the same structure as the corresponding quantum corrections [42,50,110].
In particular, this implies that the new gauge couplings α ′ K are related to the original ones by the equations where C K are some constants. Consequently, the lowest term in the expression for α ′ K contains only α 2 K , Substituting this expression into Eq. (56) and taking into account that in the one-loop approximation β K ∼ α 2 K and that the one-loop expression for (β λ ) ijk is given by a sum of terms proportional to αλ and λ 3 we obtain In the lowest approximation each β-function can be presented in the form where β 0K are numerical constants and β 2-loop,K (α, λ) contains terms proportional to α 3 and α 2 λ 2 . Substituting the expression (61) into Eq. (60) we see that both sides of the resulting equation have the same dependence of the constants C K , From this equation we conclude that the two-loop gauge β-functions are scheme-independent even for theories with multiple gauge couplings, 6 Under the finite renormalization (55) the anomalous dimension of chiral matter superfields changes as Evidently, in the one-loop approximation the anomalous dimension is scheme-independent. This implies that the NSVZ equations (23) do not depend on a renormalization prescription up to the order O(α, λ 2 ) inclusive exactly as for theories with a single gauge coupling constant. Next, let us construct an equation which describes how the NSVZ equations change under the finite renormalization (55). For this purpose we first write the transformations inverse to (56) and (64), which are obtained by making the replacement α ′ K ↔ α K , λ ′ ↔ λ, z a → (z a ) −1 . We will assume that original RGFs (without primes) satisfy the NSVZ equation (23). In this equation we express them in terms of the new RGFs (denoted by primes). This gives the equation Solving it for β ′ K (α ′ , λ ′ ) we obtain the equation relating new RGFs, where the summation over L is performed from 1 to n, and we use the notations For theories with a single gauge coupling constant this equation agrees with the results of Refs. [49,106]. From Eq. (66) we see that under arbitrary finite renormalizations the NSVZ equation does not preserve its form. However, as we will see in the next section, there is a certain class of finite renormalizations which do not break the form of the NSVZ relation.

A class of the NSVZ schemes
Even for N = 1 SQED there is a continuous set of the NSVZ renormalization schemes [39]. For non-Abelian theories with a single coupling a similar class of the NSVZ schemes was obtained in [40]. In this section we construct a generalization of these results to the case of theories containing multiple gauge couplings and Yukawa couplings. For this purpose we rewrite Eq. (23) in the equivalent form Using the definitions of the β-functions and the anomalous dimensions of matter superfields (Eqs. (16) and (18), respectively) this equation can be presented in the form [8] d d ln µ Two different renormlization schemes can always be related by a finite renormalization [107], which for the considered theory can be written in the form (55). Let us assume that the NSVZ equations (23) are valid for both original and new couplings. This implies that Eq. (70) also holds for the new couplings. Therefore, they satisfy the equation Subtracting Eq. (70) from Eq. (71) we obtain the equation which contains only functions which describe the considered finite renormalization, We will be interested only in trivial solutions of this equation, which are obtained if the expressions in the curly brackets (numerated by the index K) are constants and do not depend on couplings. Denoting them by B K we derive the system of equations describing finite renormalizations which do not take out of the class of the NSVZ schemes for the considered theory, Evidently, this class is parameterized by m arbitrary functions z a and n arbitrary constants B K , where m is the number of the chiral matter superfields φ a , and n is the number of factors in the product (6). Note that if B K defined as the left hand side of Eq. (73) are not constants, then the matrices M K L and P L (given by Eqs. (67) and (68), respectively) can be written as From these equations and Eq. (66) it is evident that the NSVZ equations really remain invariant under such finite renormalizations for that B K are independent of couplings. For theories with a single gauge coupling constant a certain subset of NSVZ schemes is produced by the HD+MSL prescription, when a theory is regularized by higher covariant derivatives, and divergences are removed by minimal subtractions of logarithms. Presumably, this is also so for theories with multiple gauge couplings. Note that the HD+MSL prescription is not unique, because minimal subtractions of logarithms can supplement various versions of the higher covariant derivative regularization. The explicit dependence of the two-loop anomalous dimension and of the three-loop β-function on regularization parameters for a general N = 1 supersymmetric gauge theory with a single gauge coupling in the HD+MSL scheme can be found in Ref. [112]. This dependence is not trivial, so that HD+MSL schemes really give a continuous set of NSVZ schemes, which, in general, constitute a certain subset of the whole NSVZ scheme class. 5 The higher covariant derivative regularization for theories with multiple gauge couplings 5.1 General features of the higher derivative regularization in supersymmetric theories with multiple gauge couplings For theories with a single gauge coupling constant some NSVZ schemes are given by the HD+MSL prescription in all loops [63][64][65]. Possibly, this is also so for theories with multiple gauge couplings. That is why it is expedient to construct the higher covariant derivative regularization for such theories. Very likely, if it is supplemented by minimal subtractions of logarithms, then Eq. (23) is valid in all orders of the perturbation theory.
The versions of the higher covariant derivative regularization which were used for supersymmetric theories in Refs. [58,59] can be introduced only for theories satisfying the anomaly cancellation condition tr(T A {T B , T C }) = 0 [72]. For theories considered in Sect. 3 this condition really holds, so that the higher covariant derivative regularization can be constructed. Although it is made in a rather standard way, some features specific for the models under consideration should be described in more detail.
To introduce the higher covariant derivative regularization, first, we need to add a higher derivative term to the classical action [53,54]. Next, one should insert into the generating functional the Pauli-Villars determinants, which cancel one-loop divergences that survived after this [57]. For quantizing the theory we will use the background field method [113][114][115] in the supersymmetric formulation [14,99]. Also it is necessary to take into account that the quantum gauge superfield is renormalized nonlinearly [116][117][118], and that the nonlinear renormalization is really needed for cancelling certain divergences [119,120] and for the renormalization group equations to be satisfied [121]. In principle, it can be reduced to a linear renormalization of an infinite set of parameters (which are very similar to the parameter in the gauge fixing term) included into a certain function F(V ). Then the background superfield method and the nonlinear renormalization are introduced by making the substitutions where F K (V K ) = e 0K F K (V K ) A K T A K and T A K are the generators of a relevant representation. Note that for Abelian gauge superfields there is no need to introduce the non-linear renormalization, so that in this case F K (V K ) = V K . However, if a subgroup G K is non-Abelian, then the corresponding function F K (V K ) is not linear. For theories with a single gauge coupling constant its explicit form in the lowest nontrivial approximation can be found in [119,120]. Taking into account that gauge superfields corresponding to different subgroups G K of the gauge group (6) commute, the resulting replacement for the gauge superfield (12) can be written as After making this replacement and adding terms containing higher derivatives the gauge part of the action corresponding to a subgroup G K is modified as where the gauge superfield strength in the right hand side is defined by the equations The expression (78) contains the regulator function R(x) rapidly increasing at infinity and satisfying the condition R(0) = 1. For simplicity, we use the same regulator function for all subgroups of the gauge group, although this is of course not necessary. In our notation the covariant derivatives present inside this function are defined as ∇ a = D a ;∇ȧ = e 2F (V ) e 2VDȧ e −2V e −2F (V ) .
The parameter Λ in the argument of the function R has the dimension of mass and plays the role of an ultraviolet cutoff. Similarly, for all matter superfields it is necessary to make the substitution where F (x) is a regulator function which should rapidly grow at infinity and satisfy the condition F (0) = 1. As in the gauge part of the action, for simplicity, we use the same regulator function for all chiral matter superfields.
In the superpotential we do not introduce any functions with higher derivatives. The regularized action obtained after the modifications described above will be denoted by S reg .
The gauge fixing terms should be added for all subgroups G K , where K(x) is another regulator function, ξ 0K are bare gauge parameters, and the background covariant derivatives are written as For non-Abelian subgroups G K it is also necessary to introduce the Faddeev-Popov and Nielsen-Kallosh ghosts. Their actions are given by expressions It is important that the replacement S → S reg removes divergences only beyond the oneloop approximation. A characteristic feature of the higher derivative regularization [122] is the presence of residual one-loop divergences, which should be regularized by inserting the Pauli-Villars determinants into the generating functional [57].
Following [58,59], to cancel the one-loop (sub)divergences originating from gauge and ghost loops, for each non-Abelian subgroup G K of the gauge group (6) we introduce three (commuting) chiral superfields ϕ 1,K , ϕ 2,K , and ϕ 3,K in the adjoint representation. Then the one-loop divergences will be regularized if we insert into the generating functional the factor where Note that in this case (88) Similar equations are also valid for the other superfields ϕ.
The (sub)divergences generated by a matter loop are regularized by inserting into the generating functional where c K are certain numbers which will be specified below (see Eq. (92)). For the non-Abelian subgroups G K it is possible to choose the commuting chiral Pauli-Villars superfields Φ K in the adjoint representation of the subgroup G K . With respect to the other subgroups these superfields are singlets. Then, for regularizing one-loop matter (sub)divergences the actions for them should be taken in the form If a subgroup G K coincides with U (1), then we introduce two commuting chiral Pauli-Villars superfields Φ K and Φ K with the charges ±1 with respect to the considered U (1) subgroup. Again, with respect to the other subgroups they are invariant. The action for these superfields is given by the expression The cancellation of the one-loop divergences and subdivergences occurs if the constants c K in Eq. (89) have the values for non-Abelian subgroups G K ; (Note that the factor 1/2 appears because in the Abelian case we introduce two Pauli-Villars superfields.) Also it is important that the masses of the Pauli-Villars superfields should be proportional to the dimensionful parameter Λ present in the higher derivative terms, the coefficients a ϕ,K and a K being independent of couplings.

The higher covariant derivative regularization for MSSM
As an example, we construct a version of the higher covariant derivative regularization for MSSM. In this case the regularized action is constructed exactly according to the prescriptions (78) and (81) described in the previous section. That is why below we will not discuss this construction in detail. As for the Pauli-Villars determinants, to cancel the one-loop divergences coming from gauge and ghost loops, we insert into the generating functional Here we introduced three commuting Pauli-Villars superfields for both non-Abelian subgroups (SU (3) and SU (2)) of the gauge group SU (3) × SU (2) × U (1). Evidently, for the U (1) subgroup there is no need to introduce such superfields. According to Eq. (87) the action for them is written as The corresponding masses are proportional to the parameter Λ in the higher derivative terms, where the coefficients a ϕ,3 and a ϕ,2 are constants independent of couplings.
To cancel the one-loop (sub)divergences which appear from a loop of chiral matter superfields, we introduce the chiral Pauli-Villars superfields Φ 3 , Φ 2 , Φ 1 and Φ 1 . The superfield Φ 3 lies in the adjoint representation of SU (3) and is invariant under SU (2) and U (1) transformations. Similarly, the superfield Φ 2 lies in the adjoint representation of SU (2) and remains invariant under SU (3) and U (1) transformations. The superfields Φ 1 and Φ 1 are singlets with respect to SU (3) and SU (2), but have nontrivial (±1 in units of e 01 ) opposite U (1) charges. The gauge invariant actions for these superfields have the form where and the coefficients a 1,2,3 do not depend on couplings. The (inverse) Pauli-Villars determinants are given by the expressions According to Eq. (92), to find the degrees of these determinants in the generating functional, we need the coefficients T K (R) defined by Eq. (21) and the constants C 2 (G K ) for the non-Abelian subgroups G K . In the case of MSSM they have the following values: Therefore, using Eq. (92) we conclude that the generating functional for the regularized theory is given by the expression (103) where S reg is the regularized action, which in particular includes the higher derivative terms, S gf is the gauge fixing action, S FP and S NK are actions for the Faddeev-Popov and Nielsen-Kallosh ghosts, respectively, S ϕ is the Pauli-Villars action (95), and S sources contains all relevant sources.

The higher covariant derivative regularization for the flipped SU(5) GUT
As another example we consider the flipped SU (5) model (for which the gauge group is SU (5) × U (1)). The regularized action is again constructed according to the prescriptions (78) and (81)  To cancel the one-loop divergences generated by the gauge and ghost superfields, we use three chiral superfields ϕ 1,2,3 in the adjoint representation of SU (5) neutral with respect to U (1) and insert into the generating functional the corresponding (inverse) determinant where The divergences coming from a loop of chiral matter superfields are regularized with the help of the (commuting) Pauli-Villars superfields Φ 5 , Φ 1 , and Φ 1 . The chiral superfield Φ 5 lies in the adjoint representation of SU (5) and is neutral with respect to the U (1) subgroup. The chiral superfields Φ 1 and Φ 1 are SU (5) singlets and have the opposite U (1) charges ±e 01 , where e 01 is normalized by Eq. (44). The actions for these superfields are given by the expressions which are used in constructing the corresponding Pauli-Villars determinents Standardly, the Pauli-Villars masses are proportional to the parameter Λ, Then, according to the prescription described in Sect. 5.1, the generating functional for the regularized theory can be written as 6 Non-renormalization of the triple gauge-ghost vertices and a new form of the NSVZ equation A very important ingredient of the all-loop derivation of the NSVZ β-function for theories with a single gauge coupling is a new form of the NSVZ equation [64]. It relates the β-function in a certain loop to the anomalous dimensions of the quantum gauge superfield, of the Faddeev-Popov ghosts, and of the matter superfields in the previous loop. It is this equation that appears after summing singular contributions [63] produced by integrals of double total derivatives which determine the β-function in the supersymmetric case [65]. In this section we construct such a form of the NSVZ β-function for theories with multiple gauge couplings. This is done using the non-renormalization theorem for the triple gauge-ghost vertices which is proved for these theories with the help of the Slanov-Taylor identities [123,124]. The proof is very similar to the one discussed in [64] for theories with a single gauge coupling, so that here we will describe it more briefly. 7 Note that for theories with multiple gauge couplings we will consider only the triple gauge-ghost vertices in which external lines correspond to the same simple subgroup G K . If the external lines correspond to different subgroups, then the divergent contributions are forbidden by the renormalizability.

Slavnov-Taylor identities for the triple gauge-ghost vertices in theories with multiple gauge couplings
In the case of using the background superfield method the original gauge invariance of a classical theory generates two types of transformations. The background gauge invariance remains a manifest symmetry of the effective action, while the quantum gauge invariance is broken down to the BRST transformations [128,129] by the gauge fixing procedure. In the supersymmetric case the BRST transformations are written as [116,130,131] where ε K are anticommuting parameters which do not depend on the coordinates. Note that for the gauge and matter superfields the BRST transformations are reduced to the quantum gauge transformations for which the parameters corresponding to non-Abelian subgroups G K are given by the chiral superfields This implies that any gauge invariant expression is also BRST invariant. The BRST invariance leads to the Slavnov-Taylor identities [123,124], which are a key ingredient for proving the finiteness of the triple gauge-ghost vertices. To obtain these identities, we make a change of variables coinciding with the transformations (112) corresponding to the subgroup G K in the generating functional Z. The total action of the theory is invariant under the transformations (112). By construction, the Pauli-Villars determinants are invariant under both background and quantum gauge transformations. As a consequence of the latter invariance, they are also BRST invariant. Therefore, only the source term (115) in the generating functional is not BRST invariant. Using the standard technique (see, e.g., [122]) a set of the generating Slavnov-Taylor identities in the considered case can be written as Each of these identities corresponds to a certain non-Abelian subgroup of the gauge group. The (super)fields in Eq. (116) are not set to 0, so that it is possible to differentiate these equalities and construct an infinite set of identities relating various Green functions. For deriving the non-renormalization theorem we need the identity which is obtained after differentiating Eq. (116) with respect toc +B K y , c C K z , and c D K w . After this we set all superfields to 0. In the theories with multiple gauge couplings the ghost number conservations hold for each simple subgroup G K of the gauge group G. Therefore, the numbers of ghosts and antighosts corresponding to any non-Abelian G K should be equal in each nontrivial Green function. Taking this into account we obtain the Slavnov-Taylor identity (Terms containing matter Green functions vanish due to the global Z 3 -symmetry φ j → e 2πik/3 φ j of the massless theory.) The Green functions entering this equation can be expressed in terms of certain functions depending on external momenta. To construct the corresponding expressions, we need to involve chirality and symmetry considerations. For instance, taking into account that the ghost superfields are chiral or antichiral we can write the corresponding two-point Green functions in the form where the functions G c K depend on couplings and ∂ 2 /Λ 2 . Evidently, all two-point functions in which ghosts correspond to different subgroups vanish due to the symmetries responsible for the ghost number conservations. Note that in theories with multiple gauge couplings the functions G c K are in general different for different K.
The derivative of δV A K x with respect to the ghost superfield can also be expressed in terms of the function G c K . For this purpose we will use the identity obtained with the help of the substitutionc +B K →c +B K + a +B K (where a +B K is an arbitrary chiral superfield) in the generating functional Z. In terms of the effective action the result can be written as the identity Again, in this equation (super)fields are not set to 0, so that it can be differentiated. This implies that Eq. (119) generates an infinite set of identities relating various Green functions. In particular, if we differentiate Eq. (119) with respect to the ghost superfield c A K y and take into account chirality consideration, then the result can be written in the form Next, using chirality considerations we write down all possible structures that can appear in the part of the effective action corresponding to the triple Vc + c vertex, where f A K B K C K denotes structure constants of the subgroup G K . The functions s K , S µ,K , and S K have the dimensions m −2 , m −1 , and m 0 , respectively. Note that in [64] similar functions were denoted by f , F µ , and F , but here (following Ref. [72]) we use other letters in order to distinguish them from the higher derivative regulator function F , see Eq. (81). From Eq. (121) we obtain the explicit expression for the three-point Vc + c vertex, where we introduced the notation Also Eq. (117) contains the second derivative of the expression with respect to the ghost superfields. It is convenient to introduce the auxiliary source term which is invariant under the BRST transformations due to their nilpotency. This term allows presenting the considered correlator as the third derivative of the effective action with respect to the source J and two ghost superfields. Involving chirality considerations the corresponding contribution to the effective action can be written in the form where H K (p, q) are dimensionless functions. Due to the symmetry with respect to permutations of the ghost superfields they satisfy the equations Then the correlator entering Eq. (117) can be expressed in terms of these functions, Substituting the expressions (118), (120), (122), and (128) into the Slavnov-Taylor identity (117) after some transformations we rewrite this equation in the equivalent form where, for simplicity, the argument of the functions G c is denoted by q instead of −q 2 /Λ 2 . We see that for theories with multiple gauge couplings the number of such equations is equal to the number of non-Abelian subgroups G K in the direct product (6). Eq. (129) is the most convenient form of the Slavnov-Taylor identity for the Vc + c vertex. In the next section it will be used for proving its finiteness in all orders of the perturbation theory.

All-loop finiteness of the triple gauge-ghost vertices
To prove the all-loop finiteness of the triple gauge-ghost vertices in which all external lines correspond to the same non-Abelian subgroup G K , we first prove the finiteness of the functions H K defined by Eq. (126). Each of these functions is contributed by superdiagrams in which one external line corresponds to the source J , and the other two correspond to the Faddeev-Popov ghost c K . The vertex with the external source is given by the expression (125). If both ghost superfields in Eq. (125) produce propagators, then the supergraph contributing to the function H K contains where two ghost propagators connect the superspace point y to the superspace points 1 and 2.
With the help of the identity the expression (130) can be presented in the form of the integral over the full superspace, Evidently, if (only) one ghost superfield in the vertex (125) corresponds to an external line, then the vertex can also be presented as such an integral. (Two external ghosts in Eq. (125) can appear only in the tree contribution to the function H K .) This implies that all vertices in the considered supergraph can be written as integrals over the full superspace. Therefore, it is possible to apply standard rules for calculating supergraphs. It is well known (see, e.g., [100]) that the resulting expression is also given by an integral over the full superspace. However, all external lines (which correspond to J K and two c K ) are chiral. Therefore, due to Eq. (131) a nontrivial result can be obtained only if at least two right spinor covariant derivatives act on the external lines. If we take into account that the function H K is dimensionless, this implies that the dimension of the remaining loop integral is m −2 . In other words, the degree of divergence is equal to (−2), so that this integral is finite in the ultraviolet region. Due to the renormalizability of the considered theory possible subdivergences are removed after the renormalization in the previous orders. Thus, the functions H K (for any K such that G K = U (1)) appear to UV finite. Next, we involve the Slavnov-Taylor identity (129) multiplied by the ghost renormalization constant Z c K . Using the finiteness of the renormalized Green function we see that Therefore, differentiating Eq. (129) multiplied by Z c K with respect to ln Λ and keeping in mind the finiteness of the function H K we obtain In this equation we set p = −q and take into account that G c K (−q) = G c K (q) because this function depends only on q 2 /Λ 2 . The result can be written in the form The function S K (−q, q) also depends only on q 2 /Λ 2 , so that from the above equation we conclude that it is finite in the ultraviolet region. In our notation the renormalization constants for the gauge coupling constants, for the quantum gauge superfields, and for the Faddeev-Popov ghosts are defined by the equations respectively. Taking into account that we see that the components of the quantum gauge superfields and of the Faddeev-Popov ghosts are renormalized as Using these equations one can obtain the renormalization constants for the triple gauge-ghost vertices from Eq. (121). Certainly they are the same for all 4 types of the considered vertices and for all functions entering Eq. (121). For the vertex corresponding to the subgroup G K the renormalization constant is equal to Z −1/2 The renormalized function in the left hand side of this equation is UV finite by definition. Certainly, this is so for the particular case p = −q. From the other side, as we discussed above, the function S K (−q, q) is not divergent in the ultraviolet region. Therefore, the corresponding renormalization constant should be finite, This implies that it is possible (although not necessary) to choose a renormalization scheme in which the renormalization constants satisfy the relation Evidently, this equation holds in the HD+MSL scheme, because in this case various Z-s include only powers of ln Λ/µ, while all finite constants vanish.

A new form of the NSVZ β-function for theories with multiple gauge couplings
Let us derive a new form of the NSVZ equation for theories with multiple gauge couplings making the transformations analogous to the ones proposed in Ref. [64]. However, here we will use RGFs defined in terms of the renormalized coupling assuming that Eq. (143) is satisfied. Certainly, for RGFs defined in terms of the bare couplings a similar result is also valid. (In this case it is sufficient to use Eq. (142), while Eq. (143) is not needed.) First, we equivalently rewrite the NSVZ equation (23) in the form It is convenient to express the function β K (α, λ) in the right hand of this equation in terms of the renormalization constant for the gauge coupling α K , This renormalization constant can in turn be related to the renormalization constants of the quantum gauge superfield and of the Faddeev-Popov ghosts corresponding to the subgroup G K with the help of equation (143). As a result, the β-function for a certain gauge coupling can be expressed in terms of the anomalous dimensions for the corresponding quantum gauge superfield and Faddeev-Popov ghosts, where Substituting the expression (146) for the function β K into the right hand side of Eq. (144) we obtain the equivalent form the NSVZ relation which generalizes a similar equation constructed in [64] for theories with a single gauge coupling constant.
Unlike the original NSVZ relation (23), this form of the NSVZ equation has a simple graphical interpretation. If we consider a vacuum supergraph (which does not contain external lines), then the corresponding contribution to the function β K is produced by the sum of all superdiagrams which are obtained by attaching two external lines of the background gauge superfield V K . From the other side, cutting internal lines we produce superdiagrams contributing to anomalous dimensions of various quantum superfields. Eq. (148) relates these contributions to the βfunction. 8 However, it is worth to note that for the quantum gauge superfields only cuts of the V K -propagators contribute to the function β K . Possibly, this can be explained with the help of the method for constructing integrals of double total derivatives proposed in [65]. We also hope that the generalization of the results obtained in [63] allows deriving Eq. (148) in a way analogous to the case of theories with a single gauge coupling constant.

Conclusion
In this paper we investigated the NSVZ equations for N = 1 supersymmetric models with multiple gauge couplings. In these theories the gauge group is a direct product of simple subgroups and U (1) factors. Certainly, there are a lot of various phenomenologically interesting examples of such theories, for instance, MSSM or flipped SU (5) GUT considered in this paper in detail. In theories with multiple gauge couplings a number of β-functions is equal to a number of (simple and U (1)) factors in the gauge group. We assume that in all loops each of these βfunctions satisfies the NSVZ equations (23) and (148) in a certain subtraction scheme. A strong evidence in favor of this is that in MSSM and flipped SU (5) GUT predictions of the NSVZ equations for the two-loop β-functions exactly coincide with the corresponding (well-known) results obtained by direct calculations made in DR-scheme. Certainly, this coincidence occurs because the two-loop gauge β-functions are scheme independent. 9 From the NSVZ equations for theories with multiple gauge couplings it is also possible to derive the exact NSVZ-like equation for the Adler D-function in N = 1 SQCD, which was first obtained in [90,91].
It is highly probable that the NSVZ scheme for theories with multiple gauge couplings is given by the HD+MSL prescription as for theories with a simple gauge group [63]. HD means that a theory should be regularized by Higher covariant Derivatives. This regularization in particular includes insertion of the Pauli-Villars determinants for removing one-loop divergences. Note that the ratios of the Pauli-Villars masses to the dimensionful regularization parameter (Λ) should be independent of couplings. MSL indicates that divergences are removed with the help of Minimal Subtractions of Logarithms, when only powers of ln Λ/µ are included into renormalization constants. Therefore, for investigating the NSVZ equations it is highly desirable to construct a version of the higher covariant derivative regularization applicable for theories multiple gauge couplings. This was also done in this paper for the general case and was illustrated by the examples of MSSM and flipped SU (5) GUT. Also the higher covariant derivative regularization is a very important ingredient needed for the all-loop perturbative derivation of the NSVZ β-function. According to [63][64][65][66], summing the perturbation series we do not obtain the original form of the NSVZ equation. Instead of it we obtain the relation between the Lloop β-function and the (L − 1)-loop anomalous dimensions of the quantum gauge superfield, of the Faddeev-Popov ghosts, and of the matter superfields. The original NSVZ equation can be derived from it with the help of the non-renormalization theorem for the triple gauge-ghost vertices. In this paper we generalized this theorem to the case of theories with multiple gauge couplings and obtained from it the new form of the NSVZ equations for theories with multiple gauge couplings, which is given by Eq. (148). We hope to finalize the all-loop perturbative derivation of these equations in the forthcoming publications.

A NSVZ equations in the order O(α, λ 2 )
In this appendix we present expressions for the (scheme-indepedent) one-loop anomalous dimensions for all chiral matter superfields in MSSM and flipped SU (5) GUT. We substitute them into the NSVZ equations and obtain two-loop expressions for the β-functions. As we discussed in Sect. 4.1, they do not depend on a renormalization prescription, so that it is possible to compare them with the known results obtained in the DR-scheme. Below we demonstrate the agreement of the results derived by these two methods. (For MSSM this was first done in [85].) This agreement can be considered as a strong evidence that NSVZ equations discussed in this paper are really valid for theories with multiple gauge couplings.

A.1 MSSM
Let us construct two-loop contributions to the β-functions of MSSM starting from the exact β-functions (41) - (43). For this purpose we need expressions for the one-loop anomalous dimensions of all chiral matter superfields. For a theory with a simple gauge group the one-loop anomalous dimension is given by the expression To generalize this result to the case of theories with multiple gauge couplings, we first discuss the structure of the matrix λ * imn λ jmn . It is evidently Hermitian and invariant under the transformations of the gauge group. If matter superfields lie in an irreducible representation of a simple gauge group, then this matrix is proportional to δ j i . However, in the general case it is not so. Nevertheless, making a proper unitary rotation in the generation space it is always possible to diagonalize the matrix λ * imn λ jmn . We will always assume that this is done. When this matrix becomes diagonal, the superfields φ a are its eigenvectors corresponding to the certain eigenvalues, which we will denote by (λ * λ) a . Then we see that for theories with multiple gauge couplings the one-loop expression for the anomalous dimension of the superfield φ a generalizing Eq. (149) can be written in the form where the factor C(R aK ) is defined by Eq. (19). These factors can be calculated as follows. If R is an irreducible representation of a simple group G, then where r is a dimension of the group G and dim R is a dimension of the representation R. In particular, for the (anti)fundamental representation of the group SU (N ) For the U (1) subgroups C(R aK ) should be replaced by q 2 a , where q a is a charge of the superfield φ a with respect to the considered U (1) subgroup.
Because Eq. (150) contains the Yukawa couplings, for calculating one-loop anomalous dimensions we also need the explicit expression for a part of the MSSM action containing the superpotential. It is written in the form where the superpotential is given by the expression In this expression the indices I, J = 1, 2, 3 numerate generations of elementary particles. The Yukawa coupling λ ijk can be written in terms of the 3 × 3 Yukawa matrices (Y U ) IJ , (Y D ) IJ , and (Y E ) IJ . The parameter µ (which should be distinguished from the normalization point µ) has the dimension of mass and is not essential for calculating the considered renormalization group functions. Now, it is possible to obtain all expressions for the anomalous dimensions of the chiral matter superfields starting from Eq. (150). Certainly, constructing them we need to take into account the factor 5/3 present in the definition of the constant α 1 , see Eq. (35). The result can be written in form (see, e.g., [86]) Note that in our conventions the matrices Y + U Y U , Y + D Y D , Y + E Y E should be made diagonal by specially constructed unitary rotations of chiral superfields U , D, and E which act on the indices numerating generations. Similarly, the matrices Y U Y + U + Y D Y + D and Y E Y + E should be made diagonal by proper unitary rotations of the superfields Q and L, respectively. 10 In Eqs. (155) -(161) we certainly assume that for calculating anomalous dimensions of the quark and lepton superfields it is necessary to take a diagonal element of a proper product of the Yukawa matrices which corresponds to the considered generation.

A.2 The flipped SU(5) model
Now let us calculate two-loop β-functions for the flipped SU (5) model by substituting oneloop anomalous dimensions of various chiral superfields into the NSVZ equations. For completeness, here we consider the model containing N G generations of the matter superfields in the representation 10(1) + 5(−3) + 1 (5), N H Higgs superfields in the representation 10(−1) + 10(1), N h Higgs superfields in the representation 5(2) + 5(−2), and N G + 1 singlets φ. The theory considered in Sect. 3.3 corresponds to Below indices of the chiral superfields 10 ij , 5 i , and E which numerate the generations will be denoted by the letters I, J = 1, . . . , N G . Similarly, the Higgs superfields h i and h i will be numerated by the letters K, L = 1, . . . , N h . For the Higgs superfields H ij and H ij we will use the letters M, N = 1, . . . , N H . The singlets will be numerated by the letters R, S, T = 1, . . . , N G + 1.
The superpotential of the model has the form [96] W (φ a ) = (λ 1 ) IJK ε ijklm (10 ij ) I (10 kl ) J ( h m ) K + (λ 2 ) IJK (10 ij ) I (5 i ) J (h j ) K + (λ 3 ) IJK (5 i ) I Note that in our notation the superfields in the representations 10 and 10 are normalized in such a way that their kinetic terms contain the coefficient 1/4. For instance, in the classical action the term quadratic in the superfield H ij (without the gauge superfields) has the form To avoid double summations, one can rescale all superfields in the representations 10 and 10 in 1/ √ 2 times, e.g., H ij → H ij / √ 2, etc. If N G , N H , and N h are not fixed, then from Eq. (23) and Table 3 (and also using Eq. (50)) we obtain the exact β-functions 2γ 10 I (α, λ) + 9γ 5 I (α, λ) + 5γ E I (α, λ) To construct the anomalous dimensions entering these equations, we use Eq.
Using these values after some calculations we obtain the anomalous dimensions for all chiral superfields of the model. In the expressions for them listed below, the summation is not performed over the indices indicated in bold letters,  Substituting these anomalous dimensions into the NSVZ equations (168) and (169) we obtain the two-loop expressions for the β-functions of the model under consideration, − N h α 1 10π + 12α 5 5π + 9 2π 2 (λ * 1 ) IJK (λ 1 ) IJK + 7 16π 2 (λ * 2 ) IJK (λ 2 ) IJK + 19 32π 2 (λ The gauge parts of these expressions exactly coincide with the ones presented in [98]. (N 5 and N 10 used in [98] are related to N h and N H used in this paper by the equations N 5 = 2N h and N 10 = 2N H .) As for the Yukawa parts, we did not manage to find them in earlier papers.