Gauge and Yukawa coupling beta functions of two-Higgs-doublet models to three-loop order

We compute the beta functions for the three gauge couplings and the Yukawa matrices of a general two-Higgs-doublet model in the modified minimal subtraction scheme to three loops. The calculations are performed using Lorenz gauge in the unbroken phase. We discuss in detail the occurence of poles in anomalous dimensions and propose practical prescriptions to avoid them. We provide explicit results for the often used $\mathbb{Z}_2$-symmetric versions of the two-Higgs-doublet model of type I, II, X and Y. Furthermore, we provide the first independent cross-check of the three-loop Yukawa coupling beta functions of the Standard Model.


I. INTRODUCTION
An appealing renormalization scheme for the couplings of the Standard Model of particle physics (SM) and of its extensions is the minimal modified subtraction (MS) scheme. As a consequence, the numerical values of the couplings depend on the renormalization scale μ, which in general is of the same order as the energy scale of the considered process. The values of the couplings at different scales are related by so-called beta functions which in perturbation theory are given as power series in all couplings of the theory.
In the SM there are three gauge couplings (g 1 , g 2 , g s ), the quartic Higgs boson coupling λ, and a Yukawa coupling for each massive fermion, where often only the third generation couplings, y t , y b , and y τ are considered as nonzero. For all couplings the three-loop beta functions have been completed recently: the gauge coupling beta functions have been computed in Refs. [1][2][3], the ones for the Yukawa couplings have been computed in Refs. [4][5][6], and λ has been considered in Refs. [7,8]. Leading terms to the fourloop QCD beta function and the Higgs self-coupling involving the top Yukawa coupling and α s have been computed in Refs. [9][10][11][12], and within QCD the beta function is even known to five loops [13][14][15].
There are a number of two-loop results that can be immediately adapted to a large class of nonsupersymmetric beyond-the-SM theories. In particular, two-loop results for gauge [16], Yukawa [17], and scalar self-couplings [18] have been known since the middle of the 1980s. Furthermore, the three-loop gauge coupling beta function for a simple gauge group has been calculated [19]. In this work we consider the so-called two-Higgs-doublet model (2HDM) and compute the gauge and Yukawa coupling beta functions to three-loop order.
The 2HDMs, where the SM Higgs sector is extended by a second SUð2Þ Higgs doublet, are attractive extensions of the SM. Although simple and probably not realized in nature in its minimal version, 2HDMs nevertheless constitute prototype extensions of the SM, which can be used to study several features of beyond-SM theories. In particular, for a certain choice of parameters it implements the Higgs sector of the minimal supersymmetric Standard Model. Further motivation and several phenomenological applications can be found in the review [20].
The most general 2HDM has many parameters and furthermore several unwanted features such as flavorchanging neutral currents (FCNCs) at tree level. Thus, often additional symmetries are imposed. For example, if CP conservation in the Higgs sector is assumed, one has five physical scalar degrees of freedom that correspond to two scalar, one pseudoscalar, and a charged Higgs boson. In these models, both Higgs doublets acquire vacuum expectation values v 1 and v 2 such that v ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v 2 1 þ v 2 2 p ≃ 246 GeV determines the W AE and Z boson masses in the same way as in the SM. The ratio v 2 =v 1 is denoted by tan β.
The scope of the present work is twofold: First, we provide the first independent cross-check of the three-loop Yukawa coupling beta functions in the SM. In this context it is particularly important to carefully investigate the scheme used for γ 5 in D ≠ 4 dimensions. Note that for the gauge couplings it is possible to choose Green's functions without external fermions. For Yukawa couplings this is not possible anymore. As a second aim, we extend both the gauge and the Yukawa beta functions to a general 2HDM. There is no change in the underlying integrals, which have to be evaluated; however, there are conceptional challenges in connection to the wave function renormalization of the scalar fields.
The remainder of the paper is organized as follows: In the next section we introduce the 2HDM that serves to fix the notation. Section III is devoted to technical details. In particular, we introduce the renormalization constants for the parameters and fields and define the beta functions and anomalous dimensions that we want to compute. The main focus of Sec. IV relies on the proper definition of the renormalization constants such that the anomalous dimensions are finite. We investigate this problem in detail and propose practical solutions. A detailed discussion of the computation of the gauge and Yukawa coupling beta functions is provided in Secs. V and VI, respectively. In these sections we also explain how one can arrive at special versions of the 2HDM and the SM results. Furthermore, we compare the Yukawa beta functions to Ref. [4]. The findings of this paper are summarized in Sec. VII.

II. TWO-HIGGS-DOUBLET MODEL
An extensive discussion of a general 2HDM model can be found in Ref. [20]. For convenience we repeat in the following the features that are important for our calculation.
The additional Higgs doublet leads to an enlarged Yukawa sector that can be written as The sum runs over the two doublets and "H.c." refers to the Hermitian conjugate part. Y u a , Y d a , and Y l a are generic 3 × 3 complex matrices containing the Yukawa couplings, and Q L , L L , u R , d R , and l R represent left-and right-handed quark and lepton fields.Φ ¼ iτ 2 Φ Ã j is the charge conjugated doublet with τ 2 being the second Pauli matrix.
The 2HDM has furthermore a more involved scalar potential that in its general form is given by [21] The parameters m 2 11 , m 2 22 , and λ 1 ; …; λ 4 are real, whereas in general m 2 12 , λ 5 , λ 6 , and λ 7 are complex. This leads to 14 degrees of freedom, 11 of which are physical as can be seen by an appropriate basis choice for Φ 1 and Φ 2 [20].
As in all multi-Higgs-doublet models, the Lagrange densities given in Eqs. (1) and (2) contain FCNCs. For example, the up-type Yukawa matrices Y u 1 and Y u 2 will not be in general simultaneously diagonalizable, and thus neutral Higgs scalars ϕ will mediate FCNCs of the form uu 0 ϕ already at the tree level, where u ≠ u 0 are two different up-type quarks. To avoid FCNCs at tree level [22,23] it is necessary that all fermions with the same quantum numbers couple to one and the same Higgs multiplet. This condition can be satisfied if all quarks couple to just one of the Higgs doublets or the right-handed up-and down-type quarks couple to different Higgs doublets. Depending on whether the right-handed leptons couple to the Higgs doublets in the same manner as the right-handed down-type quarks, or in the opposite way, further two possibilities can be identified. The resulting four models are summarized in Table I. They can be realized by imposing a Z 2 symmetry to the general model. In fact, the type I 2HDM can be obtained by enforcing an additional Z 2 symmetry under which the theory has to be invariant, namely Φ 1 → −Φ 1 and Φ 2 → Φ 2 . The type II 2HDM can be derived via the symmetries The additional discrete symmetries required for the other two models can be derived similarly. Note that the Z 2 symmetries require that In a generic quark basis as given in Eq. (1) the condition for the nonexistence of FCNCs in the up-type (down-type) quark sector is that the Yukawa matrices Y u 1 and Y u 2 (Y d 1 and Y d 2 ) commute [24]. If one of the two Yukawa matrices is zero, as it is actually the case for the four models shown in Table I, this condition is trivially fulfilled.
The most general Lagrange densities in Eqs. (1) and (2) contain several fields with the same quantum numbers that can mix. Therefore, one can rewrite the Lagrangian in terms of the new fields obtained from the original ones by simple basis transformations. In the following we will refer to these transformations as flavor transformations for both fermions and scalars. Obviously, the physical observables do not depend on such redefinitions. They can depend only on quantities that are invariant under arbitrary unitary flavor transformations. Ideally, one would be able to express the fundamental Lagrangian parameters in terms of these invariants. However, some of the Lagrangian parameters in Eq. (1) that do not take into account flavor symmetries are not physical. That is, there are Lagrangian parameters that can be expressed as linear combinations of others. This also means that there is a basis where the unphysical parameters are identically zero; i.e. one can rotate them away via flavor transformations. In other words, any coupling or mixing angle can be expressed in terms of so-called flavor invariants. This statement has been explicitly proven for the Yukawa sector of the SM [25] and for the scalar sector of the 2HDM, for example, in Refs. [26,27]. 1 In this paper, we (re)confirm the findings of [25][26][27] explicitly for the Yukawa sector of the SM and for Z 2 -symmetric 2HDMs through three loops.
The flavor transformations for fermion and scalar fields in Lagrangian Eq. (1) can be summarized as follows: where U Q , U u , and U d are unitary 3 × 3 matrices and U Φ is a unitary 2 × 2 matrix. Under these unitary basis transformations, the gauge and kinetic terms are unchanged and L Y in Eq. (1) is invariant if the Yukawa matrices transform as In a similar manner, one can derive the transformation properties of the parameters in the potential under redefinitions of the scalar fields [20]. One introduces the rank two and four tensors, K ab and λ ab;cd , so that with K ab ¼ K Ã ba ; λ ab;cd ¼ λ cd;ab ; λ ab;cd ¼ λ Ã ba;dc : ð6Þ One can match with the standard notation given in Eq. (2) and obtain the following relations: The two tensors transform under the basis change given in Eq. (3) as Since the calculation of the MS renormalization constants can be performed in the unbroken phase, the dimensionful parameters m ij are irrelevant, and thus for our calculation of the beta functions only the second transformation in Eq. (7) will be of interest.
Within the SM the physical Yukawa couplings are defined via the diagonalization of the Hermitian matrices i.e., Y d 1 is diagonalized via 1 For more details see Ref. [20] and references therein.
However, Eqs. (10) and (13) are in conflict with each other, and only one of the Yukawa matrices can be diagonalized. This leads to the definition of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, which, in the basis where the up-type Yukawa matrix is diagonal, is given by 5 unphysical phases can be eliminated from V via further quark field redefinitions.
The discussion up to now is in analogy to the SM. Within a general 2HDM the unitary transformations discussed above do not necessarily simultaneously diagonalize the other two Yukawa matrices Y u 2 and Y d 2 . We can still define the additional set of (nonphysical) Yukawa couplings as the positive square roots of the eigenvalues of the matrices To summarize, using the unitary rotations in Eq. (10), the set of Yukawa matrices transforms as where N u and N d are complex 3 × 3 matrices. Note that for the special case of a 2HDM with a Z 2 symmetry only two of the four matrices in Eq. (15) are nonzero. Their eigenvalues define the physical parameters, and their mixing matrix is defined in analogy to the CKM matrix in the SM. We want to stress that within the SM and the four Z 2 -symmetric 2HDMs (cf. Table I)  The strategy to construct flavor invariants in the Yukawa sector consists in taking products of Yukawa matrices, contracting over the internal flavor indices, and taking the trace over the external flavor indices. For example, the simplest flavor invariants that can be constructed within a 2HDM read where Tr denotes the trace over the open indices of the lefthanded fermions Q L . In a generic 2HDM the matrices are invariant under scalar flavor transformations, and one can thus construct nine other flavor invariants similar to those for the SM [25,28]. Using Sec. 3. 1 of [25] and adapting the notation (i.e., replacing U and D by M u and M d ) leads to whereM ¼ M −1 detðMÞ. All Z 2 -symmetric 2HDMs have the same 11 invariants as the SM.
In a generic 2HDM further higher rank invariants can be constructed using tensorial properties of the Yukawa matrices. For example, the simplest additional type of rank four invariant tensors are and T ð2Þ dd ¼ and similar ones where Tr is replaced by the determinant. A systematic analysis of all independent invariants for a general 2HDM is, however, beyond the scope of this article. Let us also mention at this point that the definitions for the physical Yukawa couplings and mixing matrices introduced above holds to all orders in perturbation theory.

III. TECHNICALITIES
In this work we compute the beta functions of the three gauge couplings and the Yukawa matrices in the MS scheme.
Our calculation of the beta functions are based on the Lagrange densities in Eqs. (1) and (2). The specification to the types I, II, X, and Y is straightforward. Note that the MS renormalization constants can be computed in the unbroken phase since they do not depend on the particle masses.
It is convenient to denote the gauge couplings by α 1 , α 2 , and α 3 ¼ α s , where α i ¼ g 2 i =ð4πÞ and Y f a with f ¼ u, d, l and a ¼ 1, 2 (labeling the scalar doublets). Furthermore, we introduceλ ab;cd ¼ λ ab;cd =ð4πÞ (a, b, c, d ¼ 1, 2), where λ ab;cd are the quartic coupling in the scalar potential. We define the beta functions via where ϵ ¼ ð4 − dÞ=2. Note that the dependence of the couplings on the renormalization scale is suppressed. The equations defining the beta function for Yukawa matrices have to be understood as matrix equations in flavor space. The gauge couplings are related to the fine structure constant, the weak mixing angle, and the strong coupling as follows: where the SU(5) normalization has been adopted, which leads to the factor 5=3 in the definition of α 1 . For models where the first and second generation Yukawa couplings are neglected, it is convenient to introduce α The beta functions are obtained from the renormalization constants relating bare and renormalized couplings. For the gauge couplings we have From this equation one obtains the following explicit formula for the beta functions after taking into account that the α bare i do not depend on μ The first term in the first factor of Eq. (24) originates from the term μ 2ϵ in Eq. (23) and vanishes in four space-time dimensions. Equations (23) and (24) hold for the Yukawa couplings only for models where the Yukawa matrices are diagonal, e.g., in case only the third generation Yukawa couplings are taken into account. The generalization of Eq. (24) to incorporate tensorlike couplings, such as the Yukawa matrices and quartic couplings, is straightforward to derive. However, care has to be taken when computing derivatives of renormalization constants. Furthermore, the relations between Yukawa matrix and quartic coupling beta functions and the corresponding renormalization constants take a slightly different form than in Eq. (24), since in general, due to the tensorial nature, it is not possible to compute the inverse of the renormalization constants. For more details see Ref. [29].
Another option would be to derive the scale dependence of the eigenvalues of the Yukawa matrices and quartic couplings starting from the definition in Eq. (21).
Note that the one-loop results of Z α i only contain α i , whereas at two loops all other couplings are present except for the quartic couplings. The renormalization constants of the Yukawa matrices contain all couplings except the quartic couplings already at one-loop order, while the quartic couplings enter at two loops. Therefore, it is necessary to compute the renormalization constants and beta functions of the quartic couplings only at one-loop order. 2 For our calculation we use the automated setup developed for the calculation of the SM gauge beta functions to three loops [1,2]. For convenience we repeat the flowchart that illustrates the interaction of the various program packages in Fig. 1.
In a first step we implement the unbroken version of the general 2HDM discussed in Sec. II in the package FeynRules [31], which generates a model file for FeynArts [32]. The program FeynArtsToQ2E [33] works on the model file and translates it into input files for QGRAF [34] and q2e [35][36][37]. QGRAF is used for the generation of the amplitudes that are translated by q2e and exp [35][36][37] to FORM [38] code. The latter is processed by MINCER [39] and/or MATAD [40], which compute the Feynman integrals and output the ϵ expansion of the result. For the first part of the calculation up to the generation of the input files for QGRAF and q2e no parallelization is necessary. The individual steps take at most a few minutes. However, the parallelization of the horizontal part of the flowchart (cf. Fig. 1) is essential since for some of the Green's functions we have to deal with several hundred thousands of diagrams. Once QGRAF has produced the output file all following steps can be applied in parallel to blocks of diagrams that typically contain 1000 Feynman amplitudes.
We perform the calculation in Lorenz gauge using general gauge parameters for each gauge group. It is an important cross-check that they drop out in the expressions for the renormalization constants (and beta functions) of the gauge and Yukawa couplings.
The described setup is used to compute various MS renormalization constants for fields and vertices. They are required for the construction of the renormalization constants for the gauge, Yukawa, and quartic couplings.
For the SM and 2HDM with Z 2 symmetry, one can perform the calculation in a basis where all the Yukawa matrices are diagonal and the elements of the CKM matrix are present only in the vertices containing charged currents. In such a basis, the Lagrangian parameters are physical parameters and the number of free parameters is reduced to the number of independent degrees of freedom.
In Table II we list all Green's functions, which we have considered in the course of the calculations performed in this paper, and the number of generated Feynman amplitudes up to three loops. We used the following notation for the fields: B and W i denote the gauge bosons, c x refers to the ghost fields, and Φ 0 i and Φ AE i (i ¼ 1, 2) are the neutral and charged components of the scalar doublets. In Fig. 2 we show typical Feynman diagrams contributing to individual Green's functions.
Because of Ward identities there are various choices for each gauge coupling to obtain Z α i : where we have used , Z c gcg g , Z gg , and Z c gcg . For the Yukawa matrices we are restricted to vertices that involve components of the scalar doublets as well as left-and right-handed fermion fields. The explicit definition of the Yukawa matrix renormalization constants will be postponed to the next section.
At the end of this section a comment concerning γ 5 is in order. For the computation of some of Green's functions an odd number of γ 5 matrices is present in the traces. We have checked that it is sufficient to follow the prescription provided in Ref. [2] in the context of the SM. This means that a formal replacement of expressions like is applied, whereε μνρσ is antisymmetric in all indices. In practice the product of two such objects occurs, where all indices are contracted, which we replace bỹ The square brackets denote complete antisymmetrization. This leads to the correct result in the limit d → 4. We have checked explicitly that the ambiguity of OðϵÞ in Eq. (26) is multiplied by at most simple poles in ϵ and therefore does not lead to ambiguous renormalization constants and beta functions.

A. Renormalization constants
For the computation of the renormalization constants for fields, couplings, and vertices we follow the procedure described in Ref. [41]. However, since we consider general Yukawa couplings that are nondiagonal both in flavor and in doublet space, several modifications have to be applied, in particular for the calculation of the fermion and scalar wave function renormalization constants, and the renormalization constants for the Yukawa matrices. These issues are discussed in this section.
The renormalized inverse fermion propagator can be written as where P L=R ¼ ð1 ∓ γ 5 Þ=2, Σ L=R are the left-and righthanded parts of the fermion self-energy, and Z f L=R are the renormalization constants for the left-and right-handed fermion fields. Both Σ L=R and Z f L=R are matrices in flavor space where flavor indices have been suppressed. The index f ∈ fu; d; lg indicates whether the up-, down-, or lepton matrix shall be considered.
The renormalized two-point function of the scalar fields can be written as where the wave function renormalization constant Z Φ and Πðp 2 Þ are matrices in doublet space. The corresponding indices have been suppressed. In Z 2 -symmetric models Z Φ has to be diagonal, which we have checked up to the three-loop level. From Eqs. (28) and (29) the following relations can be derived: The operator K ϵ extracts the poles in ϵ. Solving these equations recursively allows one to determine the corresponding renormalization constants. Let us stress at this point that from the equations above we can compute only the Hermitian parts of the renormalization matrices Z f L ; Z f R , and Z Φ . In the SM the anti-Hermitian parts of the quark wave functions renormalizations are related to the renormalization of the CKM matrix [42][43][44]. In the next section we will also introduce anti-Hermitian contributions to the FIG. 2. Sample Feynman diagrams contributing to Green's functions that have been used for our calculation of the gauge and Yukawa coupling renormalization constants. Solid, dashed, dotted, curly, and wavy lines denote fermions, scalar bosons, ghosts, gluons, and electroweak gauge bosons, respectively. renormalization matrices defined above. However, in our case, they should not be identified with the renormalization of any physical quantity.
Let us in the next step use this information to obtain formulas that allow one to compute the renormalization constants for the scalar-fermion vertices and the Yukawa couplings. The Yukawa vertex renormalization constants for a fermion of type f can be extracted from where for convenience the flavor (α, β, γ) and the doublet (a, b) indices are shown explicitly. The sums over β and γ run over all down-(up-)type fermions in case α is a down-(up-)type fermion. Note that in Eq. (31) Γðp; 0Þ b;βγ is the vertex function where one of the external momenta is set to zero and the external fields are Φ b and fermions with flavor β and γ. Furthermore, the Yukawa coupling in the tree-level contribution of Γðp; 0Þ b;βγ is not renormalized. Once P β P 2 b¼1 Z ffΦ ab;αβ Y f b;βα 0 is obtained from Eq. (31) the Yukawa matrix renormalization constants (ΔY f a ) can be computed from This equation has to be solved iteratively for ðY f a þ ΔY f a Þ.

B. Invariants in the quark sector
The renormalization constants introduced in Eq. (30) are used to derive the corresponding anomalous dimensions. Note, however, that the anomalous dimensions might contain poles in ϵ. This is not surprising, since in the case of general Yukawa matrices we do not take into account the invariance of the theory under unitary rotations such as those given in Eq. (3). In other words, from the 72 parameters 4 of the Yukawa sector of the general 2HDM 30 can be eliminated using flavor transformations [45]. 5 In contrast to the case of Lagrangian parameters, the beta functions of the flavor invariants [cf. Eq. (18)] are finite, because the flavor symmetry relations are by construction taken into account in such quantities.
We also want to remark that the gauge couplings are trivially invariant under unitary flavor transformations and the corresponding beta functions do not suffer from uncanceled singularities. On the other hand, in analogy to the Yukawa matrices, we expect that in the case of the quartic couplings of the scalar potential only certain combinations of them have finite beta functions.
In the SM and in Z 2 -symmetric 2HDMs there are 11 flavor invariants in the quark sector as has been discussed in Sec. II. From them the six quark masses, three CKM mixing angles, and the cosine and the sign of the CPviolating phase can be derived. Their behavior under renormalization group evolution has been studied up to two loops, for example, in Ref. [25]. We have extended the analysis and checked by an explicit calculation that all 11 invariants in Eq. (18) of the SM and Z 2 -symmetric 2HDMs have finite anomalous dimensions at three loops. From the three-loop anomalous dimensions of the 11 invariants mentioned above one can derive the beta functions for the physical couplings and the CKM mixing angles to the same order.
As already mentioned in Sec. II, we have not classified all flavor invariants in the general 2HDM. However, an explicit calculation for the invariants where we sum over a and b and the trace is taken over the fermionic indices, shows that all poles cancel.

C. Invariants for a simplified model
In this subsection we consider a simplified version of the general 2HDM. Explicitly, we study the case where the Yukawa interactions for the first and second generations are neglected. As a consequence, the Yukawa matrices reduce to complex numbers, parametrizing the Yukawa couplings for the t and b quarks and only scalar flavor symmetries occur. Following Ref. [45], one observes that from the eight parameters 6 in the Yukawa sector ðn 2 − 1Þ ¼ n¼2 3 can be rotated away. We also notice that the up-and down-type Yukawa couplings transform as vectors under unitary rotations of the scalar fields; see Eq. (4) where U Q and U u;d are replaced by the identity matrix for this simplified model. We thus rotate the scalar fields with the following matrix: 4 Four complex 3 × 3 matrices with 18 parameters each. 5 The flavor symmetry group ½Uð3Þ 3 ⊗ Uð2Þ is broken by the Yukawa sector to Uð1Þ, leading to 30 broken generators. 6 We have y f i ∈ C with f ¼ u, d and i ¼ 1, 2.
Under this transformation the scalar fields change to with I t ¼ δ ij y t i y tÃ j . δ ab and ε ab denote the Kronecker delta and Levi-Cività tensor, respectively, and the sum over the repeated indices i, j ¼ 1, 2 is assumed. Furthermore, the Yukawa couplings transform as Taking into account the tensorial properties of δ ab and ε ab and the transformations of the Yukawa couplings under unitary rotations of the scalar fields, one can easily prove that both the rotated fields and couplings are actually flavor invariants. In other words, in the new basis the Lagrangian parameters are expressed through flavor invariants and are therefore directly related to physical quantities. An explicit calculation shows that the anomalous dimensions of the new fields Φ 0 a with a ¼ 1, 2 and the beta functions of the new couplings y t0 1 ; y b0 1 , and y b0 2 are finite through three loops. This is not the case for the original basis, where the Yukawa sector contains too many parameters. The new basis makes use of the flavor symmetries and gets rid of one of the up-type Yukawa couplings; the other one is rendered real. It is also important to notice that the relation y t0 2 ¼ 0 is stable under renormalization. To verify this statement, we checked through three loops that the beta functions of the three Yukawa couplings obtained after the rotation to the new basis can be expressed only in terms of couplings present in this basis. At this stage, also the scalar quartic couplings have to be transformed according to Eq. (8). This shows that the set of couplings in the new basis is complete. Even for this simplified model the explicit three-loop results are quite lengthy. Thus, in Sec. VI we only present results for β y t0 1 .

D. Poles in anomalous dimensions
In this subsection we describe a practical method, which allows us to use the beta functions and anomalous dimensions for a general 2HDM, although they develop poles in the first place. A transformation to physical observables, which, as mentioned above, becomes quite involved, is not necessary. We follow Ref. [6], where this issue has been discussed for the case of the SM. It is argued that the poles can be eliminated by choosing the square roots of the renormalization constants to be non-Hermitian. We define ffiffiffi ffi where Z is any of the renormalization constants introduced in Eq. (30). The subscripts H and U in Eq. (37) denote the Hermitian and unitary parts.
To obtain Y u þ ΔY u one has to invert Eq. (32). This can be done either by choosing Hermitian square roots or by allowing for additional unitary factors, which leads to the following relation: where Y u þ ΔY u is calculated with Hermitian square roots. This equation resembles the transformation in Eq. (4), and therefore a choice of the unitary part of the square root of the Z factors can be interpreted as a certain choice of basis. ffiffiffi ffi Z p H is fixed by the poles of the corresponding two point functions in Eq. (30) and can be used to obtain the Hermitian part of the corresponding anomalous dimensions, i.e., the combination γ þ γ † . For the left-and righthanded fermion fields and the scalar fields considered in Eq. (30) we observe that γ þ γ † is finite, whereas the individual terms are not.
Note, however, that ffiffiffi ffi Z p U is an arbitrary unitary matrix that can be chosen such that γ is finite. This choice is not unique, and it is possible to also influence the finite parts of the anomalous dimensions (and in general the beta functions) this way. We will postpone the discussion of this apparent ambiguity and its physical significance to Sec. IV E and concentrate in the following on the discussion of the left-handed quark fields.
At one-loop order there is no possibility to construct ffiffiffiffiffiffi Z Q L q U ≠ I, and therefore ffiffiffiffiffiffi Z Q L q is purely Hermitian. At two-loop order there is one unitary combination of Yukawa matrices where a and b are arbitrary constants. A nonzero value for b enters into the finite parts of the left-handed quark field anomalous dimension and therefore into the beta functions for Y u and Y d , contributing to the mentioned ambiguity (see below Sec. IV E). One can choose a nonzero value for a to cancel possible ϵ poles in the non-Hermitian part of the two-loop anomalous dimension. However, such poles cannot appear as can be seen by the following arguments: For the anomalous dimensions we schematically write where the second equation simplifies to in the case Hermitian square roots are chosen. Thus, the anomalous dimension is Hermitian (and finite) if the commutator μ d dμ is satisfied. At two-loop order only the contribution where both terms in Eq. (42) receive one-loop contributions involving two Yukawa matrices could possibly lead to a nonvanishing commutator of the form

E. Ambiguities in the Yukawa matrix beta function
The possibility to introduce ffiffiffi ffi Z p U ≠ I introduces an ambiguity in the definition of the renormalization constants, anomalous dimensions, and beta functions. Nevertheless, let us stress that this statement only holds for the unphysical parameters, e.g., the Yukawa matrices. Once we construct flavor invariants, the unitary roots cancel and their anomalous dimensions are finite and unambiguous. Consequently, the anomalous dimensions of the physical quantities derived from them are finite and unambiguous, as expected.
We verified the cancellation of the unitary roots and consequently the poles in the beta functions for all 11 invariants of the quark sector of Z 2 -symmetric 2HDMs. Furthermore, we checked the cancellation in the general 2HDM for the invariants introduced in Sec. II as well as further invariants entering the gauge coupling beta functions, which we will present in the next section [see Eqs. (48) and (49)].
In addition, we also performed numerical checks by computing the Yukawa matrices in Z 2 -symmetric 2HDMs at a low scale and run them up to 10 16

V. RESULTS FOR THE GAUGE COUPLING BETA FUNCTIONS
In this section we present the analytical results for the gauge coupling beta functions of a general 2HDM. We notice that both the Yukawa matrices and the self-couplings occur in the gauge beta functions only through flavor invariants, 7 as expected.
We present the results keeping the full information contained in the Yukawa matrices and arrive at the following expressions for the beta functions: and In the above equations n G denotes the number of fermion generations and n D the number of scalar doublets. We sum over the indices i, j, k, l of the quartic couplings from 1 to n D . The matrix M l is defined by in analogy to Eq. (17). Other combinations of Yukawa matrices are given by as well as TrðY l i Y l † j ÞTrðY l j Y l † i Þ; T which are defined in analogy to Eqs. (19) and (20). We rescaled the Yukawa matrices in the above results such that M f ¼ M f =ð4πÞ,T We have performed a number of cross-checks on the correctness of our result. Among them is the independence on the three gauge parameters. Furthermore, we can easily take the SM limit by setting n D ¼ 1, Y f 2 ¼ 0, and λ ij;kl ¼ λ and find agreement with Refs. [1][2][3]. We also agree with the findings of Ref. [19] where results for a general theory based on a simple gauge group are presented. 8 A comment on the validity of our results for n D ≥ 3 is in order. At three-loop order, all diagrams containing at least one internal gauge boson or a closed fermion loop can only receive contributions from up to two different scalar doublets. However, diagrams containing two quartic couplings can get contributions from more than two doublets. Therefore, all contributions to the three-loop beta functions are also valid for n D ≥ 3 apart from those containing two quartic couplings.

VI. RESULTS FOR THE YUKAWA COUPLING BETA FUNCTIONS
As discussed in Sec. IV the Yukawa matrix beta functions themselves are ambiguous and one should either work in a proper basis or only consider invariants of the Yukawa sector. In general the expressions are quite lengthy at three loops. Thus, we restrict ourselves to the beta function in the simplified model discussed in Sec. IV C. For simplicity we drop the primes introduced in Eq. (36) and write y t ≡ y t 1 since β y t 2 ¼ 0. We obtain In this context one has to take into account the comments presented at the end of Sec. IV of Ref. [2].