1 Introduction

After the discovery of a Higgs boson [1, 2] with a mass around 125 GeV, intense studies were performed to reveal its nature. Although within the present experimental uncertainties the measured properties of this new boson are consistent with the expectations for the Higgs boson of the Standard Model (SM) [3, 4], it could be part of an extended model like the theoretically well motivated minimal supersymmetric Standard Model (MSSM). In the MSSM the observed particle could in principle be interpreted as one of the three neutral physical Higgs bosons. At the tree level, the physical states are given by the neutral CP-even hH and CP-odd A bosons, together with the charged \(H^{\pm }\) bosons, and can be parametrized in terms of the A-boson mass \(m_A\) and the ratio of the two vacuum expectation values, \(\tan \beta = \left. v_2\diagup v_1\right. \). An admixture of these CP eigenstates is introduced to the Higgs sector via loop contributions involving complex parameters from other supersymmetric (SUSY) sectors [5,6,7,8].

Loop corrections to the masses of the Higgs bosons are sizable and therefore phenomenologically very important. Accordingly, numerous calculations for higher-order corrections to the mass spectrum within the MSSM for the case where CP conservation has been assumed [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63] as well as for the general case of the MSSM with complex parameters (cMSSM) [5,6,7,8, 61,62,63,64,65,66,67,68,69,70,71,72] have already been performed. The largest one-loop contributions originate from the Yukawa sector due to the size of the top-quark Yukawa coupling \(h_t\), where \(\alpha _t=\left. h_t^2/(4\pi )\right. \). For large values of \(\tan \beta \) contributions of the order \(\alpha _b=\left. h_b^2/(4\pi )\right. \), with the bottom Yukawa coupling \(h_b\), can become sizable. At the two-loop level both types of contributions receive further potentially large corrections. The dominant contribution is given by the leading \({\mathcal {O}}{\left( \alpha _{t}\alpha _{s}\right) }\) terms [26,27,28, 69] which are known in the MSSM with complex parameters. Additional corrections involving the strong coupling \(\alpha _s\) are known in the special case of the CP-conserving MSSM [49,50,51,52,53]. Another important class of two-loop corrections are Yukawa-coupling enhanced contributions of the order \({\mathcal {O}}{\left( \alpha _t^2+\alpha _t\alpha _b+\alpha _b^2\right) }\) which are known in the CP-conserving MSSM as well [36, 39] (the corrections of \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) with on-shell parameters are only available in the approximation \(\tan \beta \rightarrow \infty \) and \(m_b\rightarrow 0\)). A computation of the leading corrections of \({\mathcal {O}}{\left( \alpha _t^2\right) }\) has been published for the general MSSM [70, 71]. In this article also the other pieces of the two-loop Yukawa terms are obtained for the general case of the MSSM with complex parameters.

The phases of complex parameters in the MSSM are constrained by limits on electric dipole moments (EDMs) [73,74,75,76,77,78], the impact of meson mixings and decays (see Ref. [79] and references therein), and Higgs-coupling measurements [4].

Following the usual convention, we choose to fix the phase of the mass of the electroweakinos, \(\phi _{M_2}\), to zero; then the phase of \(\mu \) from the superpotential, \(\phi _\mu \), needs to be close to zero or \(\pi \) in order to be compatible with the experimental constraints. The other relevant parameters are the phase of the gluino mass parameter, \(\phi _{M_3}\), and the trilinear soft-breaking parameters of the stops, \(\phi _{A_t}\), and sbottoms, \(\phi _{A_b}\). These phases, \(\phi _{M_3}\), \(\phi _{A_t}\) and \(\phi _{A_b}\), are less constrained; especially the bounds on the phases of the trilinear soft-breaking parameters are weaker for the third generation than for the second and first generation.

The calculation presented here extends the Yukawa-type contributions of \({\mathcal {O}}{\left( \alpha _t^2\right) }\) from Refs. [70, 71], and profits from previously developed tools [80]. For this reason the theoretical framework is just briefly outlined and only new aspects are explained in detail in Sect. 2. The numerical analysis in Sect. 3 is focused on the impact of the new contributions on the Higgs masses, showing substantial shifts of 2 GeV and more in certain regions of parameter space. In the limit of vanishing phases of the parameters, our results agree well with previous results in the MSSM for the case where CP conservation has been assumed [39]; differences are shown for the comparison with an interpolation for non-zero phases which so far has been used in FeynHiggs [28, 29, 68, 81, 82]. The new results will become part of the public code FeynHiggs.

2 Higgs masses at higher orders in the complex MSSM

In this section, we briefly outline the theoretical framework for Higgs-mass predictions at higher orders in the MSSM. We introduce our conventions, explain some details of our chosen renormalization scheme, and comment on the gauge-less limit and the bottom mass resummation.

2.1 Notation and conventions at the tree level

The two scalar SU(2) Higgs doublets are expressed in terms of their components in the following way:

$$\begin{aligned} {\mathcal {H}}_{1}= & {} \left( \begin{array}{*{20}c} {v_{1} + \frac{1}{\sqrt{2}}(\phi _{1} - {i}\chi _{1})}\\ {-\phi ^{-}_{1}}\\ \end{array}\right) ,\nonumber \\ {\mathcal {H}}_{2}= & {} \left( \begin{array}{*{20}c} {\phi ^{+}_{2}}\\ {v_{2} + \frac{1}{\sqrt{2}}(\phi _{2} + {i}\chi _{2})} \end{array}\right) . \end{aligned}$$
(2.1)

After rotation to mass eigenstates, the Higgs potential reads

$$\begin{aligned} V_{H}&= -T_h \, h- T_H \, H - T_A \, A - T_G\, G\nonumber \\&\quad + \frac{1}{2}\begin{pmatrix} h,&H,&A,&G \end{pmatrix} {\mathbf {M}}_{hHAG} \begin{pmatrix} h \\ H\\ A\\ G\end{pmatrix}\nonumber \\&\quad + \begin{pmatrix} H^{-},&G^{-}\end{pmatrix} {\mathbf {M}}_{H^\pm G^\pm } \begin{pmatrix} H^{+}\\ G^{+}\end{pmatrix} + \dots \ , \end{aligned}$$
(2.2)

with the tadpole coefficients \(T_{h,H,A,G}\), and the mass matrices

$$\begin{aligned} {\mathbf {M}}_{hHAG}&= \begin{pmatrix} m^2_{h} &{}\quad m^2_{hH} &{}\quad m^2_{hA} &{}\quad m^2_{hG} \\ m^2_{hH} &{}\quad m^2_{H} &{}\quad m^2_{HA} &{}\quad m^2_{HG} \\ m^2_{hA} &{}\quad m^2_{HA} &{}\quad m^2_{A} &{}\quad m^2_{AG} \\ m^2_{hG} &{}\quad m^2_{HG} &{}\quad m^2_{AG} &{}\quad m^2_{G} \end{pmatrix} , \nonumber \\ {\mathbf {M}}_{H^\pm G^\pm }&= \begin{pmatrix} m^2_{H^\pm } &{}\quad m^2_{H^-G^+} \\ m^2_{G^-H^+} &{}\quad m^2_{G^\pm } \end{pmatrix}. \end{aligned}$$
(2.3)

The matrices \({\mathbf {M}}_{hHAG}\) and \({\mathbf {M}}_{H^\pm G^\pm }\) are diagonal at the tree level after minimizing the potential. Explicit expressions for the entries are given in Ref. [68].

2.2 Gauge-less limit

The gauge-less limit in our calculation is defined by neglecting all couplings proportional to \(g_1\) or \(g_2\). As a consequence of this approximation the gauge-boson masses \(M_W\) and \(M_Z\) are equal to zero in the new two-loop contributions.

Accordingly, the Higgs-boson masses entering the two-loop calculation take on the values

$$\begin{aligned} m_h = m_G = m_{G^\pm } = 0,\quad m_H = m_A = m_{H^\pm }\,. \end{aligned}$$
(2.4)

In this limit, the tree-level mixing angles \(\alpha \in \left[ -\pi /2,0\right) \) and \(\beta \in \left[ 0,\pi /2\right) \) fulfill the relation

$$\begin{aligned} \alpha&= \beta - \frac{\pi }{2}\,. \end{aligned}$$
(2.5)

2.3 Higgs masses at the two-loop order

The Higgs mass matrix elements at the two-loop order receive contributions from self-energies, leading in general to mixing of all neutral states. In this article the full one-loop corrections are used, while the \({\mathcal {O}}{\left( \alpha _t\alpha _s\right) }\) and the new \({\mathcal {O}}{\left( \alpha _t^2+\alpha _t\alpha _b+\alpha _b^2\right) }\) terms are evaluated in the gauge-less limit and at zero external momentum. Therefore, the loop-corrected propagator \(\Delta _{h H A G}\) is given by

$$\begin{aligned}&\Delta _{h H A G}(p^2) \nonumber \\&\quad = {i}\left[ p^2\mathbf{1}- {\mathbf {M}}_{ h H A G}^{(0)} + \hat{\mathbf {\Sigma }}_{h H A G}^{(1)}(p^2) + \hat{\mathbf {\Sigma }}_{h H A G}^{(2)}(0)\right] ^{-1}\, . \end{aligned}$$
(2.6)

Therein, \( \hat{\mathbf {\Sigma }}_{h H A G} ^{(k)}\) denotes the matrix of the renormalized diagonal and non-diagonal self-energies for the hHAG fields at loop order k, and \({\mathbf {M}}_{ h H A G}^{(0)}\) denotes the diagonal tree-level mass matrix.

Mixing of the Goldstone boson (and of the longitudinal Z boson) with the other Higgs bosons yields negligible effects to the propagators of the physical Higgs bosons [83,84,85]. Therefore, in the following we will only consider the \(\left( 3\times 3\right) \) submatrix of \(\Delta _{h H A G}\) involving the physical Higgs bosons. Though, Goldstone–Higgs mixing is taken into account in subloop renormalization terms of the type \(\left( \text {one-loop}\right) ^2\) [71].

The neutral Higgs masses are derived from the real parts of the complex poles of the hHA propagator matrix, obtained as the zeros of the determinant of the renormalized two-point function,

$$\begin{aligned} {\text {det}}{\hat{\Gamma }}_{hHA}{\left( p^2\right) }&= 0 \nonumber \\ {\hat{\Gamma }}_{hHA}{\left( p^2\right) }&= {i}\left[ p^2 {\mathbf{1}} - {\mathbf {M}}_{ h H A}^{(0)} + \hat{\mathbf {\Sigma }}_{h H A}^{(1)}(p^2) + \hat{\mathbf {\Sigma }}_{h H A}^{(2)}(0)\right] . \end{aligned}$$
(2.7)

2.4 Counterterms

The renormalized two-loop self-energies can be written as

$$\begin{aligned} \hat{\mathbf {\Sigma }}_{h H A} ^{(2)} (p^2)&= \mathbf {\Sigma }_{h H A} ^{(2)} (p^2) - \delta ^{(2)} \mathbf {M}^{\mathbf {Z}}_{hHA} \, , \end{aligned}$$
(2.8)

with \(\mathbf {\Sigma }_{h H A} ^{(2)}\) denoting the unrenormalized self-energies at the two-loop order, and \(\delta ^{(2)} \mathbf {M}^{\mathbf {Z}}_{hHA}\) comprising all two-loop counterterms resulting from parameter and field renormalization. The notation follows [71], where the required expressions for \(\delta ^{(2)} \mathbf {M}^{\mathbf {Z}}_{hHA}\) can be found.

The Feynman-diagrammatic calculation of the self-energies has been performed with the help of FeynArts [86, 87] for the generation of the Feynman diagrams, and TwoCalc [88] for the two-loop tensor reduction and trace evaluation. The one-loop renormalization constants have been obtained with the help of FormCalc [89].

2.4.1 Genuine two-loop renormalization

The two-loop counterterms for the Higgs self-energies given in Ref. [71] also apply to the corrections described in the present article. However, there is an interesting difference for the cancelation of the divergence in the self-energy \(\Sigma ^{(2)}_{hH}(0)\). The corresponding counterterm reads

$$\begin{aligned} \delta ^{(2)}m_{hH}^{{\mathbf {Z}}}&= \delta ^{(2)}m_{hH}^2 + \tfrac{1}{2}m_{H^{\pm }}^2\delta ^{(2)}Z_{hH} + \dots \,, \end{aligned}$$
(2.9)

where terms with products of two one-loop counterterms have been omitted. In the gauge-less limit \(\delta ^{(2)}m_{hH}^2\) is the only counterterm which contains \(\delta ^{(2)}t_{\beta }\),

$$\begin{aligned} \delta ^{(2)}m_{hH}^2&= m_{H^{\pm }}^2\,c_{\beta }^2\,\delta ^{(2)}t_{\beta } + \dots \,. \end{aligned}$$
(2.10)

Here we define \(t_\beta \equiv \tan \beta \), \(s_\beta \equiv \sin \beta \) and \(c_\beta \equiv \cos \beta \). The two-loop field renormalization constant for the same matrix element is given by

$$\begin{aligned} \delta ^{(2)}Z_{hH}= & {} -c_{\beta }s_{\beta }\left[ \delta ^{(2)}Z_{{\mathcal {H}}_2} - \tfrac{1}{4}\left( \delta ^{(1)}Z_{{\mathcal {H}}_2}\right) ^2\right. \nonumber \\&\left. - \delta ^{(2)}Z_{{\mathcal {H}}_1} + \tfrac{1}{4}\left( \delta ^{(1)}Z_{{\mathcal {H}}_1}\right) ^2\right] \,. \end{aligned}$$
(2.11)

The two-loop counterterm for \(t_{\beta }\) in the \(\overline{\text {DR}}\) scheme and in the gauge-less limit can be expressed as

$$\begin{aligned} \delta ^{(2)}t_{\beta }= & {} \frac{t_{\beta }}{2}\left[ \left( \delta ^{(2)}Z_{{\mathcal {H}}_2} - \delta ^{(2)}Z_{{\mathcal {H}}_1}\right) \right. \nonumber \\&\left. - \tfrac{1}{4}\left( \delta ^{(1)}Z_{{\mathcal {H}}_2} - \delta ^{(1)}Z_{{\mathcal {H}}_1}\right) ^2\right. \nonumber \\&\left. -\left( \delta ^{(1)}Z_{{\mathcal {H}}_2} - \delta ^{(1)}Z_{{\mathcal {H}}_1}\right) \delta ^{(1)}Z_{{\mathcal {H}}_1}\right] \,. \end{aligned}$$
(2.12)

Combining Eqs. (2.9)–(2.12) yields

$$\begin{aligned} \delta ^{(2)}m_{hH}^{{\mathbf {Z}}} = \frac{c_{\beta }\,s_{\beta }\,m_{H^{\pm }}^2}{4}\left( \delta ^{(1)}Z_{{\mathcal {H}}_2} - \delta ^{(1)}Z_{{\mathcal {H}}_1}\right) \delta ^{(1)}Z_{{\mathcal {H}}_1} + \dots \,. \end{aligned}$$
(2.13)

The \(\overline{\text {DR}}\) field-renormalization constant \(\delta ^{(1)}Z_{{\mathcal {H}}_1}\) is a pure UV-divergent term, calculated as the derivative with respect to the external momentum \(p^2\) of the \(\phi _1\) Higgs self-energy. The only contribution in the gauge-less limit is a bottom loop, i.e. in the case of the previously calculated \({\mathcal {O}}{\left( \alpha _t^2\right) }\) corrections [71], \(\delta ^{(1)}Z_{{\mathcal {H}}_1}\) was equal to zero due to the approximation \(m_b=0\). The terms originating from two-loop field-renormalization and two-loop renormalization of \(t_{\beta }\) canceled each other exactly.

Now, for the more general case of a non-zero bottom mass, also \(\delta ^{(1)}Z_{{\mathcal {H}}_1}\) is non-zero and the cancelation is not complete anymore. The genuine two-loop parts of the field-renormalization constants, \(\delta ^{(2)}Z_{{\mathcal {H}}_1}\) and \(\delta ^{(2)}Z_{{\mathcal {H}}_2}\), drop out in the gauge-less limit at zero external momentum in Eq. (2.13) because of a cancelation of the contributions in \(\delta ^{(2)}t_\beta \) and \(\delta ^{(2)}Z_{hH}\). In principle, \(\delta ^{(2)}Z_{{\mathcal {H}}_1}\) and \(\delta ^{(2)}Z_{{\mathcal {H}}_2}\) could still appear as field-renormalization constants for the other Higgs-mass counterterms. However, also there they drop out exactly (see Eq. (2.23) in [71]):

  • for \(\delta ^{(2)}m_h^{{\mathbf {Z}}}\) since \(m_h^2 = 0\) in the gauge-less limit,

  • for \(\delta ^{(2)}m_H^{{\mathbf {Z}}}\) since \(m_H^2 = m_{H^\pm }^2\) and \(\alpha = \beta - \tfrac{\pi }{2}\) in the gauge-less limit,

  • for \(\delta ^{(2)}m_A^{{\mathbf {Z}}}\) since \(m_A^2 = m_{H^\pm }^2\) in the gauge-less limit,

  • for \(\delta ^{(2)}m_{hA}^{{\mathbf {Z}}}\) and \(\delta ^{(2)}m_{HA}^{{\mathbf {Z}}}\) since the Higgs sector is CP conserving at the tree level.

2.4.2 Resummation

Radiative corrections to the relation between the bottom-quark mass and the Yukawa coupling of the bottom quark \(h_b\) are proportional to \(t_{\beta }\). In order to resum the leading \(t_\beta \)-enhanced contributions, an effective bottom Yukawa coupling is used as described in Refs. [84, 90,91,92,93,94,95], leading to a UV finite and complex correction factor \(\Delta m_b\). Using a \(\overline{\text {DR}}\) renormalization for \(m_b\) in the MSSM, the largest contributions of this type are captured through an effective bottom-quark mass which is given by

$$\begin{aligned} m_b^{\overline{\text {DR}},\text {MSSM}}(m_t^{\text {os}}) \simeq m_{b,\text {eff}}&= \frac{m_b^{\overline{\text {DR}},\text {SM}}(m_t^{\text {os}})}{|1 + \Delta m_b |}. \end{aligned}$$
(2.14)

The symbol \(m_b^{\overline{\text {DR}},\text {SM}}(m_t^{\text {os}})\) denotes the bottom mass in the \(\overline{\text {DR}}\) renormalization scheme, taking into account SM-type QCD corrections, evaluated at the on-shell top mass.

We use the correction factor \(\Delta m_b\) at the one-loop order which is implemented in FeynHiggs. For illustrating the effects seen in our numerical analysis below, we give here the explicit form of the leading contributions:

$$\begin{aligned}&\Delta m_b = \frac{2\,\alpha _s}{3\,\pi }\,\mu ^*\,M_3^*\,t_\beta \,{\mathcal {I}}{\left( m_{{\tilde{b}}_1}^2,m_{{\tilde{b}}_2}^2,m_{{\tilde{g}}}^2\right) }\nonumber \\&\qquad \qquad +\frac{\alpha _t}{4\,\pi }\,\mu ^*\,A_t^*\,t_\beta \,{\mathcal {I}}{\left( m_{{\tilde{t}}_1}^2,m_{{\tilde{t}}_2}^2,|\mu |^2\right) }\,, \end{aligned}$$
(2.15a)
$$\begin{aligned}&{\mathcal {I}}{\left( a,b,c\right) } = -\frac{b\,a\,\log {\tfrac{b}{a}} + c\,b\,\log {\tfrac{c}{b}} + a\,c\,\log {\tfrac{a}{c}}}{\left( b - a\right) \left( c - b\right) \left( a - c\right) }. \end{aligned}$$
(2.15b)

Further subleading contributions involve terms \(\propto \alpha _b\) and \(\propto \alpha \). With this definition a part of the considered two-loop corrections of \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) to the Higgs-boson masses are absorbed into an effective bottom-quark mass. In order to avoid a double counting of contributions from the bottom–sbottom sector to the Higgs-boson self-energies, the bottom-mass is renormalized in the \(\overline{\text {DR}}\) scheme as specified in Eq. (2.14).

2.4.3 Subloop renormalization

One-loop counterterms for subloop renormalization enter the self-energies \(\mathbf {\Sigma }_{h H A} ^{(2)}\) in Eq. (2.8). In contrast to the previously calculated \({\mathcal {O}}{\left( \alpha _t^2\right) }\) corrections, the approximation of massless bottom quarks is dropped in the present calculation. Accordingly, new counterterms for the bottom–sbottom sector are induced, which are specified in the following.

The squark mass matrices in the \(\big ({\tilde{q}}_{\text {L}},{\tilde{q}}_{\text {R}}\big )\) bases, \(q=t,b\), in the gauge-less limit are given by

$$\begin{aligned} {\mathbf {M}}_{{\tilde{q}}}&= \begin{pmatrix} m_{{\tilde{q}}_{\text {L}}}^{2} + m_{q}^{2} &{} m_{q}\left( A_{q}^{*} - \mu \,\kappa _q \right) \\ m_{q}\left( A_{q} - \mu ^{*}\,\kappa _q \right) &{} m_{{\tilde{q}}_{\text {R}}}^{2} + m_{q}^{2} \end{pmatrix}, \nonumber \\ \kappa _{t}&= \frac{1}{t_{\beta }},\quad \kappa _{b} = t_{\beta } . \end{aligned}$$
(2.16)

SU(2)-invariance requires \(m_{{\tilde{t}}_{\text {L}}}^{2} = m_{{\tilde{b}}_{\text {L}}}^{2} \equiv m_{{\tilde{Q}}_{3}}^{2}\). The squark mass eigenvalues can be obtained by performing unitary transformations,

$$\begin{aligned} {\mathbf {U}}_{{\tilde{q}}}{\mathbf {M}}_{{\tilde{q}}}{\mathbf {U}}_{{\tilde{q}}}^{\dagger }&= \mathrm {diag}{\left( m_{{\tilde{q}}_{1}}^{2},\, m_{{\tilde{q}}_{2}}^{2}\right) }. \end{aligned}$$
(2.17)

The independent parameters entering the two-loop calculation via the quark–squark sector are the quark masses \(m_{q}\), the soft SUSY-breaking parameters \(m_{{\tilde{Q}}_{3}}\) and \(m_{{\tilde{q}}_{\text {R}}}\), \({\tilde{q}} = {\tilde{t}},{\tilde{b}}\), the complex trilinear couplings \(A_{q} = | A_{q}|\mathrm {e}^{i\,\phi _{A_{q}}}\), \(q=t,b\), the complex \(\mu \) parameter from the superpotential, and the ratio of the vacuum expectation values \(t_{\beta }\). All of them have to be renormalized at the one-loop level,

$$\begin{aligned} m_q \rightarrow m_q + \delta ^{(1)}m_q,\quad {\mathbf {M}}_{{\tilde{q}}} \rightarrow {\mathbf {M}}_{{\tilde{q}}} + \delta ^{(1)}{\mathbf {M}}_{{\tilde{q}}}. \end{aligned}$$
(2.18)

Here \(\delta ^{(1)}{\mathbf {M}}_{{\tilde{q}}}\) denotes the matrix of counterterms after applying the renormalization transformation to the parameters in Eq. (2.16). The renormalization of the top–stop sector, as well as of \(\mu \) and \(t_{\beta }\), is carried out as specified in Ref. [71].

For the renormalization of the bottom–sbottom sector, we refer to Refs. [39, 50, 96] where renormalization of \(m_b\) and \(A_b\) in the \(\overline{\text {DR}}\) scheme has been proposed to avoid numerical instabilities. Also for the applied resummation of \(\Delta m_b\) the \(\overline{\text {DR}}\) scheme for \(m_b\) is convenient, as explained above. The renormalization scale is chosen to be the on-shell top mass.

  • The bottom-quark self-energy in a Lorentz decomposition is given by

    $$\begin{aligned} \Sigma _b (p)&=\, {\not \!\!{p}}\, \omega _-\, \Sigma _b^{\mathrm {L}} (p^2) + {\not \!\!{p}}\, \omega _+\, \Sigma _b^{\mathrm {R}}(p^2) \nonumber \\&\quad + m_b \,\Sigma _b^{\mathrm {S}}(p^2) + m_b \gamma _5\, \Sigma _b^{\mathrm {PS}}(p^2), \end{aligned}$$
    (2.19)

    with the left-vector part \(\Sigma _b^{\mathrm {L}}\), right-vector part \(\Sigma _b^{\mathrm {R}}\), scalar part \(\Sigma _b^{\mathrm {S}}\), and pseudo-scalar part \(\Sigma _b^{\mathrm {PS}}\). The bottom-quark mass renormalization is fixed at the on-shell top-mass scale via

    $$\begin{aligned} \delta ^{(1)}m_{b}&= m_{b} \, \mathfrak {R}{\mathfrak {e}}{\left[ \frac{1}{2}\left( \Sigma _{b}^{\text {L}}{\left( m_{b}^{2}\right) } + \Sigma _{b}^{\text {R}}{\left( m_{b}^{2}\right) }\right) + \Sigma _{b}^{\text {S}}{\left( m_{b}^{2}\right) }\right] }_{\overline{\text {DR}}}. \end{aligned}$$
    (2.20)
  • With Eqs. (2.17) and (2.18) we define

    $$\begin{aligned} {\mathbf {U}}_{{\tilde{q}}}\, \delta {\mathbf {M}}_{{\tilde{q}}}\, {\mathbf {U}}_{{\tilde{q}}}^{\dagger }&= \begin{pmatrix} \delta ^{(1)}m_{{\tilde{q}}_{1}}^{2} &{} \delta ^{(1)}m_{{\tilde{q}}_{1}{\tilde{q}}_{2}}^{2} \\ \delta ^{(1)}m_{{\tilde{q}}_{1}{\tilde{q}}_{2}}^{2\,*} &{} \delta ^{(1)}m_{{\tilde{q}}_{2}}^{2} \end{pmatrix} . \end{aligned}$$
    (2.21)

    The renormalization of the soft-breaking parameter \(A_b\) follows from Eqs. (2.16) and (2.21) with \(q=b\), yielding

    $$\begin{aligned} \delta ^{(1)}A_{b}&= \left[ {\mathbf {U}}_{{\tilde{b}}\,11}{\mathbf {U}}_{{\tilde{b}}\,12}^{*}\frac{\delta ^{(1)}M_{{\tilde{b}}_{1}}^{2} - \delta ^{(1)}M_{{\tilde{b}}_{2}}^{2}}{m_b}\right. \nonumber \\&\quad \left. + {\mathbf {U}}_{{\tilde{b}}\,21}{\mathbf {U}}_{{\tilde{b}}\,12}^{*}\frac{\delta ^{(1)}M_{{\tilde{b}}_{1}{\tilde{b}}_{2}}^{2}}{m_b}+ {\mathbf {U}}_{{\tilde{b}}\,11}{\mathbf {U}}_{{\tilde{b}}\,22}^{*}\frac{\delta ^{(1)}M_{{\tilde{b}}_{1}{\tilde{b}}_{2}}^{2\,*}}{m_b}\right. \nonumber \\&\quad \left. \phantom {\left[ {\mathbf {U}}_{{\tilde{b}}\,11}{\mathbf {U}}_{{\tilde{b}}\,12}^{*}\frac{\delta ^{(1)}M_{{\tilde{b}}_{1}}^{2} - \delta ^{(1)}M_{{\tilde{b}}_{2}}^{2}}{m_b}\right. }- \left( A_{b} - \mu ^{*}\,t_\beta \right) \frac{\delta ^{(1)}m_{b}}{m_b} + \mu ^*\,\delta t_\beta + t_\beta \,\delta \mu ^*\right] _{\overline{\text {DR}}}. \end{aligned}$$
    (2.22)

    The divergent parts of the counterterms \(\delta ^{(1)}M_{{\tilde{b}}_{i}}^{2}\), \(i=1,2\), which are needed for the \(\overline{\text {DR}}\) renormalization in Eq. (2.22), can be computed from the corresponding sbottom self-energies \(\Sigma ^{(1)}_{{\tilde{b}}_i{\tilde{b}}_i}\), \(i=1,2\), and the divergent part for the counterterm \(\delta ^{(1)}M_{{\tilde{b}}_{1}{\tilde{b}}_{2}}^{2}\) can be obtained from the sbottom mixing \(\Sigma ^{(1)}_{{\tilde{b}}_1{\tilde{b}}_2}\), where the self-energy is defined with an incoming \({\tilde{b}}_2\) and an outgoing \({\tilde{b}}_1\). The renormalization conditions to fix the auxiliary counterterms \(\delta ^{(1)}M_{{\tilde{b}}_{i}}^{2}\) and \(\delta ^{(1)}M_{{\tilde{b}}_{1}{\tilde{b}}_{2}}^{2}\) may be chosen analogously to the stop-sector counterterms in Ref. [71]. Again, the renormalization scale in Eq. (2.22) is the on-shell top mass.

  • Invariance under SU(2) yields the following relation between the stop and sbottom sector:

    $$\begin{aligned} \delta ^{(1)}m_{{\tilde{Q}}_3}^2&\equiv \sum \limits _{i\;=\;1}^2|{\mathbf {U}}_{{\tilde{t}}\,i1}|^{2}\,\delta ^{(1)}m_{{\tilde{t}}_{i}}^{2}\nonumber \\&\quad + 2\mathfrak {R}{\mathfrak {e}}{\left[ {\mathbf {U}}_{{\tilde{t}}\,21}{\mathbf {U}}_{{\tilde{t}}\,11}^{*}\,\delta ^{(1)}m_{{\tilde{t}}_{1}{\tilde{t}}_{2}}^{2}\right] } - 2\,m_t\,\delta ^{(1)}m_t\nonumber \\&= \sum \limits _{i\;=\;1}^2|{\mathbf {U}}_{{\tilde{b}}\,i1}|^{2}\,\delta ^{(1)}m_{{\tilde{b}}_{i}}^{2}\nonumber \\&\quad + 2\mathfrak {R}{\mathfrak {e}}{\left[ {\mathbf {U}}_{{\tilde{b}}\,21}{\mathbf {U}}_{{\tilde{b}}\,11}^{*}\,\delta ^{(1)}m_{{\tilde{b}}_{1}{\tilde{b}}_{2}}^{2}\right] } - 2\,m_b\,\delta ^{(1)}m_b\ . \end{aligned}$$
    (2.23)

    We trade \(\delta ^{(1)}m_{{\tilde{Q}}_3}^2\), \(\delta ^{(1)}m_{{\tilde{t}}_{\text {R}}}^2\) and \(\delta ^{(1)}A_t\) for \(\delta ^{(1)}m_{{\tilde{t}}_i}^2\), \(i=1,2\), and \(\delta ^{(1)}m_{{\tilde{t}}_1{\tilde{t}}_2}^2\), and we apply on-shell conditions to both stop particles and the stop mixing angle (see Ref. [71]). Then we choose to make \(\delta ^{(1)}m_{{\tilde{b}}_{1}}^{2}\) a dependent quantity by the relation in Eq. (2.23). The other diagonal sbottom-mass counterterm \(\delta ^{(1)}m_{{\tilde{b}}_{2}}^{2}\) is employed instead of \(\delta ^{(1)}m_{{\tilde{b}}_{\text {R}}}^{2}\) and is fixed on-shell via

    $$\begin{aligned} \delta ^{(1)}m_{{\tilde{b}}_{2}}^{2}&= \mathfrak {R}{\mathfrak {e}}{\left[ \, \Sigma _{{\tilde{b}}_{22}}^{(1)}{\left( m_{{\tilde{b}}_{2}}^{2}\right) }\right] }. \end{aligned}$$
    (2.24)

    The quantity \(\delta ^{(1)}m_{{\tilde{b}}_{1}{\tilde{b}}_{2}}^{2}\) is the off-diagonal entry of Eq. (2.21) for \(q=b\). It is already fixed by the renormalization condition of Eq. (2.22) for the independent counterterm \(\delta ^{(1)}A_b\).

3 Numerical results for the Higgs spectrum

In the following numerical analysis the new contributions of \({\mathcal {O}}{\left( \alpha _t^2+\alpha _t\alpha _b+\alpha _b^2\right) }\) are added to the known Higgs-mass corrections in the general case of the MSSM with complex parameters which are implemented in FeynHiggs (version 2.12.0).Footnote 1 While the improvement by resummation of leading logarithms as described in Refs. [57, 58, 60] can be applied also to the case of complex parameters (via an interpolation routine), we have not included contributions of this kind in our numerical results presented below.Footnote 2 The large impact of the \({\mathcal {O}}{\left( \alpha _t^2\right) }\) terms has been investigated in Refs. [70, 71] and is not presented here again. Instead the focus is set on the new corrections induced by the finite bottom mass. If not stated otherwise, we choose the following default setting for the parameters entering through the new contributions:

$$\begin{aligned} t_\beta= & {} 50, \quad m_{H^\pm } = 1.5\,\text {TeV},\quad m_{{\tilde{Q}}_3} = 2.1\,\text {TeV}, \nonumber \\ m_{{\tilde{t}}_\text {R}}= & {} m_{{\tilde{b}}_\text {R}} = 2\,\text {TeV}, \quad m_t = 173.2\,\text {GeV}, \end{aligned}$$
(3.1a)
$$\begin{aligned} A_t= & {} \left| 1.3\,m_{{\tilde{t}}_\text {R}} + \frac{\mu ^*}{t_\beta }\right| \mathrm {e}^{{i}\,\phi _{A_t}}, \quad A_b = 2.5\,m_{{\tilde{b}}_\text {R}}\,\mathrm {e}^{{i}\,\phi _{A_b}}, \nonumber \\ M_3= & {} 2.5\,\text {TeV}\,\mathrm {e}^{{i}\,\phi _{M_3}}, \quad \mu = \text {sgn}{\left[ \mu \right] }\,1\,\text {TeV}. \end{aligned}$$
(3.1b)

The quantities in Eq. (3.1a) are real parameters. The charged Higgs mass \(m_{H^\pm }\) is chosen as an input parameter, and its value is set to ensure the compatibility of scenarios with high \(t_\beta \) with the current experimental constraints from searches for heavy MSSM Higgs bosons [97, 98]. The parameters in Eq. (3.1b) are in general complex. Their respective phases \(\phi _{A_t}\), \(\phi _{A_b}\) and \(\phi _{M_3}\) are scanned in Sect. 3.2. Thereby the gluino mass parameter \(M_3\) does not occur directly in the new Higgs self-energy contributions, but it appears in the leading term of the bottom-mass resummation. The parameter \(\mu \) is also complex in general, but its phase is constrained to be very close to zero or \(\pi \) by EDM limits (see above). We remark that the phases \(\phi _{M_3}\), \(\phi _{A_t}\) and \(\phi _{A_b}\) are also constrained by EDM limits, but scenarios with large phases are possible (see e.g. Ref. [99]). We show results for the Higgs mass when varying two phases at the same time.

The absolute value of \(A_t\) has been fixed to yield a lightest Higgs-boson mass close to 125 GeV which can then be identified with the Higgs signal discovered at ATLAS and CMS. Together with \(m_{{\tilde{Q}}_3}\) and \(m_{{\tilde{t}}_\text {R}}\) it determines the mass shift which is induced by the stop contributions. We choose different values for \(m_{{\tilde{Q}}_3}\) and \(m_{{\tilde{t}}_\text {R}}\), \(m_{{\tilde{b}}_\text {R}}\) to avoid numerical instabilities due to degeneracies. Different setups for \(m_{{\tilde{t}}_\text {R}}\) and \(X_t = A_t - \mu ^*/t_\beta \) are possible to yield a lightest Higgs mass of 125 GeV as can be seen in Fig. 1. Therein the gray bar indicates the mass range \(125.1\pm 0.21(\text {stat})\pm 0.11(\text {syst})\) GeV as measured by ATLAS and CMS [100].

Fig. 1
figure 1

Dependence of the lightest Higgs-mass \(M_h\) on the stop parameters \(m_{{\tilde{t}}_{\text {R}}}\) and \(X_t=A_t-\mu ^*/t_\beta \) as predicted by FeynHiggs-2.13.0 in the version for real parameters without contributions of resummed logarithms. We set \(m_{{\tilde{Q}}_3}=m_{{\tilde{t}}_{\text {R}}}+100\) GeV. The other parameters except \(m_{{\tilde{Q}}_3}\), \(m_{{\tilde{t}}_{\text {R}}}\) and \(A_t\) have been fixed to the values given in Eq. (3.1) with vanishing phases. Two-loop corrections of \({\mathcal {O}}{\left( \alpha _t\alpha _s+\alpha _b\alpha _s+\alpha _t^2+\alpha _t\alpha _b+\alpha _b^2\right) }\) are comprised as implemented in the version of FeynHiggs for real parameters

The absolute value of \(A_b\) is close to the upper limit

$$\begin{aligned} |A_b|^2&< 3 \left( m_{H_d}^2 + |\mu |^2 + m_{{\tilde{Q}}_3}^2 + m_{{\tilde{b}}_\text {R}}^2\right) , \end{aligned}$$
(3.2a)
$$\begin{aligned} m_{H_d}^2 + |\mu |^2&= \left( m_{H^\pm }^2 - m_W^2\right) {\sin ^2}\beta - \frac{1}{2} m_Z^2 \cos {2\beta }\ , \end{aligned}$$
(3.2b)

from the approximate bound from the requirement of vacuum stability to avoid charge- and color-breaking minima [101, 102] (see Refs. [103,104,105,106,107,108,109] for more detailed discussions of this issue).

In the following analysis we call \(\Delta M_h\) the shift of the lightest Higgs-boson mass by the new Yukawa terms of \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\), i.e. excluding the previously analyzed contributions of \({\mathcal {O}}{\left( \alpha _t^2\right) }\). In Sect. 3.1, the impact of different parameters on the lightest Higgs boson mass in the CP-conserving case is investigated.

We have also investigated the mass shifts of the heavier neutral Higgs bosons. In general, the shifts are of the same absolute size as for the lightest Higgs but with opposite sign. However, since the tree-level input value \(m_{H^\pm }\) needs to be large for high values of \(t_\beta \) (where the \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) contributions are relevant) to be in agreement with experimental constraints, the relative mass shift for the heavy Higgs bosons is only \(\approx 1\permille \). Moreover, the two heavy Higgs bosons receive nearly identical corrections; in the investigated scenarios the largest difference was \(\approx 0.1\) GeV. For this reason we do not present numerical results for the mass shifts of the heavier Higgs bosons here. It should, however, be noted that even small mass shifts can have an important impact on the resonance-type behavior that typically occurs between the two heavy neutral Higgs states in CP-violating scenarios; see Refs. [110, 111].

Fig. 2
figure 2

Comparison of the lightest Higgs-boson mass \(M_h\) as predicted with our new two-loop corrections (solid), the version of FeynHiggs for the MSSM with real parameters, i.e. including \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) corrections (dashed), and the version of FeynHiggs for the MSSM with complex parameters, i.e. without \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) corrections (dotted) for the parameters specified in Eq. (3.1) with vanishing phases and \(\text {sgn}{\left[ \mu \right] }=-1\)

3.1 Scenarios with real parameters

We start to analyze our results by performing a comparison with the previously implemented two-loop corrections in FeynHiggs. The two-loop corrections of \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) were up to now only known for the MSSM with real parameters and \(m_A\) being an input. We compare our new result with the predictions obtained so far with FeynHiggs from both the versions for real parametersFootnote 3 and for complex parameters (for the latter employing the same renormalization scheme with \(m_A\) as input and in the limit of real parameters, but without terms of \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\)). Both versions contain shifts due to \(\Delta m_b\) effects (see Eq. (2.15)), including contributions of \({\mathcal {O}}{\left( \alpha _t\alpha _b\right) }\).

In Fig. 2 the predictions of FeynHiggs (dashed: MSSM with real parameters, dotted: MSSM with complex parameters) are compared to our new result (solid) as a function of \(m_A\). The different colors correspond to different values of \(t_\beta \). The large deviations between the dashed and dotted curves for large values of \(t_\beta \) are induced by the \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) terms, which are not incorporated in the dotted curve. After adding our new contributions to the result for complex parameters the agreement with the FeynHiggs result for the case of real parameters is very good, i.e. the dashed and solid lines almost coincide with each other. Since the version of the \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) corrections which is implemented in FeynHiggs so far employs further approximations of \(t_\beta \rightarrow \infty \) and \(m_b\rightarrow 0\) (see Ref. [39]), while our new result is not simplified further, the agreement is not expected to be perfect. The largest difference of \(\approx 0.3\) GeV is found in the threshold region at \(m_A=m_{{\tilde{t}}_2}-m_{{\tilde{t}}_1}\approx 200\) GeV which enters via the renormalization in the stop sector.

Fig. 3
figure 3

Dependence of the lightest Higgs-mass shift \(\Delta M_h\) on \(t_\beta \). The parameter \(\mu \) is either positive (blue) or negative (red)

In our following analysis we choose \(m_{H^\pm }\) as an input parameter. In this case the \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) terms are new contributions. We investigate the dependence of the prediction for \(M_h\) on \(t_\beta , \mu \) and \(M_3\), whereby all parameters are still kept real. The results are depicted in Figs. 3, 4, 5 and 6.

Fig. 4
figure 4

Dependence of the lightest Higgs-mass shift \(\Delta M_h\) on \(\mu \). The parameter \(M_3\) is either positive (blue) or negative (red). The region around \(\mu =0\) is left out because of numerical instabilities

Fig. 5
figure 5

Dependence of the lightest Higgs-mass shift \(\Delta M_h\) on \(M_3\). The parameter \(\mu \) is either positive (blue) or negative (red)

Fig. 6
figure 6

Dependence of the lightest Higgs-mass shift \(\Delta M_h\) on \(|A_b|\). The sign of \(A_b\) is either positive (blue) or negative (red), and \(\text {sgn}{\left[ \mu \right] } = -1\)

Fig. 7
figure 7

Dependence of the lightest Higgs-mass shift \(\Delta M_h\) on \(\phi _{A_t}\) and \(\phi _{A_b}\), \(\text {sgn}{\left[ \mu \right] } = -1\). solid: exact calculation, dotted: interpolation in FeynHiggs, the red-dotted and orange-dotted lines are identical

As can be seen in Fig. 3 large contributions above 1 GeV are only visible at high values of \(t_\beta \). In this scenario \(M_3\) is positive, leading to a much bigger \(\Delta M_h\) if \(\mu \) is negative, which can be understood from Eqs. (2.14) and (2.15). For later analysis we fix \(t_\beta = 50\).

In Fig. 4 we investigate the dependence of \(\Delta M_h\) on the size of \(\mu \). We find very large gradients for the following two cases: positive \(M_3\) and negative \(\mu \approx -1.8\) TeV, and negative \(M_3\) and positive \(\mu \approx 2.6\) TeV, which can again be understood from Eqs. (2.14) and (2.15), where for too large values of \(|\mu |\) and opposite signs of \(\mu \) and \(M_3\) the perturbative region of parameter space is left, as \(\Delta m_{b}\rightarrow -1\). A further increase of \(|\mu |\) in the regions of large gradients leads to a very strong enhancement of the bottom Yukawa coupling and accordingly to very large negative mass shifts, yielding eventually a tachyonic Higgs boson. For the following analysis, we choose to fix \(\mu =-1\) TeV, i.e. below the problematic scale and with \(\text {sgn}{\left[ \mu \right] }=-1\). However, it should be noted that scenarios with positive \(\mu \) can lead to large shifts as well, when \(M_3\) is negative, as in both cases the bottom Yukawa coupling is enhanced. Moreover, scenarios with \(\text {sgn}{\left[ \mu \right] }=1\) are in better agreement with constraints from the anomalous magnetic moment of the muon [112,113,114]. Close to \(|\mu |= m_{{\tilde{t}}_{1,2}} - m_t \approx 1.8\) TeV one can see kinks which are induced by threshold effects from the higgsino–top–stop system.

In Fig. 5 the impact of the gluino mass parameter is depicted. This effect enters the Higgs self-energies at the investigated order purely via the employed effective bottom mass. We see a rising shift at growing \(|M_3|\) for opposite signs of \(\mu \) and \(M_3\) (yielding the same enhancement in \(\Delta m_b\), see Eqs. (2.15)). At \(|M_3|= m_{{\tilde{b}}_{1,2}} - m_b \approx 2\) TeV (nearly invisible) threshold effects from the gluino–bottom–sbottom system appear. For our following analysis we fix \(|M_3|\) above that region at 2.5 TeV.

Finally, in Fig. 6 the absolute value of \(A_b\) is varied, and the resulting mass shift is plotted for positive sign (blue) and negative sign (red) of \(A_b\). The difference between both curves, i.e. the impact of the phase \(\phi _{A_b}\), is enhanced for larger absolute values. However, as too large values of \(|A_b|\) lead to instable vacua according to the upper limit of Eq. (3.2), we set it to \(|A_b|= 2.5\, m_{{\tilde{b}}_\text {R}}\) in the scenarios of the following section.

3.2 Scenarios with complex parameters

Various phases enter the self-energies of the Higgs bosons at \({\mathcal {O}}{\left( \alpha _t^2+\alpha _t\alpha _b+\alpha _b^2\right) }\). Their impact on the Higgs sector is shown in Figs. 7, 8 and 9. Here we keep \(\mu \) negative, i.e. \(\text {sgn}{\left[ \mu \right] }=-1\), and \(M_3\) positive, but we could also have chosen the opposite signs of both parameters to see enhanced effects for the phase dependent terms as has been shown in Fig. 4.

We start with the phases \(\phi _{A_t}\) and \(\phi _{A_b}\). The results are depicted in Fig. 7, where mass shifts between 0.3 and 1.4 GeV can be seen. For \(\phi _{A_b} = 0\) the variation with respect to \(\phi _{A_t}\) is maximal; the larger the phase of \(A_b\), the flatter the dependence on the phase of \(A_t\). Similarly, variation of \(\phi _{A_b}\) yields the largest effects for \(\phi _{A_t} = 0\). Also the signs of the phases matter, e.g. the mass shifts are different for \(\phi _{A_b}=\pm \tfrac{\pi }{2}\).

In addition to the exact calculation (solid lines), FeynHiggs offers an implemented interpolation of the self-energy corrections that have been known up to now for the case of real parameters but not for the complex case. Since the \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) terms were only available in FeynHiggs for the MSSM with real parameters in the limits \(t_\beta \rightarrow \infty \) and \(m_b\rightarrow 0\), deviations from the new mass shifts can be expected even for real parameters. Besides these relatively small differences, the linear interpolation can differ by \(\approx 0.5\) GeV from the full result in the investigated scenario. Also the asymmetric behavior for the change of two phases at the same time was not described correctly by the interpolation.

Fig. 8
figure 8

Dependence of the lightest Higgs-mass shift \(\Delta M_h\) on \(\phi _{M_3}\) and \(\phi _{A_b}\), \(\text {sgn}{\left[ \mu \right] } = -1\)

Fig. 9
figure 9

Dependence of the lightest Higgs-mass shift \(\Delta M_h\) on \(\phi _{M_3}\) and \(\phi _{A_t}\), \(\text {sgn}{\left[ \mu \right] } = -1\)

Figures 8 and 9 show the influence of varying the gluino phase \(\phi _{M_3}\) and in addition either \(\phi _{A_b}\) or \(\phi _{A_t}\). These terms are induced by the correction factor \(\Delta m_b\) as the investigated class of two-loop corrections does not contain the parameter \(M_3\). Also here the largest phase dependence is found when one phase is equal to zero. In Fig. 8 the mass shift is nearly symmetric in \(\pm \phi _{M_3}\) and \(\pm \phi _{A_b}\), i.e. the red and yellow curves are lying on top of each other. Nevertheless, there are small asymmetries in the renormalized two-loop self-energies \({\hat{\Sigma }}_{hA}\) and \({\hat{\Sigma }}_{HA}\). On the contrary the mass shift \(\Delta M_h\) in Fig. 9 shows a clear asymmetry similar to Fig. 7.

In summary, phase dependent contributions of \({\mathcal {O}}\left( \alpha _t\alpha _b+\alpha _b^2\right) \) lead to mass shifts of the lightest Higgs boson of \(\approx 1\) GeV in the investigated scenarios. The sign of \(\mu \) has been chosen to be negative in the considered scenarios, but similar effects can be found at positive large \(\mu \) (and opposite sign of \(M_3\)).

4 Conclusions

The two-loop corrections of \({\mathcal {O}}{\left( \alpha _t^2+\alpha _t\alpha _b+\alpha _b^2\right) }\) to the Higgs-boson masses in the MSSM with complex parameters have been computed in the gauge-less limit at vanishing external momentum. The terms of \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) have only been known in the special case of the MSSM with real parameters before, and were incorporated in FeynHiggs in the limits \(t_\beta \rightarrow \infty \) and \(m_b\rightarrow 0\). The specific aspects related to the renormalization of these new contributions have been discussed, and their numerical impact on the Higgs spectrum has been investigated.

For the lightest Higgs boson mass at \(\approx 125\) GeV we have found shifts above 1 GeV at \(t_\beta >40\) for different scenarios: moderate \(|\mu |=1\) TeV with negative sign and positive \(M_3\), \(A_t\), \(A_b\), or with positive sign and negative \(M_3\), \(A_t\), \(A_b\). The reason for that enhancement can be found in the large correction factor \(\Delta m_b\) yielding an enhancement of the bottom Yukawa coupling. The effect of varying the phases \(\phi _{M_3}\), \(\phi _{A_t}\) and \(\phi _{A_b}\) can be as large as 1 GeV. If one phase is set close to \(\pi \), the dependence on the other phases is typically weakened; the largest effects are found when only one phase is varied with all others being zero. In FeynHiggs so far an interpolation of the corrections of \({\mathcal {O}}{\left( \alpha _t\alpha _b+\alpha _b^2\right) }\) obtained for the case of real parameters is used for the case of complex parameters. We have found deviations with a size of \(\approx 0.5\) GeV from this approximation, especially when several phases are different from zero at the same time.

Mass shifts for the heavier neutral Higgs bosons have not been depicted. They are similar to the ones of the lightest Higgs boson; however, with opposite sign. Since we used a large value for \(t_\beta \) in our scenarios, we need to choose a rather large input mass \(m_{H^\pm }=1.5\) TeV to be consistent with existing experimental bounds. Therefore, the relative size of the mass shifts is small. Moreover, the two heavy Higgs bosons receive similar corrections with a maximal difference of \(\approx 0.1\) GeV in the investigated scenarios. Nevertheless, small mass shifts can be important to correctly describe the resonance-type behavior of nearly mass-degenerate mixed states like the two heavy Higgs bosons in the MSSM with complex parameters.

The new results will be implemented in the public code FeynHiggs.