Complete two-loop QCD contributions to the lightest Higgs-boson mass in the MSSM with complex parameters

Higher-order corrections to the MSSM Higgs-boson masses are desirable for accurate predictions currently testable at the LHC. By comparing the prediction with the measured value of the discovered Higgs signal, viable parameter regions can be inferred. For an improved theory accuracy, we compute all two-loop corrections involving the strong coupling for the Higgs-boson mass spectrum of the MSSM with complex parameters. Apart from the dependence on the strong coupling, these contributions depend on the weak coupling and Yukawa couplings, leading to terms of Oααs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}{\left( \alpha \alpha _s\right) }$$\end{document} and Oαq1αq2αs\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}{\left( \sqrt{\alpha _{q_1}}\sqrt{\alpha _{q_2}}\alpha _s\right) }$$\end{document}, (q1,2=t,b,c,s,u,d\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{1,2}=t,b,c,s,u,d$$\end{document}). The full dependence on the external momentum and all relevant mass scales is taken into account. The calculation is performed in the Feynman-diagrammatic approach which is flexible in the choice of the employed renormalization scheme. For the phenomenological results presented here, a renormalization scheme consistent with higher-order corrections included in the code FeynHiggs is adopted. For the evaluation of the results, a total of 513 two-loop two-point integrals with up to five different mass scales are computed fully numerically using the program SecDec. A comparison with existing results in the limit of real parameters and/or vanishing external momentum is carried out, and the impact on the lightest Higgs-boson mass is discussed, including the dependence on complex phases. The new results will be included in the public code FeynHiggs.


Introduction
Since the discovery of a signal in the Higgs-boson searches at the LHC [1,2] with a mass around 125 GeV, it is a prime goal to reveal the detailed nature of the new particle. While with the present experimental and theoretical uncertainties the measured properties of the detected particle are compatible with the expectations for the Higgs boson of the Standard Model (SM) [3,4], other interpretations corresponding to very different underlying physics are also in agreement with the data. A crucial question in this context is in particular whether the observed particle is part of an extended Higgs sector that would be associated with a more general theoretical framework beyond the SM.
Within the theoretically well motivated Minimal Supersymmetric extension of the SM (MSSM), the observed particle can be interpreted as a light state within a richer spectrum of scalar particles. 1 The Higgs-boson sector of the MSSM consists of two complex scalar doublets leading to five physical Higgs bosons and three (would-be) Goldstone bosons. At the tree-level, the physical states are given by the neutral C P-even bosons h, H and the C P-odd state A, together with the charged H ± bosons. The Higgs sector at lowest order can be parametrized in terms of the A-boson mass m A and the ratio of the vacuum expectation values of the scalar doublets, tan β = v 2 v 1 . The MSSM with complex parameters (cMSSM) is of particular interest since it provides new sources of C P-violation in addition to the C Pviolating phase of the SM. Thereby the Higgs sector is C Pconserving at the tree level, but potentially large loop contributions involving complex parameters from other supersymmetric (SUSY) sectors can lead to an admixture of the C Peven states h, H, and the C P-odd A resulting in the mass eigenstates h 1 , h 2 , h 3 [6][7][8][9][10]. In this case m A is no longer a useful input parameter; instead the mass of the charged Higgs boson m H ± is used. Besides the input parameter m A or m H ± all other Higgs-boson masses are predicted quantities in the MSSM. The Higgs-boson masses and mixings in the neutral sector are strongly affected by loop contributions. Especially for the experimentally measured Higgs boson at about 125 GeV a sufficiently high accuracy of the theoretical computation is essential for drawing reliable conclusions on the viability of the investigated region of parameter space.
Restricting to the case of real parameters, the momentumdependent O(α t α s ) corrections [53][54][55] and the contributions for the case where all Yukawa couplings except the one of the top quark are neglected [54] are known. While the phases of the complex parameters affect the predictions for the Higgs-boson masses, production cross sections [77] and decays [78][79][80], they also induce C P-violating effects that are constrained by other experiments. These concern in particular the electric dipole moments [81][82][83][84][85][86][87][88]. For the usual convention where the phase of the mass of the 1 Within the MSSM it is usually assumed that the observed particle is associated with the lightest neutral Higgs boson of the model; see Ref. [5] for a recent update on the viability of the interpretation in terms of the next-to-lightest neutral Higgs boson of the MSSM. electroweakinos, φ M 2 , is set to zero without loss of generality, the phase of the parameter μ is constrained to be very close to zero or π . The other important phases of the gluino mass, φ M 3 , and the trilinear soft-breaking parameters of the stops, φ A t , and sbottoms, φ A b , are much less constrained. In particular, the bounds on the phases of the trilinear soft-breaking parameters are significantly weaker for the third generation than for the second and first generation.
In this article the full two-loop QCD corrections to the Higgs-boson masses are presented for the general case of the MSSM with complex parameters. They contain all previously computed results for the MSSM with real or complex parameters. The contributions are comprised of the O(αα s ) terms, involving the electroweak gauge coupling α, and the O √ α q 1 √ α q 2 α s terms, involving the Yukawa couplings α q 1 , α q 2 , where q 1,2 = t, b, c, s, u, d. Terms with mixed up-and down-type Yukawa couplings only appear in conjunction with m H ± as input parameter. Mixed contributions of O √ α √ α q 1 α s involving one gauge coupling and one Yukawa coupling do not appear in the final result. The results obtained here for the MSSM can furthermore be used as an approximation for higher-order contributions within the NMSSM, as discussed in Refs. [89][90][91]. The computation carried out below makes use of previously developed tools [53,92,93]. The momentum-dependent two-loop integrals appearing in the two-loop QCD corrections are evaluated with an adapted version of SecDec 2 [94][95][96]. For the numerical analysis the new contributions are combined with the full one-loop result [10] and the leading O α 2 t terms [73,74] in the Feynman-diagrammatic approach for complex parameters, available through the public program FeynHiggs [10][11][12]78,97]. In deriving the new contributions the renormalization scheme of Ref.
[10] at the oneloop level has been adopted and applied to the case of the O(αα s ) contributions. This ensures that the obtained analytical results for the renormalized two-loop self-energies can consistently be incorporated into FeynHiggs. In the results presented in this paper no resummation of higherorder logarithmic contributions as obtained in Refs. [56][57][58][59][60][61][62] has been included. The combination of resummed higherorder logarithmic contributions with the results obtained in the present paper will be addressed in future work. In the numerical analysis below, we show results for the masses of the three neutral Higgs bosons of the MSSM with complex parameters and their phase dependence, with a particular focus on those results which are phenomenologically most relevant. The paper is organized as follows: Sect. 2 provides the theoretical framework for the calculation and renormalization of the Feynman diagrams that is used to arrive at expressions for the dressed propagators of the Higgs sector up to the two-loop level. The calculation of the unrenormalized self-energies and the construction of the two-loop counterterms are described in Sect. 3. Details on the numerical evaluation of the momentum-dependent two-loop integrals are given in Sect. 4, whereas the impact of the new contributions on the Higgs-boson masses is discussed in Sect. 5. The conclusions are given in Sect. 6.

Tree-level relations for masses and mixing
The two scalar SU (2)-doublets are conventionally expressed in terms of their components as follows, with the relative phase ξ . The Higgs potential can be written as a polynomial in the field components, where terms of third and fourth power in the fields have been omitted, and the relations φ − have been used. Explicit expressions for the tadpole coefficients T and for the mass matrices M can be found in Ref. [10]. They are parametrized by the phase ξ , the real SUSY-breaking quantities m 2 1,2 =m 2 1,2 + |μ| 2 , and the complex SUSY-breaking quantity m 2 12 . With the help of a Peccei-Quinn transformation [98] the parameter m 2 12 can be redefined such that its phase vanishes [99], leaving only the phase ξ as a potential source of C P-violation at tree level. The requirement of minimizing V H at the vacuum expectation values v 1 and v 2 is equivalent to the requirement of vanishing tadpoles of the physical fields, which in turn implies the condition ξ = 0 at tree level. As a consequence, the Higgs sector of the MSSM is C P-conserving at lowest order. This implies in Eq. (2.2) that M φχ is equal to zero, and φ 1,2 do not mix with χ 1,2 at tree-level.
The remaining (2 × 2)-matrices M φ , M χ , M φ ± can be transformed into the mass eigenstate basis with the help of orthogonal matrices D(x), using the abbreviations s x ≡ sin x, c x ≡ cos x, (2. 3) The Higgs potential in this basis can be expressed as follows, The Goldstone bosons are massless, 3 m G ± = m G = 0. The masses m H ± , m A , m h , m H fulfil the relations including the vector-boson masses M W and M Z . Given the relation in Eq. (2.9a), both m A and m H ± can be chosen as input parameter.
At lowest order, the irreducible two-point vertex functions of the neutral Higgs sector (2.10) are diagonal, and the entries of the mass matrices in Eq. (2.6) provide the poles of the diagonal lowest-order propagators (2.11)

Masses and mixing beyond lowest order
Going beyond leading order, the irreducible two-point functions are dressed by adding the matrixˆ h H AG of the renormalized diagonal and non-diagonal self-energies for the h, H, A, G fields up to the considered order, (2.12) yielding the full renormalized two-point vertex function The latter generally contains a mixing of all fields with equal quantum numbers. The dressed propagators are obtained by inverting the matrixˆ h H AG ( p 2 ). Truncating the perturbative expansion at the two-loop level, the momentum-dependent corrections to the neutral Higgs-boson mass matrices in Eq. (2.12) are given by (2.14) For the MSSM with complex parameters, the one-loop self-energies are completely known [10], and the leading two-loop O(α t α s ) and O α 2 t + α t α b + α 2 b contributions have been obtained in the approximation of zero external momentum [72][73][74][75]. In the case of the MSSM with real parameters also the momentum-dependent corrections of O(α t α s ) are known [53][54][55]. The remaining QCD contributions at the two-loop level are completed within this paper. These contributions comprise terms of the O(α x α s ), where α x is either the gauge coupling α or the Yukawa coupling α q with q = {u, d, s, c, b, t}. We neglect CKM mixing for those contributions.
In order to obtain the physical Higgs-boson masses from the dressed propagators at the considered order, it is sufficient to explicitly derive the entries of the h H A p 2 with M (2) h H A being the corresponding (3 × 3)-submatrix of Eq. (2.14). The impact of the self-energies on the mixing and couplings of the various Higgs bosons to other (MS)SM particles can be obtained with the same formalism as described in Refs. [10,100].

Calculation of the renormalized two-loop self-energies
The renormalized two-loop self-energies can be written aŝ (2) where (2) h H A denotes the unrenormalized self-energies corresponding to the sum of genuine two-loop diagrams and one-loop diagrams with counterterm insertions. The symbol δ (2) M Z h H A comprises all two-loop counterterms resulting from parameter and field renormalization.
The contributing types of Feynman diagrams for the calculation of the full two-loop QCD corrections entering Eq.  is equal to zero in that case. The diagrammatic calculation has been performed with the help of FeynArts [101,102] in generating the Feynman diagrams, and TwoCalc [92] and Reduze [103] for the two-loop trace evaluation and tensor reduction. The one-loop renormalization constants have been computed with the help of FormCalc [104].

Two-loop counterterms
In order to obtain the renormalized self-energies in Eq. (3.1), counterterms have to be introduced up to second order in the loop expansion, for the tadpoles and for the mass matrices of Eq. (2.4) The two-loop counterterms of O(αα s ) have the same structure as the corresponding one-loop counterterms. They are listed here for completeness and to fix our notation.
In order to ensure the correct form of the counterterms for the mass matrices, the rotation angles β n and β c from Eq. (2.3) have to be distinguished from β in Eq. (2.7). Whereas no renormalization is needed for α, β n and β c , a counterterm associated with β of the form β → β + δβ is required, in accordance with the renormalization of tan β, (3.4) In the resulting expressions for the counterterm matrices, the identification β c = β n = β can be made, see Ref. [10] for details of the analogous treatment at the one-loop order (note that a different convention for the counterterm of t β is used in Ref. [10]). A complete list of the two-loop counterterms is given in the Appendix of Ref. [74]. In addition to the parameter renormalization described previously, the field-renormalization constants δ (1) Z H i and δ (2) Z H i are introduced at the one-loop and two-loop order (restricting the latter to the contributions of O(αα s )) for each of the scalar doublets of Eq. (2.1) through the transformation The field-renormalization constants in the mass-eigenstate basis of Eq. (2.3) are obtained via The matrices Z i j can be expanded as The required one-loop expressions for the entries in the δ (1) Z i j -matrices are given in Ref. [10]; the corresponding set of two-loop expressions is given in Ref. [74]. The genuine two-loop counterterms δ (2) M Z h H A of Eq. (3.1) can now be summarized as where the required two-loop mass counterterms read The entries of δ (2) M h H A that are not listed here are determined by symmetry. When replacing δ (2) → δ they are formally equal to the one-loop counterterms listed in Eqs. (53) of Ref. [10] (up to the different convention for the counterterm of t β used there). The two-loop renormalization constants of Eqs. (3.2)-(3.9) are fixed by extending the renormalization scheme of Ref. [10] from the one-loop to the two-loop order: • The tadpole counterterms δ (2) T i are fixed by requiring that the minimum of the Higgs potential is not shifted, which means that the tadpole coefficients have to vanish at each order. At the two-loop level, the condition reads where the T (2) i denote the unrenormalized one-point functions at two-loop order, see Fig. 2 for the contributing two-loop diagrams. The aforementioned relation for the mixing angles β n = β c = β is a consequence of the tadpole conditions T i = 0 at lowest order.
• The charged Higgs-boson mass m H ± is the only independent mass parameter of the Higgs sector and is used as an input quantity. Accordingly, the corresponding mass counterterm is fixed by an independent renormalization condition, chosen as on-shell, given by The renormalized charged-Higgs self-energy at the twoloop level can be expressed in terms of the unrenormalized charged self-energy and its respective countertermŝ leading to the mass counterterm when applying the on-shell condition. The contributing Feynman diagrams are shown in Fig. 3. As we neglect flavor mixing q, o,q andõ always belong to the same generation. As a consequence, the vertices with four squarks in topologies 12, 14 and 15 are only non-zero when all adjacent fields are of the same generation. • The field-renormalization constants of the Higgs mass eigenstates in Eq. (3.6) are combinations of the basic doublet-field renormalization constants δ (2) Z H 1 and  (2) Z H 2 , which are fixed by the UV-divergent parts of the derivatives of the corresponding self-energies, (3.14) • Also t β is renormalized by a purely UV-divergent counterterm, which was shown to be a convenient choice [105] (see also Refs. [106,107]). Alternative process-dependent definitions for the renormalization of t β can be found in Ref. [108]. For the class of two-loop corrections of O(αα s ) the counterterm can be written as • When neglecting momentum-dependent contributions and taking the gaugeless limit, the purely UV-divergent two-loop counterterms δ (2) Z H 1 , δ (2) Z H 2 and δ (2) t β cancel each other and are therefore not required for renormalization, compare Ref.
[75]. If one of these two limitations is dropped, δ (2) Z H 1 , δ (2) Z H 2 and δ (2) t β are necessary in order to obtain a UV-finite result. In the corrections discussed in this article these counterterms have to be taken into account as none of these approximations is used. It should also be noted that the chosen renormalization conditions for δ (2) Z H 2 and δ (2) t β are not equal to pure DR conditions, since the top-mass counterterm δ (1) m t which enters in δ (2) Z H 2 is fixed by an on-shell condition. The resulting differences between the two schemes have been discussed in [54,55]. • Renormalization of the D terms in the Higgs-squark couplings which are induced by the gauge coupling g 2 , as well as the relation between the charged and C Podd Higgs masses require counterterms for the Z -and W -boson masses, δ (2) M 2 Z and δ (2) M 2 W , respectively. We treat M W and M Z as independent input parameters and fix their renormalization constants by the on-shell conditions leading to Here (2) Z ,W denote the transverse parts of the two-loop self-energies of Z and W , respectively. Most of the Feynman diagrams contributing to the twoloop self-energies (2) Z ,W differ from the Higgs selfenergies depicted in Figs. 1 and 3 only by the external fields. Pictorially, they can be obtained by replacing the neutral external Higgs fields by the Z , and the charged Higgs field by the W boson field. All additional topologies are depicted in Fig. 4.

Sub-loop renormalization
Apart from the genuine two-loop diagrams, the lowest-order QCD contributions to the self-energies and tadpoles involve one-loop diagrams with insertions of one-loop counterterms. This subrenormalization concerns masses and mixing of the colored particles.
The required one-loop counterterms for subrenormalization arise from the quark q and scalar quarkq sectors. The squark mass matrices in the q L ,q R bases are given in lowest order by with Q q and T 3 q denoting charge and isospin of q ∈ {u, c, t, d, s, b}. For the sake of convenience we suppress repeating the indices of the first and second generation in the following since renormalization is analogous to the third generation.
The squark mass eigenvalues can be obtained from unitary transformations, Since A q and μ are complex parameters, the unitary matrices Uq can be described by the mixing angle θq and an additional phase ϕq . The independent parameters which enter the two-loop calculation through the quark-squark sector are: the quark masses m q , the soft SUSY-breaking parameters mQ i and mq R , and the complex trilinear couplings A q = |A q | e iφ Aq . These parameters have to be renormalized at the one-loop level, thus defining transformations Mq → Mq + δ (1) Mq for the mass matrices in Eq. (3.18). The other free parameter μ, which is related to the Higgsino sector, enters the selfenergies as well. However, the renormalization of μ does not receive one-loop corrections of O(α s ) and is therefore not part of the contributions considered in this calculation. The individual renormalization conditions for the colored sector are formulated as follows: • Renormalization of the top quark mass is carried out in the on-shell scheme, i.e.
where the quark self-energy is given in terms of its Lorentz decomposition with the left-, right-handed projectors ω −,+ = 1 2 (1 ∓ γ 5 ). The bottom mass is renormalized in the DR scheme (see Refs. [41,42,109]) at the scale m os t . The counterterm can be obtained by using the expression in analogy to the counterterm for the top quark mass in Eq. • In order to fix the renormalization constants of the stop sector, we employ the relation Thus we derive (3.24c) The counterterm δ (1)  are fixed by onshell conditions for the top-squarks, is fixed by the renormalization condition (see Ref. [72]) which involves the non-diagonal squark self-energy shown in Fig. 5 with incomingt 2 and outgoingt 1 . • Between the gauge and mass eigenstates of the bottom squarks we employ an analogous relation to Eq. (3.23). We derive As indicated by the subscript, we choose to renormalize A b in the DR scheme, which has been shown to be convenient for reasons of numerical stability [41,42,109]. The scale of A b is chosen to be m os t . As a consequence of SU (2) invariance the counterterm δ (1) is not independent, but a derived quantity which is fixed by the renormalization of the top-stop sector in Eq. (3.24a), since (3.28) Inserting Eq. (3.27a) and solving for δ (1)  are fixed analogously as for the stops: Therefore in our scheme only mb 2 is renormalized onshell, while the counterterm δ (1) is a derived quantity according to Eq. (3.29).

Resummation of tan β-enhanced terms
The Yukawa coupling of the bottom quark h b receives radiative corrections proportional to tan β. Those tan βenhanced contributions can be resummed as described in Refs. [80,[110][111][112][113][114][115]. The resummed contributions b are UV finite and generally yield complex numerical results. For the numerical evaluation in Sect. 5, we use the version for b at the one-loop order which is implemented in FeynHiggs and outlined in the following. The largest tan β-enhanced contributions can be absorbed by using an effective bottomquark mass, which is related to the DR-renormalized bottom quark mass in the MSSM as follows, (3.32) The expression b at the one-loop order contains contributions from gluinos, charginos and neutralinos (ordered in decreasing numerical importance) and reads The couplings α s and α are running parameters and are evaluated at the scale Q = mb 1 mb 2 . The symbol C depicts the CKM matrix, and u g ,ũ g are the gth generation up-type quarks and squarks, whereas B 0 (0, m 1 , m 2 ) and B 1 (0, m 1 , m 2 ) are one-loop functions. 5 As mentioned above, we otherwise neglect CKM mixing in the two-loop contributions that we evaluate. The renormalization scale μ r from the loop integrals drops out. The coefficients c L,R and n L,R are given by 5 These one-loop functions are given by In order to obtain a full conversion of the bottom mass between the on-shell scheme and the DR scheme in Eq. At the one-loop order they read The parameters entering in b and δ b are computed in the limit of large tan β. The chargino and neutralino masses and mixing matrices are then obtained as where we use (3.37) Thereby the matrices Uχ and Vχ yield a singular value decomposition for X, and the matrix Nχ yields Takagi's factorization [116] on Y. The sbottom masses in this limit are computed from the matrix given in Eq. (3.18) at A b = 0. Since the bottom mass itself also enters that matrix, the final solution for m b,eff is found iteratively.
By using Eq. (3.31) for the bottom mass in the one-loop contributions to the Higgs masses, the leading higher-order corrections to the Higgs masses from the bottom-sbottom sector are generated. The contributions of the bottomsbottom sector to the two-loop self-energies presented in this article add further subleading shifts. It should be noted that the expression given in Eq. (3.31), which employs the DR scheme in the MSSM, is chosen such that no double counting of the terms contained in m b,eff occurs at the twoloop level.

Numerical evaluation of the self-energies
The renormalized two-loop self-energies are expressed in terms of two-loop two-point multi-scale integrals with nonzero external momenta. With the help of TwoCalc [92] and Reduze [103] all integrals can be reduced to the four irreducible scalar two-loop topologies depicted in Fig. 6, and products of analytically well-known one-loop one-and twopoint functions.
The scalar two-loop integrals are defined as The irreducible two-loop integrals of Fig. 6 may depend on up to five different internal mass scales taken from the following set, in addition to a non-zero external momentum, taking the values of p 2 = M 2 h 1 , M 2 h 2 , M 2 h 3 when entering the unrenormalized self-energies, or p 2 = m 2 H ± , m 2 W , m 2 Z when entering the self-energies through two-loop renormalization constants. Recently, a lot of progress has been made towards Fig. 6 Irreducible two-loop topologies resulting from tensor reduction, calculated numerically with the program SecDec. Some of the internal lines may also be massless describing and evaluating integrals of this class analytically [117][118][119][120][121][122][123][124][125][126]. However, to the best of our knowledge, an implementation of the analytical results for all topologies in Fig. 6 is not publicly available. We have therefore calculated these integrals numerically using the program SecDec [94][95][96].
For the evaluation, the resulting new contributions to the neutral Higgs-boson self-energies have been added to FeynHiggs via its interface to external programs, see section 2.4 of Ref. [53] for details. We have extended the existing interface to the program SecDec in FeynHiggs to deal with the 177 mass configurations of which 88 are computed at four different kinematic points, 72 at two and the rest at one kinematic point. The parameters entering the integrals are evaluated by FeynHiggs and passed on to SecDec. It should be noted that the heavy growth of mass configurations with respect to non-electroweak corrections is due to an increase in the number of mass scales involved in the renormalized self-energies.
We constructed two independent integration setups to allow for consistency checks of the numerical result. The two-point one-loop topologies entering the self-energies up to O(ε) are known analytically. The bulk of their implementation was previously tested in Ref. [53] and compared with the authors of Ref. [54]. Additional mass configurations were newly implemented and checked against SecDec. The increase in two-loop mass configurations by more than a factor five with respect to the previous setup in Ref. [53] calls for a higher precision of the integrals to avoid numerical instabilities due to cancellations. With the integral reduction, unphysical thresholds can be introduced which cancel in the sum of all contributing diagrams. Numerically, due to round-off errors, the cancellation might however not always be exact, leading to numerical instabilities. The latter are cured by introducing a small imaginary part to the denominators of the coefficients arising from the integral reduction. We have verified that the numerical dependence of the selfenergies on this technical regularization parameter is negligible.
The fact that we take a non-zero value of the bottom quark mass into account leads to a large hierarchy among the different mass scales. Numerical convergence at the desired accuracy is therefore difficult to accomplish. On the other hand, we have analyzed the influence of the quark masses of the first and second generation on the two-loop integrals in the self-energies. For the second generation and tan β 1 a negative shift in the Higgs-boson mass correction of only about 20 MeV can be observed when neglecting the light quark masses. The effect is even smaller for the quark masses of the first generation. The terms which involve the light quark masses in couplings are negligible, too. It is due to this reason that we will assume the first and second generation quarks to be massless throughout the rest of our numerical analysis. The numerical impact of the gauge contributions of the light quarks will be discussed below.
In order to achieve a relative precision of at least 10 −7 for each integral, we use the deterministic integrator Cuhre included in the Cuba library [127,128] but have optimized the integration parameters for each integral topology and mass configuration individually.
As a further crosscheck of our computation, we have compared the O(αα s ) contribution by the top-stop and bottom-sbottom particles to the Z -boson self-energy which is required for renormalization of the Higgs sector. Since [54] uses massless bottom quarks, we have reevaluated our result for the Z self-energy in the limit m b = 0. In order to avoid a dependence on the renormalization scheme of the quarksquark sectors, the Z -boson self-energy has been evaluated in the DR scheme by both groups for this comparison. Overall we have found a very good agreement with discrepancies at the level of 0.3 GeV 2 .
We find an overall uncertainty of the self-energies entering the light Higgs-boson mass of maximally 0.2% by adding all uncertainties on the numerical evaluation of the twoloop integrals in quadrature. Given the resulting size of our newly computed corrections analyzed in the next section, the absolute uncertainty on the light Higgs boson mass is maximally 0.4 MeV, which is well below the shift coming from neglecting light quark masses from the first and second generation.
The total of 513 integrals have been computed numerically on the fly before passing the resulting two-loop self-energies back to FeynHiggs, where they are added to the corresponding matrix elements just before the determination of the propagator poles.

Numerical results for the Higgs mass spectrum
In the following we analyze the numerical impact of the newly computed corrections. We start with a comparison with earlier results in the literature and then discuss our results in three different scenarios: an m mod h -like scenario (based on Ref. [129]), a scenario with a particularly large value of tan β where contributions from the bottom and sbottom sector are enhanced, and a low-m H scenario (inspired by Refs. [5,129]). For better readability of the results, we define three different Higgs-boson masses resulting from different higher-order contributions As mentioned earlier, the quark masses and Yukawa couplings of the first and second family are neglected. Thus, the first and second generation contributes only at O(αα s ) by Dterm contributions of the sfermions. We focus our numerical discussion on the fixed-order result up to the two-loop level, i.e. no combination with resummed higher-order logarithmic contributions as discussed in Refs. [59][60][61] is employed.
Using the definitions of Eq. (5.1), we assign The size of the effects of our newly computed contributions is contained in M h i , since all the previously known terms are subtracted. So far, the two-loop terms of O(α b α s ) were only known in the MSSM with real parameters and m A as input parameter. M h i shows our new contributions without these terms, if m A is chosen as input parameter.
Below we will discuss our results for non-zero phases of complex parameters. We investigate in particular the variation of the phases φ M 3 , φ A t and φ A b , which are much less constrained by experimental bounds on EDMs than the phases of μ, M 1 (in the usual convention where the parameter M 2 is chosen to be real) and the phases of the trilinear couplings of the first and second generation. As discussed e.g. in Ref. [130], scenarios with relatively large phase values are possible. In order to demonstrate the possible impact of the phase variations on the Higgs spectrum, below we display the phase dependences over the whole range (−π, π].

Comparison with earlier results
In a first step, in Table 1 we show a comparison of the results for the light Higgs-boson mass including our new contributions with the results of Ref. [54], where in the MSSM with real parameters the corrections of O(α t α s p 2 ) and the full corrections of O(αα s ) have been evaluated, and with the results up to O(α t α s p 2 ) in the MSSM with real parameters from Ref. [53]. The comparison is carried out for the benchmark scenarios m max h , m mod+ h , m mod− h defined in Ref. [129] and for a modified light-stop scenario used in Ref. [131]. We find overall good agreement with the results of Ref. [54].   [129] using M SUSY = 2, 3 TeV and m A = 500 GeV, tan β = 20. The results are compared with those provided by the authors of Ref. [54] for two different wave-function renormalization schemes

This publication
Ref. [54] δ (2) Z Ref. [53] H i δ (2) Z Ref. [53] H i δ (2)  yield a small downward shift. The numerical differences between the results for the contributions of O(α t α s p 2 ) from Refs. [53,54], which amount up to 0.3 GeV for the examples considered here, result from different renormalization scheme choices of δ (2) Z H i , see the discussion in Refs. [53][54][55]. Those differences in the renormalization schemes also affect the comparison between our results for M new h and the results for M old h + O(α t α s p 2 ) + O(αα s ) from Ref. [54] in Table 1.
The differences in the renormalization schemes and the dependence on the parameters in the stop sector are further investigated in Table 2 Table 2. Two versions of the results from Ref. [54] are shown, one using the renormalization scheme adopted in Ref. [54] with δ (2) H i , and the other using the renormalization scheme of Ref. [53], which we have adopted in the present work, with δ (2) H i . 6 It can be seen in Table 2 that there is very good agreement, at the level of about 10 MeV, between our results and the results from Ref. [54] using the renormalization scheme of Ref. [53]. The different choices of renormalization schemes in the result of Ref. [54] amount to mass shifts of up to 150 MeV for the displayed examples. The 6 We are very grateful to S. de Vita for providing us with those results. Compared to the original m mod h scenario we choose larger bilinear soft-breaking parameters for the sfermions, and also larger absolute values for μ (see below) and M 2 . Thereby mQ 3 is slightly different from m {t,b} R in order to avoid numerical instabilities by degeneracies. However, the general feature of this scenario is kept: it allows for a wide range of X t = A * t − μ/ tan β to be in agreement with experimental bounds. With our choice of parameters, A t and A b are not expected to be affected by constraints from chargeand color-breaking minima [132][133][134][135][136][137][138][139]. As A τ has negligible impact on the Higgs mass prediction, we set it to zero.
First, the dependence of the lightest Higgs-boson mass M h 1 on tan β is analyzed for different values of the μ parameter. Setting all phases of the parameters that can be complex to zero, our result can be compared to previous ones in the MSSM with real parameters where the corrections evaluated in the present paper were not included. In the considered scenario, it is possible to choose either m A or m H ± = m 2 A + M 2 W as an input parameter which is renormalized on-shell accordingly. The chosen input mass for Fig. 7 is m A . A comparison of the predicted mass from  Values for tan β above the depicted range and large negative μ lead to a further rapid decrease of M h 1 , eventually yielding a tachyonic Higgs boson. This is due to the large bottom Yukawa coupling with resummed tan β-enhanced terms which can become nonperturbative in that region of the parameter space. The rise of the red dotted curve at large tan β reflects that this decrease happens for larger values of tan β once our new corrections are taken into account. In Fig. 8 the charged Higgs mass m H ± is used as an input parameter. The latter implies the occurrence of terms of O √ α q √ α o α s and corresponds to the renormalization scheme compatible with both the MSSM with real and complex parameters. On the left-hand side of Fig. 8, the blue (μ = 500 GeV) and red (μ = −1500 GeV) lines show the prediction for the lightest Higgs mass with (solid) and without (dashed) our new contributions. In addition, the solid and dashed curves of Fig. 7 are indicated again as grey (μ = 500 GeV) and black (μ = −1500 GeV) lines. In this way, the influence of the two different renormalization schemes on the Higgs-mass prediction can be seen. While the blue and grey lines lie on top of each other over the whole range of tan β, deviations of up to 1.5 GeV can be observed between the red and black curves in the region of large tan β. Since the slope of the red curves for large tan β is smaller than for the black curves, the renormalization scheme with m H ± as input parameter is better suited for this particular region in parameter space. On the right-hand side of Fig. 8

the mass shifts M h 1 and
M h 1 resulting from our new contributions are depicted. The color coding is the same as described before. The size of the shifts is almost invariant under the exchange of m A and m H ± as input parameter, since only small differences between the two renormalization schemes can be noticed. We note that setting μ = 1500 GeV and using m H ± as input, the same qualitative behavior as for the lower positive μ value can be observed, with the new contributions being of the same size as for μ = −1500 GeV in the low and intermediate tan β region. Furthermore, the size of the mass shift M h 1 in Figs. 7 and 8 shows a similar tendency with respect to the chosen sfermion masses as depicted in Table 2, i.e. larger scales increase the size of the new corrections. However, for stop-and sbottom masses larger than ≈ 2 TeV logarithmic contributions of higher order also become important. Then, a resummation of these logarithms should be taken into account for an accurate Higgs-mass prediction. The gluino mass can have a sizable impact due to its appearance in the threshold correction of O(α t α s ).   Fig. 11. Parameters are as described in Eq. (5.4) Accordingly, the difference between the two curves is given by the pure gauge contributions of O(αα s ) from the first and second generation. They are rather small, amounting to about 30 MeV.
The variation of M h 1 with the gluino-mass parameter M 3 = |M 3 | exp i φ M 3 is shown in Fig. 10. Close to |M 3 | ≈ 1.9 TeV, thresholds of the gluino-fermion-sfermion system can be observed, which are introduced by one-loop integrals entering via the subloop-renormalization and resummation of the bottom Yukawa coupling. The effect of varying the absolute value of the gluino-mass parameter |M 3 | on M h 1 is strongest for φ M 3 = 0 and successively weakened as φ M 3 In Figs. 11, 12 and 13 the dependence on the three phases φ M 3 , φ A t and φ A b is displayed, respectively. The impact of the new (solid) corrections in comparison with the ones implemented so far in FeynHiggs (dashed) are shown for the lightest Higgs-boson mass on the left-hand side of each figure, while the differences M h 1 are shown on the right-hand side. Comparing to the MSSM with real parameters, where the phases are equal to zero or π , sizable differences for the prediction of the lightest Higgs-boson mass are visible. Concerning the total variation of M h 1 including all now available corrections, the impact of the phases φ A t and φ M 3 is seen to be rather large with effects that can exceed 2 GeV, while varying the phase φ A b yields only rather small shifts of ≈ 0.2 GeV.
The prediction for M h 1 as function of φ M 3 shown in Fig. 11 is symmetric with respect to the sign of φ M 3 . The variation of M h 1 with φ M 3 is shown on the right-hand side of Fig. 11. The pronounced dependence on the absolute value of |M 3 | seen in Fig. 10 can be observed again. The variation of φ M 3 changes M h 1 by up to 250 MeV for an |M 3 | value around the gluino-fermion-sfermion threshold, while for |M 3 | = 2.5 GeV M h 1 is shifted only by up to 70 MeV.
The phase dependence of M h 1 on φ A t and φ A b is shown on the right-hand side of Figs. 12 and 13, respectively. The variation of M h 1 with φ A t and φ A b is seen to be rather small. It reaches up to 150 MeV for the phase φ A t and up to 50 MeV for φ A b . It should be noted that the results for φ M 3 = ± π 2 lie on top of each other in Fig. 13. While the variation with φ A b is rather small for any non-zero φ M 3 , the variation with φ A t is minimal for φ M 3 = 0 and maximal for φ M 3 = π .
Using different values of φ A b (and keeping φ M 3 fixed) has only a small effect on the variation of M h 1 with φ A t . The corresponding plot is therefore not shown here.

Scenario 3: low M H
In the low-M H scenario the observed SM-like Higgs boson with a mass of about 125 GeV can be identified with the next- l ∈ e, μ. (5.5) Compared to the original scenario in [5] we had to choose a smaller value of μ in order to avoid a tachyonic lightest Higgs boson for a charged Higgs mass m H ± ≈ 160 GeV. Our value for tan β is chosen such that the scenario is valid according to Fig. 26 of [5]. In Fig. 14 the three neutral Higgs-boson masses are depicted, varying the charged Higgs-boson mass m H ± which is used as an input parameter. The light green band illustrates the mass range of 125 ± 3 GeV; it should be interpreted as a rough indication of the mass range which is theoretically in agreement with the detected Higgs boson. Up to m H ± 188 GeV the heavier Higgs h 2 could be associated with the discovered Higgs-like particle; however, as can be seen in the low-M alt+ H scenario in Fig. 26 of [5], our choice of μ and tan β is already excluded for a charged Higgs mass m H ± = 185 GeV. Yet, scenarios with values of m H ± closer to or below m t are still allowed. In this region the new corrections presented here have a negligible impact on M h 2 , but lead to a downward shift of about 1 GeV for both M h 1 and M h 3 .
As shown in Fig. 15, using a non-zero value of the gluino phase of φ M 3 = π/2 or φ M 3 = π shifts all three neutral Higgs masses to larger values as compared to the case φ M 3 = 0. For better comparison, the results of Fig. 14 are underlaid in grey. The numerical impact of the new contributions presented here rises with increasing φ M 3 . For φ M 3 = π all neutral Higgs masses can receive large corrections of up to 5 GeV.

Conclusions
We have computed the full two-loop QCD corrections to the lightest Higgs-boson mass in the MSSM with complex parameters. Compared to previous works, this primarily involves going beyond the gaugeless limit, and including a finite bottom-quark mass; furthermore the momentum dependence of loop integrals is taken into account. On the technical side, this involves the computation of 177 different mass topologies evaluated at different kinematical configurations, amounting to a total of 513 two-point two-loop integrals with up to five mass scales. These integrals have been computed numerically with the program SecDec. The results of Fig. 14 with φ M3 = 0 are depicted in grey for reference. Parameters are as described in Eq. (5.5) In the first part of our numerical analysis, we have compared our results with earlier result in the literature taking the appropriate limit of real parameters and/or vanishing external momentum of our results. We have found very good agreement with the existing results in the appropriate limit if the same renormalization scheme is employed. The contributions evaluated in this paper yield a shift in the lightest Higgsboson mass at the level of 1 GeV, where the impact has been seen to be more pronounced for an increasing mass scale of the stops.
We have furthermore investigated the dependence of the new corrections on tan β choosing different values of the μ-parameter as well as different renormalization schemes. For a large negative μ the corrections are generally larger and amount to around 0.9 GeV in M h 1 . The corrections are largest for 10 < tan β < 30, decrease by 3% for lower values and by about 20% beyond tan β = 30.
We find non-vanishing mixed up-and down-type Yukawa corrections in the charged Higgs-boson self-energy correction entering the mass predictions for the neutral Higgs bosons as renormalization constant if the charged Higgs mass m H ± instead of the neutral CP-odd mass m A is chosen as an input parameter. We have compared the mass prediction for the lightest Higgs boson in both schemes and have found good agreement in general. However, using the charged Higgs mass as an input parameter yields better numerical stability at large tan β and large negative μ.
The Yukawa contributions scale according to their Yukawa couplings, leading to much smaller contributions from the first and second generation quarks and squarks. The pure gauge terms of O(αα s ) in the limit of massless quarks are found to be of similar small size, below 20 MeV for one generation.
Analyzing the dependence on the gluino mass, we have found maximal shifts of ≈ 900 MeV in M h 1 . The correc-tions show a sensitive dependence on the gluino-fermionsfermion threshold, which enters via the counterterms of our renormalization scheme, and the gluino phase. For the μparameter a mass shift of the lightest Higgs by ≈ 850 MeV is found over large regions of parameter space.
Concerning the impact of the three phases φ M 3 , φ A t and φ A b , we find significant effects in our new corrections from varying the gluino phase and the pase of A t . For φ M 3 the phase dependence becomes particularly pronounced in the threshold region of the gluino-fermion-sfermion system, as mentioned above.
Besides scenarios where the lightest neutral Higgs boson in the spectrum of the MSSM is the SM-like state that can be identified with the detected Higgs signal, we have also analyzed the impact of the newly computed contributions on the Higgs-mass predictions for the three neutral Higgs bosons within the low-M H scenario for different values of the gluino phase φ M 3 . We have found mass corrections of ≈ 1 GeV for φ M 3 = 0 and up to ≈ 5 GeV for φ M 3 = π in this case.
Accordingly, we have found that the subleading two-loop contributions that we have evaluated in this paper yield a shift in the prediction for the mass of the light SM-like Higgs boson of the MSSM of up to the level of 1 GeV. The size of the correction sensitively depends on the mass scales of the stops and sbottoms, on the absolute value and phase of the gluino mass parameter, as well as on the absolute value and phase of the trilinear coupling in the stop sector (and to a lesser extent on the trilinear coupling in the sbottom sector). While these findings of course have an impact on the remaining theoretical uncertainties from unknown higher-order corrections, we do not attempt to provide an improved estimate of the remaining uncertainties here. Such an improved estimate should be based on a combination of the fixed-order result considered here with a resummation of higher-order loga-rithmic contributions. We leave such an analysis to future work.
It should be noted in this context that our results for the corrections of O(αα s ) beyond the gaugeless limit cannot be used directly to infer the possible size of the corresponding contributions of O α 2 to the Higgs-boson spectrum, which are unknown up to now. This is due to the fact that the requirement of a strong coupling in the corrections that we have evaluated significantly constrains the structure of the contributing Feynman diagrams, while additional classes of contributions will have to be taken into account for a full calculation of the corrections of O α 2 .
The new contributions evaluated in this paper will be made publicly available in the program FeynHiggs.