On three-loop RGE for the Higgs sector of 2HDM

We discuss renormalization group equations (RGE) for the parameters of the Higgs sector in general Two-Higgs-Doublet Model (2HDM). We present the three-loop results but consider only contributions due to self-couplings of the Higgs doublets. We study the structure of RGE and express beta-functions in terms of reparametrization invariants with respect to higgs-basis rotations. The Cayley-Hamilton theorem is utilized to reduce both the number of independent tensor structures in matrix RGE and the number of invariants to a minimal set. As a by-product of our calculation we discovered that two-loop RGE of the scalar sector in general QFT with multiple higgses were not properly implemented in a number of public packages. The latter give rise to a wrong result when mixing in the scalar sector is allowed.


Introduction
The Standard Model (SM) was established in mid-1970s. Its success is incredible: even after almost half a century, no significant deviations from the SM predictions were found. Given a minimal set of parameters, the SM provides a very precise description of different phenomena in Modern Particle Physics. To confront its predictions with ongoing and future experiments, one is forced to take various radiative corrections into account and, in many cases, carry out certain kind of re-summation. A convenient tool to deal with high-order terms in perturbative expansion is dimensional regularization [1] accompanied by modified minimal (MS) subtractions of infinities. The latter appear in loop integrals and manifest itself as poles in ǫ = (4 − d)/2.
In MS-renormalization scheme the model parameters depend on auxiliary scale µ and their numerical values at different scales are related by differential renormalization group equations (RGE). While boundary conditions should be extracted from experiment, the RG functions (beta functions and anomalous dimensions) can be calculated order-by-order in perturbation theory. Solution of RGE allows one to improve the precision of finite-order predictions by re-summing certain logarithmic corrections into redefinition of the model parameters.
For the parameters of the SM Lagrangain three-loop RG functions are known from recent computations: the gauge coupling are considered in refs. [2][3][4], the results for Yukawa couplings can be found in refs. [5][6][7], and refs. [8,9] are devoted to the SM Higgs-potential parameters. There are also partial four-loop results available in literature (see., refs [10][11][12][13]. Recently, five-loop RG functions in pure QCD have been calculated [14][15][16]. Among other things, all the results were immediately applied to state-of-the-art studies [17,18] of the vacuum-stability problem in the SM. In spite of the above-mentioned success of the SM, there are well-known issues (related to dark matter, fine-tuning, etc.) that prevent us from treating the SM as the most fundamental theory of particle interactions (see, e.g., ref. [19]). Among different possibilities to go beyond the SM (BSM) one can consider an extension with an additional Higgs doublet -the so-called Two-Higgs-Doublet Model (2HDM) (for review see refs. [20,21]). The model predicts new scalar states in the spectrum -two neutral H, A and one charged H ± higgs bosons. Being (linear combinations of) components of the SU(2) doublets, their interactions with vector fields are fixed by postulated gauge symmetry, but there is a freedom in selfinteractions and fermion Yukawa couplings.
Recently, three-loop beta-functions for the gauge and Yukawa sector of general (Type-III) 2HDM were found in ref. [22]. In this paper, we continue the study of the RG functions in 2HDM and calculate certain three-loop contributions to the beta functions of the higgs self couplings and anomalous dimensions of the higgs mass parameters entering general Higgs potential.
We restrict ourselves to the corrections due to the scalar self-interactions only and for the moment we neglect both gauge and Yukawa couplings of the Higgs fields. Nevertheless, we consider different parameterization of the scalar sector. In addition, we compute the scale dependence of reparametrization invariants (see, e.g, ref. [23]), which are constructed from the Higgs potential parameters, but contrary to the latter, do not depend on the choice of Higgs basis. For convenience, all the RG functions considered in this paper 1 are available as ancillary files of the arXiv version of the paper.
It is worth mentioning that we have tried to compare the two-loop beta-functions obtained by direct calculations with the RG-functions extracted from the well-known results for a general renormalizable QFT model [24,25]. We have found that application of the available general result to the case with many scalar fields requires some care. We discovered that, e.g., current versions of SARAH [26] and PyR@TE [27], when applied to the case of Type-III 2HDM, give rise to a wrong result at two loops (see Section 6 for details).
The paper is organized as follows. In section 2 we introduce the 2HDM Higgs potential and discuss various parametrizations of the Higgs sector. In section 3 the structure of the RG functions in one particular parametrization, which involve scalars Λ 00 , M 0 , 3vectors Λ, M and a symmetric 3 × 3 tensor Λ, is elaborated. section 4 is devoted to the description of the renormalization procedure for the above-mentioned quantities. The corresponding three-loop RG functions can be found in section 5. We discuss subtleties in the interpretation of the well-known two-loop expressions [24,25] and present our results for self-couplings λ i and masses m 2 ij in section 6. Our conclusions can be found in section 7. In a series of appendices we provide some details on the reparametrization-invariant counting via Hilbert Series (A) and present useful identities for Λ (B). In addition, the RG-functions of the reparametrization invariants are given in appendix C.

The scalar potential of 2HDM
The most general renormalizable Higgs potential can be written in the following form with Φ 1,2 being SU (2) doublets. The self-couplings λ 1−4 and the mass parameters m 2 11 , m 2 22 are real, while λ 5−7 , and m 2 12 can be complex. Not all of these fourteen (real) parameters are physical due to the freedom in redefinition of Higgs basis by a unitary rotation where a, b = 1, 2 enumerate the doublets. It is easy to see that the overall phase of U does not impact the change for the couplings and masses so in what follows we restrict ourselves to U ∈ SU(2). The three parameters of SU(2) rotation can be used to get rid of three out of 14 parameters of the potential and, thus, we are left only with 11 independent quantities. There is an alternative notation [28] which can be used as an intermediate step to rewrite the self-couplings (see refs. [20,21]) in terms of a scalar Λ 00 , a vector Λ and a symmetric matrix Λ, where σ µ ≡ (1, σ), µ, ν = 0, 1, 2, 3, i, j = 1, 2, 3 and the euclidean metric is used both for four-and three-dimensional indices. The same trick can be used for the mass term: We can also decompose the tensor in terms of a singlet r 0 and a vector r. By means of eqs.(2.4),(2.6), and (2.8) one can rewrite the potential (2.3) as Under a Higgs-basis change Φ a → U ab Φ b , U ab ∈ SU(2), Λ 00 and M 0 transform as singlets, while Λ and M transform as triplets under the corresponding SO(3) rotation The symmetric 3 × 3 matrix Λ ≡ {Λ ij } can be decomposed 2 into a singlet trΛ and a five-pletΛ ij ≡ Λ ij − 1 3 trΛδ ij .

The structure of RG functions
The parametrization of the quartic couplings in terms of Λ 00 , Λ and Λ turns out to be very convenient for calculation of RGE in the scalar sector. The main advantage of the approach is that we need to deal with at most two indices instead of four. In addition, the Cayley-Hamilton theorem, which states that the square 3 × 3 matrix Λ satisfies its own characteristic equation can be used to get rid of high powers Λ n (n ≥ 3) appearing at the intermediate steps of calculation. Due to eq. (3.1) we can enumerate possible structures that can appear in beta-functions for the components of Λ µν (t = ln µ 2 , h = (16π 2 ) −1 ):

2)
Since Λ 00 is an invariant (singlet w.r.t higgs-basis transformations), only reparametrization invariants can enter β Λ 00 . Given Λ 00 , Λ, and Λ one can introduce the following independent invariants 3 I i,j (c.f. [23]): 3c) The first index in I i,j corresponds to the order (or degree) of the invariant, i.e., the total power of the Λ µν components entering I i,j . There is also an invariant of order six which is related (up to a sign) to those presented in eq. (3.3). A convenient tool to enumerate the invariants is the so-called Hilbert Series (see appendix A and references therein). The beta-function β Λ can be cast into the general form 4 with a i being polynomials in invariants (3.3). The beta-function β Λ looks like where Λ ⊗ Λ ≡ Λ i Λ j , etc. and b i are again expressed in terms of invariants. By means of Cayley-Hamilton theorem one can also show that other symmetric tensors constructed from Λ and Λ are not independent (see appendix B for details).
Since in the MS scheme counter-terms are polynomial in momenta and masses [1], it is clear from dimensional analysis that RGE for M µ can only involve first powers of the latter. As a consequence, scalars (invariants) and vectors involving high powers of M and M 0 will not contribute to the mass anomalous dimensions, which we define here as Indeed, the RG equation for M 0 should be a linear combination of the following reparametrization invariants with coefficients being polynomials in the invariants built from Λ µν only. The anomalous dimension γ M must be a linear combination of the vectors where I M denotes one of the invariants from eq. (3.8). The results of direct evaluation of Feynman graphs (see, e.g., eqs. (5.4), and (5.5)) confirm this structure.

Renormalization procedure
In order to find RGE for dimensionless couplings we generate diagrams (self-energies Γ b a , and four-point functions Γ bd ac ) with external Φ a ,Φ † b , etc., but rewrite the quartic vertex in terms of Λ µν by means of eq. (2.4). We heavily rely on automatic index-summation algorithms of FORM [29,30] to deal with indices of different dimensions in diagrams generated by DIANA [31]. To extract the corrections to Λ 00 , Λ, and Λ from the considered Green functions we apply projectors, which imply summation over external higgs indices. The form of the projectors can be deduced from eqs. (2.4) and (2.8).
Let us briefly discuss counter-terms originating from the Lagrangian in the notation of eq. (2.9). It is convenient to consider additive renormalization of the parameters, i.e., The bare bilinear Φ a combinations (r µ ) bare are given by where z 0 and z come from the decomposition of the hermitian field renormalization constant The counter-term Lagrangian is obtained by expressing the bare fields and parameters in terms of renormalized ones by means of the above-mentioned equations. The expressions , and z = O(Λ 2 µν ) are determined order by order in perturbation theory.
The renormalization constants in the MS scheme are extracted from divergent terms of the corresponding loop integrals. Due to this, we made use of the well-known infrared rearrangement (IRR) tricks [32], which allow us to modify the infrared structure 5 of the considered integrals and convert them to fully massive bubbles. A modern version 6 of the MATAD [33] package written in FORM was used to compute the vacuum integrals.
Given δΛ µν and δM µ we find beta-functions and mass anomalous dimensions via differentiation of the bare parameters (4.2) w.r.t. the scale t = ln µ 2 : Both δΛ µν and δM µ involve higher poles in ǫ. However, the corresponding contribution to the RG functions is canceled due to the so-called pole equations [34]. As a consequence, the finiteness of (4.6) in the limit ǫ → 0 serves as a cross-check of the correctness of our final results. It turns out that one needs to utilize various relations (see appendix B) to prove that the pole equations are satisfied. In order to find δM µ in the MS scheme, it is sufficient to treat the mass term as a perturbation to the massless theory. The corresponding (bare) Lagrangian can be rewritten as where renormalized operators [Φ † σ µ Φ] are related to the bare bilinears via From (4.7) one can see that the renormalization constants Z 00 , Z i ,Z i and Z ij also enter mass-parameter renormalization (4.2) Due to this, we extract the mass-parameter counter-terms not from massive self-energies with external Φ † and Φ, but from divergences of auxiliary three-point functions with an additional (Φ † σ µ Φ)-operator insertion at zero momentum. The latter are computed by means of the above-mentioned IRR trick.

Three-loop RGE for Λ µν and M µ
The procedure discussed in the previous section was used to find RG functions for the Higgs potential parameters (2.9). The one-, two-and three-loop results for Λ µν are given by the expressions: The mass-parameter anomalous dimensions can be cast in the following form  By means of simple algebra one can easily convert these results to the expressions for RG functions β λ i , γ m 2 ij of the initial Higgs potential (2.1) (see next section) or to the beta-functions for reparametrization invariants (see appendix C) .

RGE for λ i and m 2 ij
We define the RG functions of the parameters λ i and m 2 ij from the potential given in eq. (2.1) as: Having in mind that 3) one can obtain the three-loop results for β λ i (6.1). For brevity we present here only oneand two-loop contributions 7 :
The anomalous dimensions of the mass parameters (2.1) can be cast into Again, γ m 2 22 can be deduced from γ m 2 11 if we substitute λ 1 ↔ λ 2 , λ 6 ↔ λ 7 and m 2 11 ↔ m 2 22 . It is worth noting that the one-loop results presented here coincide with that given in ref. [20,35]. In the case of real λ 1−5 , m 2 12 and vanishing λ 6 = λ 7 our two-loop expressions coincide with those presented in ref. [36] (2HDM with soft Z 2 breaking).
Before going to conclusions, let us briefly comment on the peculiar fact about wellknown two-loop result [24,25] for beta-function of quartic coupling λ abcd (see eq. (4.2) of ref. [24] or eq. (37) from ref. [25]) in a general renormalizable QFT model: Here, γ abcd is the anomalous dimension of the operator φ a φ b φ c φ d , while γ S i is said to be "the anomalous dimension of the scalar field i". The subtlety we encountered is the interpretation of the last term in eq. (6.20) that comes from the renormalization of the fields φ a , which suppose to carry a (gauge) index a. In general, the anomalous-dimension matrix of scalar fields is non-diagonal γ S ab and the mixing due to dimension-4 operators has to be taken into account. For example, for non-zero λ 6 , λ 7 the two-loop propagator corrections in 2HDM give rise to the mixing between Φ 1 and Φ 2 (see also ref. [37]) . It is interesting to note the in the Yukawa beta-functions possible mixing is taken into account in the general formula (see, e.g,. eq. (32) of ref. [25]). However, in ref. [25] γ S i seems to be interpreted as eigenvalues 8 of γ S ab , which in our opinion leads to an incorrect result, when scalar indices a, b, c, d are not related to a gauge group 9 , i.e., the expression (6.20) should be rewritten as Obviously, if γ ab is diagonal, eq. (6.21) leads to eq. (6.20). The same problem can appear in calculation of mass-parameter RGE (c.f, eq.(90) of ref. [25]). One explicit argument for eq. (6.21) comes from the following fact. We tried to use eq. (6.20) to compute two-loop β λ i . From the definition of Λ µν (2.5) one can easily find β Λµν , but again written in terms of λ i . It turns out that we were not able to convert the obtained result for β Λµν to the form (see eqs. (3.5) and (3.6)), which involve only invariants and certain vectors/tensors constructed from Λ µν . On the contrary, the expression (6.21) gives rise to the same results (5.1), (5.2) and (5.3) that we derived via explicit computation of Feynman graphs.

Conclusions
We considered the three-loop RGE for the scalar sector of general 2HDM in the limit of vanishing gauge and Yukawa couplings. In spite of the fact that the obtained result is obviously incomplete and can not be used in phenomenological analysis of the model, it can be treated as a necessary step towards full three-loop beta-functions.
A convenient parameterization of the scalar sector was utilized to deal with combinatorics of (tensor) self-couplings and to restrict the general form of the beta-functions in the considered limit. Scalar coefficients turn out to be polynomials in a finite set of invariants. The latter give rise to a basis-independent parameterization of the Higgs sector.
The approach can be easily extended to the case of gauge interactions since the corresponding couplings are bilinear in Higgs fields. However, the extension to the case of Yukawa interactions is not straightforward. Due to this, we will study these peculiarities elsewhere. It is also worth noting that one can make use of the public code FMFT [38] to generalize the obtained result to the four-loop case.
An important by-product of the paper is that one should be careful, when interpreting the two-loop beta-functions of scalar self-couplings presented in refs. [24,25]. We discovered that public software [26,27], which can be used to generate RGEs for any renormalizable model, do produce incorrect results when scalar sector with multiple higgs doublets is considered and mixing between the scalar states by dimension-4 operators is allowed. We expect that physical analyses based on RGE obtained by means of these codes may be inaccurate.
Note added: Slightly after the results of the present work were made public, a paper with comprehensive study [39] of the two-loop RGEs in general renormalizable QFT appeared on the arXiv. The authors of the reference confirmed our findings and extended the analysis to the case of cubic scalar couplings and scalar mass terms 10 . The correct expressions involving non-diagonal anomalous dimensions of the scalar fields will be incorporated in forthcoming versions of SARAH and PyR@TE in the near future.

Acknowledgements
The author would like to thank Veronika Rutberg, Andrey Pikelner, Andrey Onischenko, Mikhail Kalmykov, and Dmitri Kazakov for fruitful discussions. The work is supported in parts by the RFBR grant No. 17-02-00872-a and the Heisenberg-Landau Programme.

A Hilbert Series and the number of Reparametrization Invariants
A convenient way to enumerate quantities invariant under some group is the so-called Hilbert Series (see, e.g., refs. [40][41][42][43] for various applications in BSM Physics, Flavour Physics and Supersymmetric gauge theories). The series are defined as where c n gives the number of invariants of degree n, and c 0 = 1. The expression of H(t) can be constructed from pure group-theoretical considerations. We develop a simple MATHEMATICA code based on the LieART package [44] to derive the series for the invariants that can be constructed by contracting different representations of SU(2) group.
In the case of invariants, built from quartic couplings only, we obtain .
The number of factors in the denominator (p = 7) gives us the number of independent parameters corresponding to quartic interactions. They are encoded in the independent invariants. The order of the invariants corresponds to the power α in a denominator factor (1 − t α ). There are two of them of the order one (3.3a), two of the order two (3.3b), two of the order three (3.3c), and one of the order four (3.3d). The products of these independent invariants give rise to higher-degree invariants, the number of which can be obtained via expansion of the denominator in t. The presence of numerator different from one tells us about an additional degree-six invariant, which, however, can be eliminated when raised to the second power. One can also consider multi-graded Hilbert series for the self-couplings and introduce a separate variable t i for each irreducible representation used to construct an invariant quantity 11 : .

(A.3)
This time not only the degree of invariants can be read off the series, but also their composition in terms of different representations. Again, the denominator corresponds to the basic invariants I i,j (3.3). From the numerator of (A.3) one deduces that an invariant 11 The expression (A.2) is recovered from (A.3) in the limit ti → t.
of order six (3.4) should be built from three instances of Λ and three instances of (the traceless part of) Λ. One can also prove the following relation expressing the square of the degree-six invariant in terms of the invariants from the denominator (3.3) From the denominator one can immediately deduce the number of physical parameters (p = 11) of the scalar potential of 2HDM.

B Useful identities
Here we list useful identities that were utilized to obtain compact expressions for the betafunctions and to check the corresponding pole equations. From the Cayley-Hamilton theorem (3.1) valid for a 3 × 3 matrix Λ = A + h B with A and B also being matrices, at order h we have A 2 B + ABA + BA 2 = trB · A 2 + trA · (AB + BA) This identity can be used to prove that the degree-four tensor and the degree-five tensor can be reduced to the structures given in eq. (3.6). Moreover, the order-6 matrix is also reducible. In addition, from (3.1) the following identities that were utilized in due course of our calculations can be derived Finally, it is interesting to note the relation valid in 3d for arbitrary vectors a and b and arbitrary symmetric matrix C a × b · trC = C · ( a × b) + (C · a) × b + a × (C · b). (B.8) The relation was used to simplify the result for the mass-term RGE.

C Scale dependence of the invariants
The RG functions for reparametrization invariants are defined as i,j , t = ln µ 2 , h = 1 16π 2 . (C.1) We present here the expressions for β (l) i,j up to the two-loop order 12