1 Introduction

Supersymmetric extensions of the Standard Model (SM) can help overcome shortcomings of the SM, providing for example dark matter candidates and additional sources of CP violation. While these models are theoretically attractive, none of the predicted supersymmetric (SUSY) particles has been found thus far. In order to guide experimentalists in their continued hunt for these yet-undiscovered states and to ascertain the continued theoretical viability of these theories, precise theoretical predictions of experimentally observable quantities are needed.

The discovery of the Higgs boson in 2012 [1, 2] has provided theorists with a whole new set of results against which these theories may be tested. The best measured property of the Higgs boson is its mass \(m_h\) with a value of [3]

$$\begin{aligned} m_h = 125.09 \pm 0.24 \text { GeV} . \end{aligned}$$
(1)

This mass is a free parameter in the SM, and must be determined from experiment. In supersymmetric extensions, however, the mass of the discovered Higgs boson can be predicted using the additional parameters of the theory. Demanding that the theoretical prediction match the experimentally-measured value for each point in the parameter space constrains the theory. Since the experimental result for the Higgs mass is very precise, it is necessary to obtain a theoretical prediction that is of similar accuracy in order to fully exploit the information of the experiment and to yield stringent constraints. This entails incorporating quantum corrections in the Higgs mass calculations.

The necessity of including quantum corrections of higher order in the Higgs-mass prediction has been recognized for a long time. In the context of the Minimal Supersymmetric Standard Model (MSSM), radiative corrections have been shown to be important to lift the mass of the lightest Higgs boson above the tree-level upper limit, which is given by the Z-boson mass [4,5,6,7]. An enormous effort has been devoted to improve the theoretical prediction: At fixed order, one-loop [8,9,10,11,12,13], two-loop [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] as well as three-loop corrections [41,42,43,44,45,46,47] to the Higgs-boson masses have been calculated. As it stands, the theoretical uncertainty of the lightest Higgs boson has been estimated to be of the order 3 GeV [48], an estimate which is still applied in phenomenological studies. While the assessment of the theoretical uncertainty is still an ongoing discussion (see Refs. [49,50,51,52,53,54,55]), the general consensus is that it is challenging to reduce the uncertainty below 1 GeV.

The method by which one includes these quantum corrections is dependent on the mass of the SUSY particles. The aforementioned fixed-order approach is particularly useful for SUSY particles with masses up to the TeV scale where the mass of the SUSY partners of the top quark is most important. SUSY particles at this scale have been motivated by naturalness and grand-unification arguments among others, but as there are still no signs of low-energy SUSY at the LHC, heavy SUSY particles have begun to attract more interest. Within these scenarios, the fixed-order calculations provide a less accurate result due to the emergence of large logarithms of ratios of masses of SUSY particles and the energy scale at which the calculation is performed. These large logarithms spoil the fixed-order perturbation series. In order to take these large logarithms into account, an effective-field theory approach has been applied. In this approach, the full MSSM governs the interaction behaviour at the high-energy scale, and the effects of the heavy SUSY particles on the low-energy physics are encoded into the couplings of a viable low-energy theory such as the SM or the Two-Higgs-Doublet Model (2HDM) via matching at the matching scale. This matching entails calculating the couplings in both the high and low energy theories up to fixed order such that the physics described by the MSSM and the 2HDM at the matching scale is the same. The resulting couplings are then evolved down to the low-energy scale at which the observables are investigated with the help of the corresponding renormalization group equations (RGE). Solving the renormalization group equations leads to a resummation of the large logarithms. Assuming the SM as low-energy theory, leading logarithms (LL) have been resummed in Refs. [56, 57] and next-to leading logarithms (NLL) in Refs. [58,59,60,61] taking into account top-Yukawa corrections. Allowing all Higgs bosons to be light, LL corrections have been calculated in Refs. [62,63,64,65]. Analytical expressions for the Higgs-boson mass have been obtained taking LL contributions into account up to two-loop order [66]. A hierarchical stop-mass spectrum with one stop-quark mass much heavier than the other has been considered in Ref. [67] and the resulting two-loop LL and NLL corrections of the Higgs-boson mass have been calculated. This approach has also been used to obtain one- and two-loop LL corrections to the Higgs-boson mass spectrum in a CP-violating scenario [68, 69]. In Ref. [70], this approach has been performed taking into account an intermediate step with light electroweak fermionic SUSY partners and the SM as low-energy theory. The results have been improved further allowing also for light Higgs bosons in Refs. [71, 72].Footnote 1 The calculations of Refs. [70, 71] have been implemented into the tool MhEFT. A further tool was presented with SUSYHD [49], which includes threshold corrections to a higher order. The threshold correction to the quartic Higgs coupling, assuming the SM as low-energy theory, has been calculated taking into account two-loop QCD [74], third-generation Yukawa [75], both in the gaugeless limit, and the full QCD corrections [76]. Effects of higher-dimensional operators have been studied in Refs. [75, 77] where Ref. [77] exploits a different method to obtain the one-loop threshold corrections. The prediction of the Higgs-boson mass has been improved by resumming logarithms of fourth logarithmic order (\(\text {N}^3\)LL) [78] in the case that the SM is the low-energy theory. Furthermore, a combination of both approaches, the fixed-order and the RGE approach, has been performed, in particular to improve the intermediate regime with particles not very heavy but too heavy for the existing fixed-order results [79,80,81,82,83]. These results are implemented in FeynHiggs [13, 17, 48, 79, 80, 84, 85]. Similarly, FlexibleSUSY [86, 87] as well as a version of SARAH/SPheno [88,89,90,91,92,93] have implemented such a combination [94, 95]. Furthermore, in Ref. [96], \(\text {N}^3\)L0-fixed-order and \(\text {N}^3\)LL results have been combined to give a precise prediction for the mass of the lightest Higgs boson. A review about the different Higgs-mass calculations can be found in Ref. [53].

The particular theory we wish to explore in this work is the MSSM with complex parameters where the SUSY particles are heavy and the additional Higgs bosons have masses at an intermediate scale. In addition to studying the mass of the lightest neutral Higgs boson in this scenario, we analyze the size of this boson’s CP-odd component, which is induced by quantum corrections in the presence of complex parameters of the high-energy theory. We also look into the mixing and masses of the heavy Higgs bosons. For this, we assume a type-III complex 2HDM as the low-energy effective theory, where both Higgs doublets can couple to right-handed up- as well as right-handed down-type fermions and use different approaches to connect to the Standard Model.

Other studies of a CP-violating MSSM have been done. The scan of the phenomenological MSSM in Ref. [97] did not find a measurable size of the CP-odd component of the lightest Higgs boson, while the study presented in Ref. [98] found some more promising results with some scenarios that could be measurable at least at future runs of the LHC. In Ref. [99] a scenario with heavy superpartners has been explored using complex MSSM parameters to include CP-violating effects, taking into account the finding of Ref. [71].

In this paper, we further improve the results for scenarios of an MSSM with complex parameters and heavy superpartners: For the first time, we exploit two-loop RGEs for the 2HDM type III allowing for complex parameters including all Yukawa couplings of the third generation of quarks. Correspondingly, we perform a one-loop matching of the MSSM to the 2HDM type III including corrections proportional to the top Yukawa, the bottom Yukawa and the strong coupling assuming a common soft-SUSY breaking mass scale for the superpartner particles. In the final result, the corresponding NLL are resummed.Footnote 2

We will start out with setting our conventions for the MSSM in Sect. 2 and continue to describe the considered low-energy theory in Sect. 3. In Sect. 4 we discuss the exploited matching procedure and in Sect. 5, we present the details of the determination of the Higgs-boson masses as well as their mixing. The numerical results are discussed in Sect. 6, and, finally, we conclude in Sect. 7.

2 The Minimal Supersymmetric Standard Model

As mentioned in Sect. 1, we consider a scenario in which all the superpartner particles are much heavier than the SM particles and the Higgs bosons. Hence, the full MSSM is active at a high-energy scale and the effects of the superpartners enter the low-energy theory via threshold effects. The aim of this section is to set up the notation needed for the calculation of the threshold corrections.

The superpotential of the MSSM is given as

$$\begin{aligned} W_{\text {MSSM}}&= - \epsilon _{ij}\left[ h_u^{\text {MSSM}} {\check{H}}_u^i {\check{Q}}^j {\check{U}}^c - h_d^{\text {MSSM}} {\check{H}}_d^i {\check{Q}}^j {\check{D}}^c\right. \nonumber \\&\quad \left. -h_{e}^{\text {MSSM}}{\check{H}}_d^i {\check{L}}^j {\check{E}}^c + \mu {\check{H}}_d^i{\check{H}}_u^j\right] , \end{aligned}$$
(2)

where all fields denoted with a \(\check{}\) are left-chiral superfields. \({\check{Q}}\) and \({\check{L}}\) denote the quark and lepton superfield doublets, \({\check{U}}\), \({\check{D}}\) and \({\check{E}}\) denote the up-type quark, down-type quark, and lepton charge-conjugate superfield singlets, the two Higgs-doublet superfields are denoted by \({\check{H}}_u\) and \( {\check{H}}_d\), and \(\epsilon _{12} = 1\). The corresponding Yukawa couplings are \(h_u^{\text {MSSM}}\), \(h_d^{\text {MSSM}}\), and \(h_{e}^{\text {MSSM}}\), which are complex \(3\times 3\) matrices in general. However, in our calculation, we will neglect the first two generations, and the matrices collapse to the the top-, bottom-, and tau-Yukawa coupling \(h_t^{\text {MSSM}}\), \(h_b^{\text {MSSM}}\), and \(h_{\tau }^{\text {MSSM}}\), which can be chosen to be real [102]. The strength of the mixing of the two Higgs doublets is described by the complex parameter \(\mu \).

We do not explicitly give the vector part of the Lagrangian with the kinetic and interaction terms for the gauge bosons, but refer the reader to e.g. Refs. [103, 104].

Since supersymmetry cannot be exact, it is explicitly broken in the MSSM by soft-SUSY breaking terms:

$$\begin{aligned} {\mathcal {L}}_{\text {MSSM}}^{\text {soft}}&= - m_{H_d}^2 |H_d|^2 - m_{H_u}^2 |H_u|^2 \nonumber \\&\quad - M_{L_Q}^2 |{\tilde{Q}}|^2 - M_{R_U}^2 |{\tilde{u}}_R|^2 - M_{R_D}^2 |{\tilde{d}}_R|^2 \nonumber \\&\quad - M_{L_L}^2 |{\tilde{L}}|^2 - M_{R_E}^2 |{\tilde{e}}_R|^2 \nonumber \\&\quad + \epsilon _{ij} (m_{H_dH_u}^2 H_d^i H_u^j + h_u^{\text {MSSM}} A_u H_u^i {\tilde{Q}}^j {\tilde{u}}^*_R\nonumber \\&\quad - h_d^{\text {MSSM}} A_d H_d^i {\tilde{Q}}^j {\tilde{d}}^*_R \nonumber \\&\quad - h_e^{\text {MSSM}} A_e H_d^i {\tilde{L}}^j {\tilde{e}}^*_R + \mathrm {h.c.}) \nonumber \\&\quad - \frac{1}{2} (M_1 {\tilde{B}} {\tilde{B}} + M_2 {\tilde{W}}_i {\tilde{W}}_i + M_3 {\tilde{G}} {\tilde{G}} + \mathrm {h.c.}). \end{aligned}$$
(3)

where \({\tilde{Q}}\), \({\tilde{L}}\), \({\tilde{U}}\), \({\tilde{D}}\), \({\tilde{E}}\), \(H_d\) and \(H_u\) denote the scalar component of the corresponding superfields. The \(\tilde{}\) indicates a superpartner field. The gaugino fields corresponding to U(1), SU(2), and SU(3) are denoted by \({\tilde{B}}\), \({\tilde{W}}\), and \({\tilde{G}}\), respectively. We assume colour indices to be implicit. The gaugino soft-SUSY breaking parameters \(M_1\), \(M_2\), and \(M_3\) as well as the Higgs mixing parameter \(m_{H_dH_u}^2\) are complex numbers, while the soft Higgs mass breaking parameters \(m_{H_d}^2\), \(m_{H_u}^2\) are real. In general, the sfermion mass parameters \(M_{L_Q}^2\), \(M_{R_U}^2\), \(M_{R_D}^2\), \(M_{L_L}^2\), \(M_{R_E}^2\) are \(3\times 3\) Hermitian matrices, but reduce to diagonal matrices with three independent real entries when generation mixing is ignored (The generation indices are implicit.) Finally, the trilinear couplings \(A_u\), \(A_d\), and \(A_e\) are general \(3 \times 3\) complex matrices, but reduce to complex numbers if they are assumed to be proportional to the SM Yukawa matrices, as is done in this paper.

In this paper, we are concerned with the Higgs sector of the MSSM. Equations 2 and 3 together with the D-terms of the above mentioned vector part of the Lagrangian give rise to a Higgs potential of the form

$$\begin{aligned} V_H= & {} \frac{1}{8}(g^2+g_y^2)(|H_d|^2-|H_u^2|)^2 \nonumber \\{} & {} + \frac{1}{2}g^2|H_d^\dagger H_u|^2+|\mu |^2(|H_d|^2+|H_u|^2) \nonumber \\{} & {} + m_{H_d}^2|H_d|^2 + m_{H_u}^2|H_u^2| \nonumber \\{} & {} - m_{H_dH_u}^2(\epsilon _{ab}H^a_d H^b_u + \mathrm {h.c.}) \end{aligned}$$
(4)

for the two Higgs doublets \(H_u\) and \(H_d\) of hypercharge 1 and − 1, respectively. The SU(2) and the U(1) gauge coupling are denoted by g and \(g_y\), respectively. Finally, the squark mass matrices are given by

$$\begin{aligned} {\mathcal {M}}_{{\tilde{q}}}&= \begin{pmatrix} M_{L_Q}^2 + m_q^2 + M_Z^2 c_{2 \beta } (T_q^3 - Q_q \sin ^2 \theta _W) &{}\quad m_q X_q^* \\ m_q X_q &{}\quad M_{R_F}^2 + m_q^2 +M_Z^2 c_{2 \beta } Q_q \sin ^2 \theta _W \end{pmatrix} \end{aligned}$$
(5)

with

$$\begin{aligned} X_q&= A_q - \mu ^*\kappa ~, \quad \kappa = \{\cot \beta , \tan \beta \} \quad \text {and } \nonumber \\ F&= \{U, D\} \quad \textrm{for} \quad q = \{t, b\}~. \end{aligned}$$
(6)

Here, we introduce the gauge-boson mass \(M_Z\), the electroweak mixing angle \(\theta _W\), the quark masses \(m_q\), as well as \(\beta \), which is defined via the ratio of the Higgs vacuum expectation values of the MSSM, \(\tan \beta \equiv v_u/v_d\). The charge and the third component of the isospin of the squarks are denoted by \(Q_q\) and \(T_q^3\), respectively.

3 The effective low-energy theory

The resulting low-energy theory is a Two-Higgs-Doublet Model (2HDM) with the following Higgs potential V

$$\begin{aligned} V&= m_{11}^2\Phi _1^{\dagger }\Phi _1 + m_{22}^2\Phi _2^{\dagger }\Phi _2 - [m_{12}^2\Phi _1^{\dagger }\Phi _2 + \text {h.c.}]\nonumber \\&\quad + \frac{1}{2}\lambda _1(\Phi _1^{\dagger }\Phi _1)^2 + \frac{1}{2}\lambda _2(\Phi _2^{\dagger }\Phi _2)^2 \nonumber \\&\quad + \lambda _3(\Phi _1^{\dagger }\Phi _1)(\Phi _2^{\dagger }\Phi _2) + \lambda _4(\Phi _1^{\dagger }\Phi _2)(\Phi _2^{\dagger }\Phi _1)\nonumber \\&\quad + \Big \{\frac{1}{2}\lambda _5(\Phi _1^{\dagger }\Phi _2)^2 + [\lambda _6(\Phi _1^{\dagger }\Phi _1) \nonumber \\&\quad + \lambda _7(\Phi _2^{\dagger }\Phi _2)]\Phi _1^{\dagger }\Phi _2 + \text {h.c.}\Big \}. \end{aligned}$$
(7)

Here, the mass parameters \(m_{11}^2\) and \(m_{22}^2\) are real, \(m_{12}^2\) is complex, the quartic couplings \(\lambda _{1\ldots 4}\) are real, and \(\lambda _5\), \(\lambda _6\), and \(\lambda _7\) are in general complex. The two Higgs doublets \(\Phi _1\) and \(\Phi _2\), both having hypercharge \(Y=1\), can be decomposed into

$$\begin{aligned} \Phi _1&= \begin{pmatrix} \phi ^+_1 \\ \frac{1}{\sqrt{2}} (v_1 + \phi _1 + \text {i} a_1) \end{pmatrix}, \nonumber \\ \Phi _2&= \begin{pmatrix} \phi ^+_2 \\ \frac{1}{\sqrt{2}} (v_2 + \phi _2 + \text {i} a_2) \end{pmatrix} \end{aligned}$$
(8)

where \(v_1\) and \(v_2\) are the vacuum expectation values, \(\phi ^+_1\), \(\phi ^+_2\) two complex Higgs fields, and \(\phi _1\), \(\phi _2\), \(a_1\), \(a_2\) the neutral Higgs fields.

In the matching conditions, we take into account loop-induced couplings of the“wrong” Higgs doublet to the corresponding quarks, which renders the 2HDM a type III instead of the tree-level type II version, where one Higgs doublet couples to the right-handed up-type quarks and the other Higgs doublet couples to the right-handed down-type quarks and charged leptons. The Yukawa Lagrangian for the third generation is accordingly

$$\begin{aligned} {\mathcal {L}}_{\text {Yukawa}}&= h_t' \epsilon _{ij} \Phi _1^i t_c Q^j + h_t \epsilon _{ij} \Phi _2^i t_c Q^j - h_b \delta _{ij}\Phi _1^{*i} b_c Q^j \nonumber \\&\quad - h_b'\delta _{ij}\Phi _2^{*i} b_c Q^j + \text {h.c.}\,. \end{aligned}$$
(9)

Here, we follow the SUSY conventions and write all fields as left-handed fields. Q is the quark doublet, \(\epsilon _{12} = 1\), and \(t_c\) and \(b_c\) are the left-handed top- and bottom-quark charge-conjugate fields, respectively. The Yukawa couplings \(h_t\) and \(h_b\) are the top- and bottom-Yukawa couplings also present in the type II 2HDM case, while \(h_t'\) and \(h_b'\) denote the coupling to the “wrong” Higgs doublet only existing in the type III case. We neglect contributions from Yukawa couplings from the first two generations and neglect the \(h_{\tau }\) and \(h'_{\tau }\) Yukawa couplings.

The tree-level mass matrices, parameterized in terms of charged Higgs boson mass \(M_{H^\pm }\), have the entries

$$\begin{aligned} {\mathcal {M}}^{2}_{11}&= v^2 \left( c_{\beta }^2 \lambda _1 +\frac{1}{2} s_{\beta }^2 \left[ \lambda _4+{\text {Re}}(\lambda _5)\right] \right. \nonumber \\&\quad \left. + 2 c_{\beta }s_{\beta } {\text {Re}}(\lambda _6)\right) + s_{\beta }^2 M_{H^\pm }^2, \nonumber \\ {\mathcal {M}}^{2}_{12}&= v^2 \left( c_{\beta } s_{\beta } \lambda _3 + \frac{1}{2} c_{\beta } s_{\beta } \left[ \lambda _4+{\text {Re}}(\lambda _5)\right] \right. \nonumber \\&\quad \left. +c_{\beta }^2 {\text {Re}}(\lambda _6) + s_{\beta }^2 {\text {Re}}(\lambda _7)\right) -c_{\beta } s_{\beta }M_{H^\pm }^2,\nonumber \\ {\mathcal {M}}^{2}_{22}&= v^2 \left( s_{\beta }^2 \lambda _2 + \frac{1}{2} c_{\beta }^2 \left[ \lambda _4+{\text {Re}}(\lambda _5)\right] \right. \nonumber \\&\quad \left. +2 c_{\beta }s_{\beta }{\text {Re}}(\lambda _7) \right) +c_{\beta }^2 M_{H^\pm }^2, \nonumber \\ {\mathcal {M}}^{2}_{33}&= \frac{1}{2} s_{\beta }^2 \{v^2 \left[ \lambda _4-{\text {Re}}(\lambda _5)\right] +2 M_{H^\pm }^2\}, \nonumber \\ {\mathcal {M}}^{2}_{34}&= -\frac{1}{2} c_{\beta } s_{\beta } \{v^2 \left[ \lambda _4-{\text {Re}}(\lambda _5)\right] +2 M_{H^\pm }^2\},\nonumber \\ {\mathcal {M}}^{2}_{44}&= \frac{1}{2} c_{\beta }^2 \{v^2 \left[ \lambda _4-{\text {Re}}(\lambda _5)\right] +2 M_{H^\pm }^2\},\nonumber \\ {\mathcal {M}}^{2}_{13}&= \frac{1}{2} s_{\beta } v^2 \left[ s_{\beta } {\text {Im}}(\lambda _5) + 2 c_{\beta } {\text {Im}}(\lambda _6)\right] = -{\tan \beta } {\mathcal {M}}^{2}_{14}, \nonumber \\ {\mathcal {M}}^{2}_{23}&= \frac{1}{2} s_{\beta } v^2 \left[ c_{\beta } {\text {Im}}(\lambda _5)+2 s_{\beta } {\text {Im}}(\lambda _7)\right] =-{\tan \beta }{\mathcal {M}}^{2}_{24}, \end{aligned}$$
(10)
$$\begin{aligned} {\mathcal {M}}^{2}_{+_{11}}&= s_{\beta }^2 M_{H^\pm }^2, \nonumber \\ {\mathcal {M}}^{2}_{+_{12}}&= -c_{\beta } s_{\beta }M_{H^\pm }^2 = {\mathcal {M}}^{2}_{+_{21}},\nonumber \\ {\mathcal {M}}^{2}_{+_{22}}&= c_{\beta }^2 M_{H^\pm }^2 \end{aligned}$$
(11)
$$\begin{aligned} \text {with} \quad v^2\equiv \sqrt{v_1^2+v_2^2}, \end{aligned}$$

where \({\mathcal {M}}^2\) is the neutral Higgs mass matrix in the \((\phi _1, \phi _2, a_1, a_2)\) basis and \({\mathcal {M}}^{2+}\) is the charged Higgs mass matrix in the \((\phi _1^+, \phi _2^+)\) basis.

4 Matching the MSSM to the 2HDM

The complex MSSM Higgs potential of Eq. (4) is matched to the general type-III 2HDM given above at the one-loop level at the scale \(M_{\text {s}}\). The doublets \(H_u\) and \(H_d\) from Eq. (4) are related to those of Eq. (7) by

$$\begin{aligned} \Phi _1\equiv -i\sigma _2H_d^*, \qquad \Phi _2 \equiv H_u. \end{aligned}$$
(12)

At tree level, the matching conditions for the quartic Higgs couplings are

$$\begin{aligned} \lambda _1&= \frac{1}{4}(g^2+g_y^2),\nonumber \\ \lambda _2&= \frac{1}{4}(g^2+g_y^2),\nonumber \\ \lambda _3&= \frac{1}{4}(g^2-g_y^2),\nonumber \\ \lambda _4&= -\frac{1}{2}g^2, \nonumber \\ \lambda _5&= \lambda _6 = \lambda _7 = 0. \end{aligned}$$
(13)

For the Yukawa couplings, one obtains

$$\begin{aligned} h_t^{\textrm{2HDM}}= & {} h_t^{\textrm{MSSM}},\nonumber \\ h_b^{\textrm{2HDM}}= & {} h_b^{\textrm{MSSM}},\nonumber \\ h_t'^{\textrm{2HDM}}= & {} h_b'^{\textrm{2HDM}}=0. \end{aligned}$$
(14)

Here we are assuming that we know all the parameters of the MSSM and match the MSSM to the 2HDM. This means all couplings given below are assumed to be MSSM couplings (unless we explicitly state otherwise) and we drop the superscript \({}^{\textrm{MSSM}}\).

The one-loop threshold corrections are calculated under three assumptions. Firstly, all soft-SUSY breaking mass parameters are assumed to share the common mass scale \(M_{\text {s}}\), in particular \(M_{\text {s}}= M_{L_Q} = M_{R_U} = M_{R_D} = M_3\). Secondly, all Yukawa couplings are assumed to vanish except the top and bottom Yukawa couplings, and only contributions proportional to powers of these Yukawa couplings or the strong gauge coupling are included. This amounts to only including third-generation squarks in loops when deriving the thresholds and neglecting terms of \({\mathcal {O}}(g^a g_y^b)\) with \(a+b = 4\). Lastly, all loop functions are evaluated in the limit of zero external momenta. The diagrams were evaluated using FeynArts and FormCalc [105, 106], and the loop functions are evaluated using ANT [107].

The results given here agree with the complex results given in Ref. [99] up to gauge contributions to the couplings \(\lambda _6\) and \(\lambda _7\), which, however, are found in the work of Refs. [65] and [81] in the real case.

Fig. 1
figure 1

Two sample diagrams contributing to the one-loop threshold of the quartic Higgs couplings

Box diagrams like those of Fig. 1a lead to the following corrections to the quartics

$$\begin{aligned} \Delta \lambda _{1}^{(4)}&= -\frac{\kappa }{2 } \Big \{|{\hat{A}}_b|^4 h_b^4+h_t^4 |{\hat{\mu }}|^4\Big \}, \end{aligned}$$
(15)
$$\begin{aligned} \Delta \lambda _{2}^{(4)}&= -\frac{\kappa }{2 } \Big \{|{\hat{A}}_t|^4 h_t^4+h_b^4 |{\hat{\mu }}|^4\Big \}, \end{aligned}$$
(16)
$$\begin{aligned} \Delta \lambda _{3}^{(4)}&= \frac{\kappa }{2 } \Big \{-|{\hat{A}}_b|^2 |{\hat{A}}_t|^2 h_b^2 h_t^2 \nonumber \\&\quad -|{\hat{A}}_b|^2 h_b^4 |{\hat{\mu }}|^2-|{\hat{A}}_t|^2 h_t^4|{\hat{\mu }}|^2-h_b^2 h_t^2 |{\hat{\mu }}|^4 \nonumber \\&\quad +({\hat{A}}_b {\hat{A}}_t^* + {\hat{A}}_b^* {\hat{A}}_t) h_b^2 h_t^2 |{\hat{\mu }}|^2\Big \}, \end{aligned}$$
(17)
$$\begin{aligned} \Delta \lambda _{4}^{(4)}&= \frac{\kappa }{2 } \Big \{|{\hat{A}}_b|^2 |{\hat{A}}_t|^2 h_b^2 h_t^2-|{\hat{A}}_b|^2 h_b^4 |{\hat{\mu }}|^2\nonumber \\&\quad -|{\hat{A}}_t|^2 h_t^4|{\hat{\mu }}|^2+h_b^2 h_t^2 |{\hat{\mu }}|^4\nonumber \\&\quad - ({\hat{A}}_b {\hat{A}}_t^* + {\hat{A}}_b^* {\hat{A}}_t) h_b^2 h_t^2 |{\hat{\mu }}|^2 \Big \}, \end{aligned}$$
(18)
$$\begin{aligned} \Delta \lambda _{5}^{(4)}&= -\frac{{\hat{\mu }}^2}{2 }\kappa \Big \{{\hat{A}}_b^2 h_b^4 +{\hat{A}}_t^2 h_t^4\Big \}, \end{aligned}$$
(19)
$$\begin{aligned} \Delta \lambda _{6}^{(4)}&= \frac{{\hat{\mu }}}{2 } \kappa \Big \{|{\hat{A}}_b|^2 {\hat{A}}_b h_b^4 + {\hat{A}}_t h_t^4 |{\hat{\mu }}|^2\Big \}, \end{aligned}$$
(20)
$$\begin{aligned} \Delta \lambda _{7}^{(4)}&= \frac{{\hat{\mu }}}{2 } \kappa \Big \{|{\hat{A}}_t|^2 {\hat{A}}_t h_t^4 +{\hat{A}}_b h_b^4 |{\hat{\mu }}|^2\Big \}, \end{aligned}$$
(21)
$$\begin{aligned} \kappa&\equiv \frac{1}{16\pi ^2}. \end{aligned}$$
(22)

All hatted parameters above and following in the rest of the paper are normalized to the scale \(M_\text {s}.\)

The triangle diagrams like those of Fig. 1b give

$$\begin{aligned} \Delta \lambda _{1}^{(3)}&= \frac{3}{4 }\kappa \Big \{-|{\hat{A}}_b|^2 h_b^2 \Big (g^2+g_y^2-8h_b^2\Big )\nonumber \\&\quad +\Big (g^2+g_y^2\Big ) h_t^2 |{\hat{\mu }}|^2\Big \}, \end{aligned}$$
(23)
$$\begin{aligned} \Delta \lambda _{2}^{(3)}&= \frac{3}{4 }\kappa \Big \{-|{\hat{A}}_t|^2 h_t^2 \Big (g^2+g_y^2-8h_t^2\Big )\nonumber \\&\quad +\Big (g^2+g_y^2\Big ) h_b^2 |{\hat{\mu }}|^2\Big \}, \end{aligned}$$
(24)
$$\begin{aligned} \Delta \lambda _{3}^{(3)}&= -\frac{3}{8}\kappa \Big \{h_t^2|{\hat{A}}_t|^2 (g^2-g_y^2 -4 h_b^2)\nonumber \\&\quad +h_b^2\ |{\hat{A}}_b|^2 \Big (g^2-g_y^2-4 h_t^2\Big )\nonumber \\&\quad -h_b^2 |{\hat{\mu }}|^2\left[ g^2-g_y^2+4 (h_b^2 -h_t^2)\right] \nonumber \\&\quad -h_t^2 |{\hat{\mu }}|^2\left[ g^2-g_y^2 + 4( h_t^2 - h_b^2)\right] \nonumber \\&\quad - 4 h_b^2 h_t^2({\hat{A}}_b {\hat{A}}_t^* + {\hat{A}}_b^* {\hat{A}}_t) \Big \}, \end{aligned}$$
(25)
$$\begin{aligned} \Delta \lambda _{4}^{(3)}&= \frac{3}{4 }\kappa \Big \{+h_t^2 |{\hat{A}}_t|^2 \left( g^2 - 2 h_b^2\right) \nonumber \\&\quad +h_b^2 |\hat{A_b}|^2 \left( g^2- 2 h_t^2\right) \nonumber \\&\quad - |{\hat{\mu }}|^2(h_t^2+h_b^2)\left[ g^2- 2 (h_b^2 +h_t^2)\right] \nonumber \\&\quad -2 h_b^2 h_t^2 (\hat{A_b} {\hat{A}}_t^* + {\hat{A}}_b^* {\hat{A}}_t) \Big \}, \end{aligned}$$
(26)
$$\begin{aligned} \Delta \lambda _{5}^{(3)}&= 0, \end{aligned}$$
(27)
$$\begin{aligned} \Delta \lambda _{6}^{(3)}&= \frac{3 {\hat{\mu }}}{8 }\kappa \Big \{\hat{A_b} h_b^2 \Big (g^2+g_y^2-8 h_b^2\Big )\nonumber \\&\quad - {\hat{A}}_t h_t^2 \Big (g^2+g_y^2\Big )\Big \}, \end{aligned}$$
(28)
$$\begin{aligned} \Delta \lambda _{7}^{(3)}&= -\frac{3 {\hat{\mu }}}{8 }\kappa \Big \{\hat{A_b} \Big (g^2+g_y^2\Big ) h_b^2 \nonumber \\&\quad -{\hat{A}}_t h_t^2 \Big (g^2+g_y^2-8 h_t^2\Big )\Big \}. \end{aligned}$$
(29)

There are also contributions coming from the redefinition of the Higgs doublets. Squark loops induce mixing between the scalar fields, which must be accounted for in order to preserve canonically normalized kinetic terms for the scalar fields in the Lagrangian. This is done by redefining the Higgs doublet fields in the following manner

$$\begin{aligned} \begin{pmatrix}\Phi _1 \\ \Phi _2 \end{pmatrix} \rightarrow \begin{pmatrix}\Phi _1 \\ \Phi _2 \end{pmatrix} - \frac{1}{2}\begin{pmatrix} \Delta Z_{\Phi _1 \Phi _1} &{}\quad \Delta Z_{\Phi _1 \Phi _2} \\ \Delta Z_{\Phi _2 \Phi _1} &{}\quad \Delta Z_{\Phi _2 \Phi _2} \end{pmatrix}\begin{pmatrix}\Phi _1 \\ \Phi _2 \end{pmatrix}. \end{aligned}$$
(30)

The SU(2) invariance ensures that the corrections can be applied to the complete Higgs doublets and not only the component fields. The expressions for the wave-function-correction factors \(\Delta Z_{\Phi _i \Phi _j}\) can be derived via the finite parts of the derivatives of the self energies in the electroweak interaction basis \(\Sigma '_{\phi _i \phi _j}\) with \(\phi _{\{i, j\}} = \{\phi _1, \phi _2, a_1, a_2\}\), corresponding to a \(\overline{\textrm{MS}}\) renormalized self energy,

$$\begin{aligned} \Delta Z_{\Phi _1 \Phi _1}&= \frac{1}{2}\left( \Sigma '_{\phi _1 \phi _1} + \Sigma '_{a_1 a_1} \right) \nonumber \\&= \kappa \frac{h_t^2 |{\hat{\mu }}|^2 + h_b^2 |\hat{A_b}|^2}{2}, \end{aligned}$$
(31)
$$\begin{aligned} \Delta Z_{\Phi _1 \Phi _2}&= \frac{1}{2}\left( \Sigma '_{\phi _1 \phi _2} + \text {i} \Sigma '_{\phi _1 a_2} - \text {i} \Sigma '_{\phi _2 a_1} + \Sigma '_{a_1 a_2} \right) \nonumber \\&=-\kappa \frac{{\hat{\mu }} (h_t^2 {\hat{A}}_t + h_b^2 \hat{A_b})}{2 }, \end{aligned}$$
(32)
$$\begin{aligned} \Delta Z_{\Phi _2 \Phi _1}&= \Delta Z_{\Phi _1 \Phi _2}^*, \end{aligned}$$
(33)
$$\begin{aligned} \Delta Z_{\Phi _2 \Phi _2}&= \frac{1}{2}\left( \Sigma '_{\phi _2 \phi _2} + \Sigma '_{a_2 a_2} \right) \nonumber \\&= \kappa \frac{h_t^2 |{\hat{A}}_t|^2 + h_b^2 |{\hat{\mu }}|^2}{2}. \end{aligned}$$
(34)

In Ref. [81], it was shown for the CP-even Higgs boson fields that this choice for the wave-function-correction factors together with an appropriate choice of correction of the mixing angle leads to the physical fields being the same in the MSSM and the 2HDM at the matching scale as required. The field redefinitions lead to the following threshold corrections

$$\begin{aligned} \Delta \lambda _{1}^{(2)}&= -\frac{g^2+g_y^2}{4 } \kappa \Big (h_b^2 |\hat{A_b}|^2 + h_t^2 |{\hat{\mu }}|^2\Big ), \end{aligned}$$
(35)
$$\begin{aligned} \Delta \lambda _{2}^{(2)}&= -\frac{g^2+g_y^2}{4 } \kappa \Big (h_t^2 |{\hat{A}}_t|^2 + h_b^2|{\hat{\mu }}|^2\Big )), \end{aligned}$$
(36)
$$\begin{aligned} \Delta \lambda _{3}^{(2)}&= -\frac{g^2-g_y^2}{8 } \kappa \Big ( h_t^2 \left( |{\hat{A}}_t|^2 +|{\hat{\mu }}|^2\right) \nonumber \\&\quad + h_b^2 \left( |\hat{A_b}|^2+|{\hat{\mu }}|^2\right) \Big ), \end{aligned}$$
(37)
$$\begin{aligned} \Delta \lambda _{4}^{(2)}&= \frac{g^2\kappa }{4}\Big (h_b^2(|\hat{A_b}^2| + |{\hat{\mu }}|^2) + h_t^2(|{\hat{A}}_t^2| + |{\hat{\mu }}|^2)\Big ), \end{aligned}$$
(38)
$$\begin{aligned} \Delta \lambda _{5}^{(2)}&= \Delta \lambda _{6}^{(2)} = \Delta \lambda _{7}^{(2)} = 0. \end{aligned}$$
(39)

Finally, the Yukawa couplings receive the one-loop corrections resulting in the following 2HDM Yukawa couplings at the matching scale (including the tree-level contribution):

$$\begin{aligned} h_t^{\text {2HDM}}&= h_t\bigg \{1 - \kappa \bigg [\frac{4}{3}g_s^2\left( {\hat{A}}_t{\hat{M}}^*_3 - 1\right) \nonumber \\&\quad + \frac{1}{4}\left( h_t^2 |{\hat{A}}_t|^2 + h_b^2 |{\hat{\mu }}|^2\right) \nonumber \\ {}&\quad {{- h_b^2 |{\hat{\mu }}|^2 F_3(|{\hat{\mu }}|^2) - \frac{1}{8} (h_b^2 + 3 h_t^2)F_1(|{\hat{\mu }}|^2)}} \bigg ]\bigg \}, \end{aligned}$$
(40)
$$\begin{aligned} h_b^{\text {2HDM}}&= h_b\bigg \{1 - \kappa \bigg [\frac{4}{3}g_s^2\left( \hat{A_b}{\hat{M}}^*_3 - 1\right) \nonumber \\&\quad + \frac{1}{4}\left( h_b^2 |\hat{A_b}|^2 + h_t^2 |{\hat{\mu }}|^2\right) \nonumber \\&\quad {{- h_t^2 |{\hat{\mu }}|^2 F_3(|{\hat{\mu }}|^2) - \frac{1}{8} (h_t^2 + 3 h_b^2)F_1(|{\hat{\mu }}|^2)}} \bigg ]\bigg \}, \end{aligned}$$
(41)
$$\begin{aligned} h_t^{'\text {2HDM}}&= \kappa h_t\left\{ \frac{4}{3}g_s^2{\hat{\mu }}^*{\hat{M}}^*_3 + \frac{1}{4} \left( h_b^2 {\hat{A}}_b^* {\hat{\mu }}^* + h_t^2 {\hat{A}}_t^*{\hat{\mu }}^*\right) \right. \nonumber \\&\quad \left. {{ - h_b^2 {\hat{A}}_b^* {\hat{\mu }}^* F_3(|{\hat{\mu }}|^2)}} \right\} , \end{aligned}$$
(42)
$$\begin{aligned} h_b^{'\text {2HDM}}&= \kappa h_b\left\{ \frac{4}{3}g_s^2{\hat{\mu }}^*{\hat{M}}^*_3 + \frac{1}{4} \left( h_b^2 {\hat{A}}_b^* {\hat{\mu }}^* + h_t^2 {\hat{A}}_t^*{\hat{\mu }}^*\right) \right. \nonumber \\&\quad \left. {{ - h_t^2 {\hat{A}}_t^* {\hat{\mu }}^* F_3(|{\hat{\mu }}|^2)}}\right\} \end{aligned}$$
(43)

with

$$\begin{aligned} F_1(x)&= -\frac{1 - 4 x + x^2 \left[ 3 - 2 \ln (x)\right] }{(1 - x)^2} , \qquad&F_1(1)&= 0, \end{aligned}$$
(44)
$$\begin{aligned} F_3(x)&= -\frac{1 - x \left[ 1 - \ln (x)\right] }{(1 - x)^2},&F_3(1)&= \frac{1}{2}. \end{aligned}$$
(45)

It should be noted that, since the absolute value of the gluino mass parameter is \(|M_3| = M_s\), \({\hat{M}}_3\) is just a phase factor \({\hat{M}}_3= \text {e}^{\text {i} \varphi _{M_3}}\) where \(\varphi _{M_3}\) is the phase of the gluino mass parameter.

We do not calculate threshold corrections to parameters such as \(\tan \beta \), since they do not enter in the MSSM threshold corrections and, hence, are only needed as 2HDM parameters.

5 Calculating the Higgs-mass spectrum and the mixing

The Higgs masses are determined completely once all the parameters of the MSSM at the scale \(M_s\) are given. These are the necessary boundary conditions for solving the RGEs and obtaining the values for the relevant couplings at the scale where the masses are calculated. However, not all the relevant input parameters are given at the same scale. The soft-SUSY breaking parameters of the MSSM are given as user-defined input at the scale \(M_s\), \(\tan \beta \) is given at the scale \(M_{H^+}\) (hence treated as a 2HDM parameter) while all SM couplings relevant to the calculation are fixed at the electroweak scale.

There are different ways to approach this mixed-scale boundary-value issue. The “bottom up” approach starts with the low-energy scale values from the SM and evolves the parameters up to the high-energy scale taking matching effects into account on the way up to the high-energy scale and guessing the values of the first parameters such as \(\lambda _i\). In an iterative procedure, evolving the parameters up and down, the complete set of parameters at a single energy scale is found. With the “top down” approach, which has been exploited also in Refs. [80, 81], one guesses initial values for the high scale MSSM parameters and evolves all the couplings down to \(M_t\). Here, the couplings calculated from the EFT procedure are compared to the experimentally fixed values, and the high scale parameters are adjusted to minimize the differences using a numerical algorithm. This way, evolving parameters up to the high scale can be avoided. We adapt the “top down” approach, employing an upwards evolution of the parameters just for the initial guess. The process is sketched in Fig. 2 and the single steps are described in the following:

  1. 1.

    First, initial values for the high scale MSSM couplings as a first guess have to be found. To obtain these,

    1. (a)

      we start at the scale \(M_t\) (with \(M_t\) being the top pole mass) where the SM couplings are fixed, guess a value for the SM quartic Higgs coupling of \(\lambda _{\text {SM}} = 0.25\), and evolve the SM couplings up to the intermediate scale \(M_{H^+}\) using SM RGEs obtained from SARAH [90, 91].

    2. (b)

      At the scale \(M_{H^+}\), it is assumed that all the 2HDM quartic couplings \(\lambda _1\),..., \(\lambda _7\), “wrong” Yukawa couplings \(h_t'\) and \(h_b'\), and the phases of the Yukawa couplings \(\varphi _{h_t}\), \(\varphi _{h_t'}\), \(\varphi _{h_b}\), and \(\varphi _{h_b'}\) are zero, and the 2HDM Yukawa couplings \(h_t\) and \(h_b\) are calculated accordingly via the tree-level matching of the Yukawa couplings,

      $$\begin{aligned} h_t^{\textrm{2HDM}}= & {} \frac{1}{\sin \beta }\, y_t^{\textrm{SM}}, \end{aligned}$$
      (46)
      $$\begin{aligned} h_b^{\textrm{2HDM}}= & {} \frac{1}{\cos \beta }\, y_b^{\textrm{SM}}. \end{aligned}$$
      (47)

      Then, the 2HDM couplings are evolved up to the scale \(M_s\) using the full two-loop 2HDM RGEs including complex phases. For the gauge, Yukawa, and quartic couplings, we have calculated these implementing the general prescription first developed by Refs. [108,109,110,111] and expanded upon in Refs. [112, 113] to account for kinetic mixing of scalar fields in the presence of multiple Higgs doublets. For the running vevs, we use the formulae from Refs. [114, 115]. We have checked our results for the couplings with the findings of the authors of Refs. [116, 117] and find agreement.Footnote 3

    3. (c)

      As our first guess, the values of the 2HDM gauge and Yukawa couplings emerging from the previous step are taken to determine the initial values of the MSSM gauge and Yukawa couplings,

      $$\begin{aligned} c^{\text {MSSM}} = c^{\text {2HDM}} \quad \text {with} \quad c = g_y, g, g_s, h_t, h_b. \end{aligned}$$
      (48)

      It should be noted that in the MSSM, the “wrong” Yukawa couplings are purely loop-induced and that the Yukawa phases can be absorbed into the fields. Hence, we only have the real parameters \(h_t^\textrm{MSSM}\) and \(h_b^\textrm{MSSM}\).

  2. 2.

    Now, the MSSM parameters are given by the gauge and Yukawa couplings obtained in step 1 (or adapted in the minimization procedure) and the soft-SUSY breaking parameters \(A_t\), \(A_b\), \(\varphi _{M_3}\) as well as the parameter \(\mu \) defined as input at the scale \(M_s\). They are used in the following steps:

    1. (a)

      With the MSSM parameters, the 2HDM couplings are calculated using the matching conditions given in Sect. 4. These MSSM threshold corrections give the non-vanishing values for the 2HDM quartic couplings \(\lambda _1, \dots , \lambda _7\), the “wrong” Yukawa couplings \(h_t^{\prime \textrm{2HDM}}\), \(h_b^{\prime \textrm{2HDM}}\), and the Yukawa phases of the 2HDM \(\varphi _{h_t}\), \(\varphi _{h_t'}\), \(\varphi _{h_b}\), and \(\varphi _{h_b'}\). The couplings are then run down to the scale \(M_{H^+}\).

    2. (b)

      Then, the 2HDM is matched to the SM. The tree-level matching conditions for the SM Yukawa couplings \(y_t\) and \(y_b\) to the 2HDM ones \(h_t\), \(h_b\), \(h_t'\) and \(h_b'\) are

      $$\begin{aligned}&\sqrt{|h_t|^2\sin ^2\beta + 2|h_t||h_t^\prime |\cos \beta \sin \beta \cos (\varphi _{h_t} - \varphi _{h_t^\prime }) + |h_t^\prime |^2\cos ^2\beta } = y_t\nonumber \\&\sqrt{|h_b|^2\cos ^2\beta + 2|h_b||h_b^\prime |\cos \beta \sin \beta \cos (\varphi _{h_b} - \varphi _{h_b^\prime }) + |h_b^\prime |^2\sin ^2\beta } = y_b \end{aligned}$$
      (49)

      where \(\varphi _\chi \) is the phase of the coupling \(\chi \) with \(\chi = h_t, h_t^\prime , h_b, h_b^\prime \). The vacuum expectation value v in the SM, we obtain via the relation

      $$\begin{aligned} v^2 = v_1^2 +v_2^2. \end{aligned}$$
      (50)

      The quartic Higgs coupling in the SM \(\lambda _{\textrm{SM}}\) can be calculated at tree-level via

      $$\begin{aligned} \lambda _{\text {SM}}= & {} c_\beta ^4 \lambda _1 + 4 c_\beta ^3 s_\beta {\text {Re}}(\lambda _6) + 2 c_\beta ^2 s_\beta ^2\nonumber \\{} & {} \times \left[ \lambda _3 + \lambda _4 + {\text {Re}}(\lambda _5)\right] \nonumber \\{} & {} + 4 c_\beta s_\beta ^3 {\text {Re}}(\lambda _7) + s_\beta ^4 \lambda _2. \end{aligned}$$
      (51)

      The one-loop threshold correction to \(\lambda _{\text {SM}}\) is obtained by integrating out the heavy Higgs bosons. In the real case where all phases are set to zero, the answer is known in closed form [81].Footnote 4 In the complex case, on the other hand, the calculation is complicated by the \(4\times 4\) neutral mixing and mass matrices. We have evaluated the full threshold corrections numerically for the complex case, which leads to the problem that the result includes contributions of order \({\mathcal {O}}(v/M_{H^+})\) that are ignored elsewhere in the calculation. Comparing the results for the Higgs masses using only the tree-level matching, the one-loop threshold of Ref. [81], and the full one-loop threshold including \({\mathcal {O}}(v/M_{H^+})\) terms leads to very small differences, so we can neglect the one-loop threshold entirely. Similarly, the one-loop 2HDM threshold corrections to the SM Yukawa couplings are numerically negligible.

    3. (c)

      In the next step, the SM couplings are evolved from the scale \(M_{H^+}\) down to \(M_t\) and checked against the experimental values forFootnote 5\(g_y\), g, \(g_s\), \(y_t\), \(y_b\). We repeat this procedure, adjusting the high scale MSSM couplings each time to minimize the differences between the SM couplings at the low-scale and the experimental values until good agreement is found.

  3. 3.

    Via the minimization procedure in step 2, we obtained a final set of values for all MSSM high scale parameters. These are evolved down one last time to \(M_{H^+}\). At this stage, all the low-scale 2HDM parameters necessary for computing the Higgs masses are determined, and one could in principle calculate the eigenvalues of the loop-corrected mass matrix at the scale \(M_{H^+}\) and determine the pole masses. However, this will lead to terms containing potentially large logarithms of \(\ln (M_{H^+}/m_t)\) with \(m_t\) being the running top quark mass, which are additionally enhanced by factors of the large top-Yukawa coupling. These terms originate from the one-loop corrections in the conversion of the \(\overline{\text {MS}}\) mass to the pole mass of the SM-like Higgs boson. Therefore, we considered three conceptually different methods to calculate the Higgs-boson masses: In the case that the charged Higgs boson is sufficiently light, the 2HDM can be used as the low-energy theory (options (a) and (b) below). If the charged Higgs boson is heavy, then the SM is the appropriate low-energy theory and a matching procedure for the 2HDM and the SM is performed at the scale \(M_{H^+}\) (option (c)). Finally, we apply an approximation that interpolates between the two results (option (d)). In the following, we list the options and include some details about the calculation:

    1. 1.

      The parameters are taken at the scale \(\mu _{\text {ren}} = M_{H^+}\) and the on-shell Higgs masses are calculated via the zeros of the determinant

      $$\begin{aligned} \det \left[ p^2 - {\mathcal {M}}^2 (\mu _{\text {ren}}) + {\hat{\Sigma }}(\mu _{\text {ren}}, p^2) - {\hat{T}}\right] = 0 \end{aligned}$$
      (52)

      expanded up to one-loop order where \({\hat{\Sigma }}(\mu _{\text {ren}}, p^2)\) denotes the top and bottom Yukawa contributions to the self energy matrix in the \(\overline{\text {MS}}\) renormalization scheme at momentum \(p^2\). To ensure the proper minimum of the effective potential, tadpole contributions \({\hat{T}}\) originating from top and bottom loops have to be taken into account. The entries of the matrix \({\hat{T}}\) are given in the Appendix B in terms of tadpole contributions in the interaction basis. The mass matrix \({\mathcal {M}}^2 (\mu _{\text {ren}})\) has the form of the tree-level mass matrix given in Eqs. (10) with the parameters evaluated at the scale \(\mu _{\text {ren}} = M_{H^+}\), where the charged Higgs mass is given in the \(\overline{\text {MS}}\) scheme.

    2. 2.

      In this option, the low-energy theory is still the 2HDM, however, the parameters are evolved down to the running top-quark mass \(m_t\) calculated in terms of 2HDM parameters, and Eq. (52) is evaluated at the scale \(\mu _{\text {ren}} = m_t\) where the \(\overline{\text {MS}}\) mass of the charged Higgs boson \(M_{H^+}\) is interpreted as given at the scaleFootnote 6\(m_t\). Using this scale choice, the logarithms \(\ln (\mu _{\text {ren}}/m_t)\) in the self energies vanish, since we evaluate the self-energies using the running top and bottom masses. This method is only valid as long as \(M_{H^+}\) is not much larger than \(m_t\) since the 2HDM RGEs are not the correct RGEs for evolving the couplings below the scale \(M_{H^+}.\)

    3. 3.

      In this method, the SM is decoupled completely from the 2HDM and treated as the low-energy theory. This method applies when \(M_{H^+}\gg m_t\). In this case, the heavy Higgs bosons of the 2HDM are decoupled by matching the 2HDM to the SM in the same way as step 2b, and the SM couplings are evolved down to \(m_t\). In this case, however, \(m_t\) is calculated in terms of SM parameters only. The \(\overline{\text {MS}}\) mass of the lightest Higgs boson is taken to be \(v^2(m_t) \lambda _{\textrm{SM}}(m_t)\), which is converted to the pole mass via

      $$\begin{aligned} p^2 - v^2(m_t) \lambda _{\textrm{SM}}(m_t) + {\hat{\Sigma }}^\text {SM} - {\hat{T}}^\text {SM} = 0 \end{aligned}$$
      (53)

      where \( {\hat{\Sigma }}^\text {SM}\) and \({\hat{T}}^\text {SM}\) are the self-energy and tadpole contributions of the SM-like Higgs boson of \({\mathcal {O}}(\alpha _t)\) and \({\mathcal {O}}(\alpha _b)\) with \(\alpha _{\{t,b\}} = y_{\{t,b\}}^2/(4\pi )\) evaluated in the \(\overline{\text {MS}}\) scheme. In this option, since the heavy 2HDM Higgs bosons are decoupled, all information about the heavy Higgs bosons at the scale \(m_t\) is encoded in the size of the SM couplings that are obtained via the tree-level matching of the 2HDM couplings to the SM at the scale \(M_{H^+}\), where the former result from the evolution of the 2HDM parameters from the scale \(M_s\) to scale \(M_{H^+}\). The masses of the heavy Higgs bosons can be estimated to be of order \(M_{H^+}\). However, at the scale \(M_{H^+}\), one can still obtain direct information about the heavy Higgs bosons.

    4. 4.

      Finally, in this option, an approximation is exploited that allows one to resum large logarithms in the scenario of \(M_{H^+}\gg m_t\) while still retaining information about the heavy Higgs bosons and the mixing between them and the lightest Higgs boson. In order to do so, firstly, we only consider top-Yukawa effects. We do not attempt to resum logarithms proportional to other couplings in the 2HDM. This approximation is good in the low-\(\tan \beta \) regime where \(h_b^{\textrm{2HDM}}\) is small, which is also the most phenomenologically relevant regime, especially for low values of the mass of the charged Higgs boson. The second assumption is that for the resummation of these logs, our 2HDM can be considered a classic type-II CP-even 2HDM where only one Higgs doublet couples to the right-handed top quarks. This is because the logarithms we wish to resum arise from one-loop corrections, and CP-violating and “wrong-type” Yukawa couplings are already loop-suppressed, so any effect these couplings will be suppressed by an extra loop factor. With this regime in mind, we wish to incorporate the effect of running from \(M_{H^+}\) to \(m_t\) (with \(m_t\) evaluated with parameters of the 2HDM) into the full 2HDM neutral mass matrix at the scale \(M_{H^+}\). To begin, we evaluate \(\lambda _{\text {SM}}\) and v at \(M_{H^+}\) and \(m_t\) according to steps 2b and 2c, but we take \(y_b = 0\) into account when running the SM couplings down to \(m_t\). Then, at the scale \(M_{H^+}\), we rotate into the so-called “Higgs Basis” [118,119,120], defined by

      $$\begin{aligned}{} & {} H_1 = c_{\beta }\Phi _1 + s_{\beta }\Phi _2, \nonumber \\{} & {} H_2 = c_{\beta }\Phi _2 - s_{\beta }\Phi _1,\nonumber \\{} & {} H_1 = \begin{pmatrix} h_1^+\\ \frac{1}{\sqrt{2}}(v+h_1+ib_1)\end{pmatrix},\nonumber \\{} & {} H_2 = \begin{pmatrix} h_2^+\\ \frac{1}{\sqrt{2}}(h_2 +ib_2)\end{pmatrix}\nonumber \\{} & {} s_\beta \equiv \sin \beta , \qquad c_\beta \equiv \cos \beta ,\nonumber \\{} & {} \tan \beta \equiv \frac{v_2}{v_1},\qquad v^2\equiv v_1^2 + v_2^2, \end{aligned}$$
      (54)

      where \( h_j^+\), \(h_j\), and \(b_j\) with \(j =1, 2\) are the charged, the CP-even, the CP-odd Higgs fields in the Higgs basis, respectively. In this basis, only Higgs doublet \(H_1\) gets a vev v and can therefore be identified with the SM Higgs doublet. The mass matrix in the Higgs basis can be obtained via \({\mathcal {M}}^\text {Higgs} = {\mathcal {U}} {\mathcal {M}}^2 {\mathcal {U}}^\dagger \), where \({\mathcal {M}}^2\) is given in Eq. (10) and

      $$\begin{aligned} {\mathcal {U}} = \begin{pmatrix} U &{}\quad 0 \\ 0 &{}\quad U\end{pmatrix} \quad \text {with} \quad U = \begin{pmatrix} c_\beta &{}\quad s_\beta \\ -s_\beta &{}\quad c_\beta \end{pmatrix}, \end{aligned}$$
      (55)

      and the (1,1) component of \({\mathcal {M}}^\text {Higgs}\) can be identified with \(v_{\textrm{SM}}^2\lambda _{\textrm{SM}}\) with \(v_{\textrm{SM}} = v\) leading to the threshold condition given in Eq. (51). The submatrix of \({\mathcal {M}}^\text {Higgs}\) given by the second and third row and column describe the heavy Higgs bosons. Before continuing, we note two relevant facts. First, in the decoupling limit [119]Footnote 7 where \(M_{H^+}\gg v\), the mixing angle \(\alpha \) (which diagonalizes the CP-even neutral mass matrix in the CP-conserving 2HDM) can be approximated by \(\beta -\frac{\pi }{2},\) and the Higgs basis is up to a minus sign the mass basis,

      $$\begin{aligned}{} & {} \begin{pmatrix} h\\ H\end{pmatrix} = \begin{pmatrix} -s_{\alpha } &{}\quad c_{\alpha } \\ c_{\alpha } &{}\quad s_{\alpha }\end{pmatrix} \begin{pmatrix} \phi _1 \\ \phi _2 \end{pmatrix} \xrightarrow {\tiny M_{H^+}\gg v} \begin{pmatrix} c_{\beta } &{}\quad s_{\beta }\\ s_{\beta } &{}\quad - c_{\beta } \end{pmatrix}\nonumber \\{} & {} \quad \begin{pmatrix}\phi _1\\ \phi _2\end{pmatrix} \end{aligned}$$
      (56)

      which according to Eq. (54) shows that

      $$\begin{aligned} h=h_1\qquad H=-h_2. \end{aligned}$$
      (57)

      Second, since at tree-level in the gauge-eigenstate basis only the 2HDM doublet \(\Phi _2\) couples to right-handed top quarks according to our above assumptions, one can approximate the self-energy corrections to the CP-even Higgs bosons by

      $$\begin{aligned} \Sigma ^{\textrm{tops}}_{h_1h_1} = \Sigma ^{\textrm{tops}}_{hh}&= s_\beta ^2\Sigma ^{\textrm{tops}}_{\phi _2\phi _2}\nonumber \\ \Sigma ^{\textrm{tops}}_{h_1h_2} = -\Sigma ^{\textrm{tops}}_{hH}&= c_\beta s_\beta \Sigma ^{\textrm{tops}}_{\phi _2\phi _2} = \frac{1}{\tan \beta }\Sigma ^{\textrm{tops}}_{hh} \nonumber \\ \Sigma ^{\textrm{tops}}_{h_2h_2} = \Sigma ^{\textrm{tops}}_{HH}&= c_\beta ^2 \Sigma ^{\textrm{tops}}_{\phi _2\phi _2} = \frac{1}{\tan ^2\beta }\Sigma ^{\textrm{tops}}_{hh}. \end{aligned}$$
      (58)

      Now we consider \(v^2(m_t)\lambda _{\textrm{SM}}(m_t)\) to be the leading-log resummation of the one-loop-leading-log contribution coming from the \(\Sigma ^{\textrm{tops}}_{h_1h_1}\) self-energy correction. Therefore, using Eq. (58), we incorporate this into the Higgs basis matrix as follows

      $$\begin{aligned}{} & {} {\mathcal {M}}^{\textrm{Higgs}}_{\text {approx}} = \mathcal M^{\textrm{Higgs}} \nonumber \\ {}{} & {} \quad + \begin{pmatrix}\gamma ^{\textrm{resum}} &{}\quad \frac{1}{\tan \beta }\gamma ^{\textrm{resum}} &{}\quad 0 &{}\quad 0\\ \frac{1}{\tan \beta }\gamma ^{\textrm{resum}} &{}\quad \frac{1}{\tan ^2\beta }\gamma ^{\textrm{resum}}&{}\quad 0 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad 0&{}\quad 0 \\ 0 &{}\quad 0&{}\quad 0 &{}\quad 0\end{pmatrix}, \end{aligned}$$
      (59)

      where \(\gamma ^{\textrm{resum}}\equiv v^2(m_t)\lambda _{\textrm{SM}}(m_t) - v^2(M_{H^+})\lambda _{\textrm{SM}}(M_{H^+}).\) Then, the eigenvalues of the matrix

      $$\begin{aligned} {\mathcal {M}}^{\textrm{Higgs}}_{\textrm{approx}} - {\hat{\Sigma }}(p^2,\mu ) + {\hat{T}} \end{aligned}$$
      (60)

      are calculated, where \({\hat{\Sigma }}\) is the \(4\times 4\) matrix of one-loop self-energies of the neutral Higgs bosons in the gauge basis and \({\hat{T}}\) is the matrix of one-loop neutral tadpole corrections in gauge basis given in Appendix B, both in the \(\overline{\textrm{MS}}\) scheme and using the parameters at \(m_t\). To calculate these eigenvalues, it is necessary to choose the renormalization scale \(\mu \) and the external momentum \(p^2\). The renormalization scale is set to \(m_t\), and the external momenta are chosen to be the tree-level masses. We choose to calculate the eigenvalues one a time, with the mass matrix diagonalized once for each tree-level mass. For example, the squared mass of the lightest neutral Higgs boson corresponds to the lightest eigenvalue of the loop-corrected mass matrix evaluated with the external momentum set to the tree-level neutral Higgs mass. In order to proceed along the same lines as in the pure 2HDM case, see step 3a, the matrix \(\mathcal M^{\textrm{Higgs}}_{\text {approx}}\) can be rotated back to the interaction eigenstates using the transformation matrix (55). This results in \({\mathcal {M}}^2\) with an additional contribution to the (2,2) element,

      $$\begin{aligned} \mathcal (M^2_{22})^{\textrm{resum}} = {\mathcal {M}}^2_{22} + \frac{1}{\sin ^2 \beta } \gamma ^{\textrm{resum}}. \end{aligned}$$
      (61)

      Using this new matrix \(({\mathcal {M}}^2)^{\textrm{resum}}\) and replacing \({\mathcal {M}}^2\) by \(({\mathcal {M}}^2)^{\textrm{resum}}\) in Eq. (52), one can calculate the Higgs masses by finding the zeros of the resulting equation up to one-loop order. We find that both approaches lead to nearly the same result.

    For options (b) to (d), an evaluation at the scale of the running top-quark mass \(m_t (m_t)\) is performed. In these cases, the running top-quark mass is calculated iteratively.

Fig. 2
figure 2

Pictorial description of the mass calculation

6 Numerical results

6.1 Choice of input parameters

The list of relevant input parameters for the calculation is

$$\begin{aligned}{} & {} \{\underbrace{y_t,\;y_b,\;g_y,\;g_3,\;g_2,\;v,}_{\text {at } M_t}\;\underbrace{\tan \beta ,\;M_{H^+},}_{\text {at } M_{H^+}} \nonumber \\{} & {} \quad \underbrace{A_t,\;A_b,\;\varphi _{M_3},\;\mu ,\;M_s}_{\text {at } M_s}\}. \end{aligned}$$
(62)

As mentioned in Sect. 5, these parameters are fixed at different scales. The soft-SUSY breaking parameters of the MSSM and the Higgs mixing parameter \(\mu \) are defined at the scale \(M_s\), \(\tan \beta \) at the scale \(M_{H^+}\), and the SM input parameters are fixed at the low-scale \(M_t\) by current experimental results. The relevant SM observables needed to define the SM couplings, taken from [71, 121], are

$$\begin{aligned} \alpha _s(M_z)= & {} 0.1184,\qquad M_t = 173.34\;\textrm{GeV}, \nonumber \\ M_W= & {} 80.384\;\textrm{GeV},\qquad M_Z = 91.1876 \;\textrm{GeV}, \nonumber \\ m_b(m_b)= & {} 4.18\;\textrm{GeV},\nonumber \\ v^2_{G_F}\equiv & {} \frac{1}{\sqrt{2}G_F} = 246.21971 \;\textrm{GeV}, \end{aligned}$$
(63)

with \(G_F\) being the Fermi constant. The values for \(y_b,\;y_t,\;g_y,\;g_2,\;g_3\) are extracted from these observables in [71, 121] and given below as running parameters at the scale \(M_t\)

$$\begin{aligned} g_3= & {} 1.1666,\qquad g_2 = 0.64779, \; g_y = \sqrt{\frac{3}{5}} g_1 = 0.35830 \nonumber \\ y_t= & {} 0.94018,\qquad y_b = 0.0156 \end{aligned}$$
(64)

where in the conversion the value of the SM Higgs pole mass \(M_h^{\text {SM}} = 125.15\) GeV was used according to Ref. [121]. One must also determine the running vev \(v_{\,\overline{\textrm{MS}}} \equiv v\) from the vev \(v_{G_F}\), which is experimentally determined via the Fermi constant, which is measured via the muon lifetime. We derive the \(\overline{\textrm{MS}}\) vev from the on-shell vev using

$$\begin{aligned} v^2_{\,\overline{\textrm{MS}}} = v^2_{\text {OS}} + \delta {v}^2_{\text {OS-finite}}, \end{aligned}$$
(65)

defining the on-shell vev by

$$\begin{aligned} v^2_{\text {OS}} \equiv \frac{4 M^2_W s^2_W}{e^2} = v^2_{G_F}(1+\Delta r) \end{aligned}$$
(66)

with the counterterm \(\delta {v}^2_{\text {OS-finite}}\) given in Eq. (85) in the appendix. Here, \(s^2_W = 1-\frac{M_W^2}{M_Z^2}\) denotes the sine squared of the weak mixing angle and \(\Delta r\) parameterizes the one-loop radiative corrections to the muon decay in the Fermi Model [122], the process by which \(v_{G_F}\) is defined. The different formulas needed for the conversion are collected in the Appendix C. The value for \(v_{\,\overline{\textrm{MS}}}\) at the scale of the top pole mass \(M_t\) is then

$$\begin{aligned} v_{\,\overline{\textrm{MS}}}(M_t) = 247.3897 \text { GeV} \end{aligned}$$
(67)

employing again the Higgs-boson mass \(M_h^{\text {SM}} = 125.15\) GeV. In the numerical evaluation of the masses of the Higgs bosons and mixings, we use the numbers given in Eqs. (64) and (67) as input values for the SM parameters.Footnote 8

For the remaining input parameters \(M_{H^+}\), \(\tan \beta \), \(A_t\), \(A_b\), \(M_s\), \(\mu \), no measured values exist, but the experimental searches and measurements constrain the viable parameter space. The exclusion bounds from searches of further Higgs bosons [123, 124] constrain in particular the region of high \(\tan \beta \) and light “heavy” Higgs bosons. Taking these results together with further studies of different parameter scenarios [125, 126], it is clear that scenarios with \(\tan \beta >10\) and \(M_{H^+},M_A <500\) GeV are strongly disfavoured by LHC data. Flavour observables support these constraints, as discussed in Ref. [127] for different types of the 2HDM. In our numerical analysis, we therefore favour values for the mass of the charged Higgs boson of \(M_{H^+} \ge 500\) GeV and \(\tan \beta = 5\). To show specific features of the results of our calculation, we will however partly take into account scenarios that do not fulfill these constraints.

It should be noted that, in the MSSM, some of the non-vanishing parameter phases can be eliminated by symmetry transformations. Hence, only certain combinations of phases are physical, i.e. can change the value of a physical observable. Important constraints of these phase combinations come from electric-dipole moment (EDM) measurements (see e.g. Refs. [128,129,130] for more recent studies of the constraints of the MSSM phases due to EDM). Since in our calculation the phases of the U(1) and SU(2) gaugino-mass parameters \(M_1\) and \(M_2\) do not enter, we can assume that they are chosen such that the effect on the EDM is minimized. Furthermore, larger masses of the SUSY particles tend to relax the constraints coming from the EDM.

It should be noted that, in this paper, we refrain from explicit checks whether a certain parameter point is viable, since we focus on specific features of the results. In particular, using our results at the low scale to study the constraints from the EDMs for the MSSM phases at the high scale will be interesting but is left for future work.

Our default scenario is

$$\begin{aligned} A_t&= A_b \equiv A; \quad |A| = |\mu | = 3M_s; \quad M_{H^+} = 500 \text { GeV}; \nonumber \\ \tan \beta&= 5; \nonumber \\ \varphi _{A_t}&= \varphi _{A_b} = 2.1 \approx 2 \pi /3 = 120^\circ ; \quad \varphi _\mu = 0; \quad \varphi _{M_3} = 0. \end{aligned}$$
(68)

where \(M_s\) is varied. The choice \(|A| = |\mu | = 3M_s\) leads to large threshold corrections to the 2HDM couplings and maximizes the amount of CP-violation introduced into the theory. This way we can give an estimate of the largest effects that can occur. Similarly, we observe that \(\varphi _{A_t} = \varphi _{A_b} = 2.1 \approx 2 \pi /3 = 120^\circ \) maximizes roughly the size of the CP-odd component of the lightest neutral Higgs boson, see Sect. 6.5. We will however deviate from this default scenario in order to study the different characteristics and state that explicitly.

6.2 Influence of the running of complex parameters

Fig. 3
figure 3

The quartic couplings’ dependence on \(M_s\) for the default scenario: \(\tan \beta =5\), \(\varphi _{A}=2.1 \approx 120^\circ \), \(\varphi _{\mu } = \varphi _{M_3}=0\), \(|\mu |=|A|=3M_s\), \(M_{H^+} = 500\) GeV. The red curves are the result of employing complex RGEs, and the blue real RGEs. Dashed lines are the couplings at the high scale \(M_s\), and solid lines are couplings at the scale \(M_{H^+}\)

In this work, we exploit the two-loop RGEs for the 2HDM with each Higgs doublet coupling to right-handed up- as well as right-handed down-type fermions including all phases, see Sect. 5 for the details of our calculation of the RGEs. The first numerical results exemplify the effect of including this phase dependence versus the “real RGE” approximation where the phase dependence is taken into account only via the threshold effects and the phases are assumed to be unaffected by the running.

In Fig. 3, we show the effect of including these phases on the running of the 2HDM quartic couplings. The values of the couplings are plotted against the scale \(M_s\), demonstrating how this dependence changes when the running of both the real and the imaginary part of the couplings is taken into account. The coupling values shown are the ones that either enter the final evolution at the scale \(M_s\) or are calculated via the final evolution to the scale of the charged Higgs mass in step 3. The red lines represent the values obtained with RGEs taking phases into account while the blue ones represent values obtained using only real parameters, determining the sign via the argument of the corresponding parameter at the scale \(M_s\). The dashed lines are the parameters at the scale \(M_s\), while the solid lines are those at \(M_{H^+}\). One can clearly see a dependence on whether the running of the phases is taken into account or not. It changes the resulting absolute value of the couplings \(\lambda _5\), \(\lambda _6\) and \(\lambda _7\) as well as the phases themselves. The dependence of the running on the phases is relatively small; only the phase of \(\lambda _7\) shows a change of up to a couple of degrees. Hence, the overall dependence of the phases determined using the RGEs for the complex case on the value of \(M_s\) is also relatively small. Comparing the phases at \(M_s\) (dashed lines) with the ones obtained with using the real RGEs (blue), one can double-check that indeed the phases do not change when exploiting the real RGEs. The absolute values of \(\lambda _5\) and \(\lambda _6\) change more when the complex RGEs are applied compared to a result with only real RGEs. The opposite is true for the absolute value for \(\lambda _7\), which changes less if the complex RGEs are applied.

Since the phase values at the low scale enter the prediction of the EDM, they can be relevant for checking the exclusion of high-scale CP-violating MSSM scenarios due to the measurements of the EDM.

6.3 Comparison of methods for computing masses

In Sect. 5, the calculation procedure was explained, and in step 3 of this procedure we discussed several possibilities for calculating the mass of the Higgs boson \(m_h\) at the low-energy scale: a) exploiting the 2HDM at the scale \(M_{H^+}\), b) employing the 2HDM at the scale \(m_{t}\), c) matching the 2HDM to the SM and using the SM at the scale \(m_t\), or d) approximating the effects of the matching to the SM and running down to scale \(m_t\). Using the first two approaches, one keeps information about all the Higgs masses and their mixings with one another. This is very important, as we are also interested in the size of the CP-odd component of the lightest Higgs boson. However, if the scale \(M_{H^+} \gg v \sim m_t\), we again encounter large logarithms. The masses of the heavy Higgs bosons will not pose a problem, since the logarithms are not equally enhanced by large prefactors as the ones appearing in the calculation of light Higgs mass and since the relative shifts due to the large tree-level masses are smaller, so we can trust the perturbative results for these masses without an additional resummation of logarithms. These large logarithms are more important for the lightest Higgs, on the other hand.

Fig. 4
figure 4

The mass of the light Higgs boson is shown in dependence on the scale \(M_{H^+}\) for the scenario \(|A|=|\mu |=3M_s\), \(M_s = 3\) TeV, \(\varphi _A = \varphi _\mu = \varphi _{M_3} = 0\), \(\tan \beta = 5, 10, 20\) employing different approximations. The curves labeled “resummed” present the result where the logarithms of the form \(\ln (M_{H^+}/m_t)\) have been resummed by matching to SM, but the Higgs mass is still calculated using the mass matrix of the 2HDM. Those labeled with “SM” are calculated by decoupling completely from the 2HDM, and calculating a purely SM mass. Those labeled with “2HDM” refer to those calculated treating the 2HDM as the low-energy theory. For the “tree-level” curves, one-loop corrections have not been included in the calculation of the pole mass

Figure 4 shows the result of calculating the lightest Higgs-boson mass using the different approaches for real parameters. The calculation of the mass of the lightest Higgs boson using the 2HDM at the scale \(M_{H^+}\) and \(m_t\) is denoted with “\(m_h(M_{H^+})\) 2HDM (a)” and “\(m_h(m_t)\) 2HDM (b)”, respectively. The result where the SM is the low-energy theory is called “\(m_h\) SM pole (c)”. A variant of this result is “\(m_h\) SM pole tree level” where the self-energy contribution \({\hat{\Sigma }}^\text {SM}\) as well as the one-loop tadpole contribution \({\hat{T}}^{\text {SM}}\) in Eq. (53) are neglected. Finally, the results labeled “\(m_h\) resummed (d)” correspond to the approximation where parts of the logarithms \(\ln \left( M_{H^+}/m_t\right) \) are resummed. The result “\(m_h\) resummed tree level” neglects the self-energy and the tadpole contributions in Eq. (60).

The “\(m_h\) SM pole” (c) is expected to be the most “correct” answer for large \(M_{H^+}\), while the result obtained by using the 2HDM as the EFT should be the best result when \(M_{H^+}\sim v \sim m_t\). The two 2HDM results agree for low \(M_{H^+}\) but start to deviate quickly with rising \(M_{H^+}\). In the self-energy contribution in Eq. (52), logarithms of \(\ln \left( M_{H^+}/m_t\right) \) arise if the self energy is evaluated at the scale \(M_{H^+}\). Enhanced via large top Yukawa couplings, this result differs quickly from the other 2HDM result where these logarithms vanish in the self energy due to the scale choice and is taken into account via the running of the parameters. Therefore, it is preferable to calculate the Higgs mass \(m_h\) at the scale \(m_t\). If \(M_{H^+}\) however is large, \(M_{H^+} \gg m_t\), the heavy Higgs-bosons have to be decoupled from the running of the parameters and the SM is the correct low-energy theory. For \(M_{H^+} = 1000\) GeV, the deviation of the “\(m_h\) SM pole (c)” result from the “\(m_h(m_t)\) 2HDM (b)” result is roughly 500 MeV for all values of \(\tan \beta \). That means that, for \(M_{H^+} < 1000\) GeV, the 2HDM can still be used as a reasonable low-energy theory but with increasing theoretical uncertainty for increasing values of \(M_{H^+}\). The difference between both results “\(m_h\) SM pole (c)” and “\(m_h(m_t)\) 2HDM (b)” first decreases and then starts growing. For \(\tan \beta = 10\) and \(\tan \beta = 20\), the increase sets already in for lower values of \(M_{H^+}\).

Fig. 5
figure 5

The upper row presents the mass of the lightest Higgs boson depending on the phase \(\varphi _A\) while the lower row shows the dependence on the absolute value of A for a scenario with \(M_{H^+}=500\) GeV, \(M_s = 5\) TeV, \(\varphi _{M_3}=0\), and \(\tan \beta = 5, 10, 15, 20\)

For small \(M_{H^+}\), the mixing of the Higgs bosons becomes relevant, which leads to a decrease of the Higgs-boson mass in the 2HDM results. The “\(m_h\) SM pole (c)” result does not take the mixing of the Higgs bosons into account due to the decoupling of the heavy Higgs bosons. The “\(m_h\) resummed (d)” result follows nicely the 2HDM result for low values and the SM result for large values of \(M_{H^+}\). Hence, it interpolates well between the two options. Therefore, we have chosen this as default option.

In addition, the effect of the one-loop self energies in the calculation of the pole masses can be read off when comparing the “\(m_h\) resummed (d)” and “\(m_h\) resummed tree level” result (or similar the corresponding SM pole mass results). The top and bottom quark contributions that we take into account lead to approximately a 0.6 GeV rise of the result.

6.4 The mass of the light Higgs boson

In this section, we discuss the dependence of the mass of the lightest Higgs boson on the MSSM input parameters \(\mu \), \(A_t=A_b=A\), \(\tan \beta ,\) and \(\varphi _{M_3}\).

The first row of Fig. 5 shows the dependence on the common phase \(\varphi _A\). Firstly, the plots demonstrate that the sensitivity of the mass of the lightest Higgs boson to the phase is highly dependent on \(\tan \beta \), as the mass fluctuates more with \(\varphi _A\) for low values of \(\tan \beta \). For \(\tan \beta = 5\), the variation of \(\varphi _A\) leads to a change of the Higgs-boson mass of up to almost 20 GeV in the quite extreme case of \(|A|/M_s =3\), while for \(\tan \beta = 20\) it is only about 5 GeV for otherwise the same parameters. Secondly, it can be seen that the quantitative features are strongly dependent on the ratio \(r_{\textrm{inputs}}\equiv \frac{A_b,A_t,\mu }{M_s}\). The qualitative dependence on the phase is in fact opposite for the two cases \(r_{\textrm{inputs}} =2\) and \(r_{\textrm{inputs}} =3\), where the mass has either a minimum for \(r_{\textrm{inputs}} =3\) or near maximum value for \(\varphi _A = 180^\circ \) for \(r_{\textrm{inputs}} = 2\). The different dependence on the phase \(\varphi _A\) with respect to the ratio \(|A|/M_s\) can also be read off the second row of Fig. 5 where the lightest Higgs mass is plotted against the magnitude of the common trilinear |A| with \(\mu \) fixed to \(\mu = M_s\). Here, it is seen that the mass peaks for \(|A| \approx 12.5 \text { TeV } = 2.5 M_s \) for \(\varphi _A = 0\) (left plot), and increasing \(\tan \beta \) shifts the peak slightly to lower values of \(|A|/M_s\). For \(\varphi _A = 2.1 \approx 120^\circ \), the peak is shifted to a lower value \(|A|/M_s\approx 2.2\) while increasing \(\tan \beta \) leads to a shift to slightly higher values in this case. Comparing the left and right plot of the lower row of Fig.  5, one can read off for which \(|A|/M_s\), the values of the Higgs-boson mass are smaller for \(\varphi _A = 0\) than for \(\varphi _A = 180^\circ \) and vice versa resulting in the changed maxima in the upper row of Fig. 5.

The effect of changing \(|\mu |\) can be seen by comparing the upper rows from Fig. 5 with the lower rows. The value of \(\mu \) changes from \(\mu = 2 M_s\) (\(\mu = 3 M_s\)) to \(\mu =M_s\) from the top left (right) figure to the bottom row. Comparing values for \(m_h\) for \(\varphi _A = 0^\circ \) in the upper left (right) plot with the ones for \(|A| = 10\) TeV (\(|A| = 15\) TeV) in the lower left plot leads to a rise of the Higgs-boson mass by roughly 1 GeV. For \(\varphi _A = 2.1 \approx 120^\circ \), the same change leads to a shift of up to 6 GeV in the example scenarios, which is reached for \(\tan \beta = 5\).

Fig. 6
figure 6

Dependence of the mass of the lightest Higgs boson on the phase \(\varphi _{M_3}\) for a scenario with \(M_{H^+}=500\) GeV, \(\varphi _{A} = 0\), \(M_s = 5\) TeV, and \(|A|=|\mu |=3M_s\)

Fig. 7
figure 7

The other parameters of the scenario are \(\tan \beta =5\), \(M_s=30\) TeV, \(\varphi _\mu = \varphi _{M_3} = 0\), and \(|A| = |\mu | = 3M_s\)

Figure 6 shows the dependence on \(\varphi _{M_3}\). The overall dependence on \(\varphi _{M_3}\) is not very large and varying \(\varphi _{M_3}\) leads to changes of the order of 1 GeV. The largest shift of 1.3 GeV in the considered scenarios can be found for \(\tan \beta = 5\) and similar for \(\tan \beta = 20\). For \(\tan \beta = 10, 15\), the changes are slightly smaller. It should be noted that in our approximation the phase \(\varphi _{M_3}\) only enters via the one-loop threshold contributions to the Yukawa couplings. The Yukawa couplings, however, appear only at the one-loop level in the Higgs-mass calculation. Further contributions come from two-loop threshold contributions to the quartic couplings and enter the Higgs-mass calculation already at the tree level. We, however, do not include these, since in a full next-to-next-to leading log (NNLL) resummation, not only the two-loop threshold contributions but also three-loop RGE for the 2HDM would be required, which would change the dependence of the SM-like Higgs mass on the phase \(\varphi _{M_3}\).

6.5 A CP-odd admixture to the light Higgs boson

The detection of a CP-odd component in the SM-like Higgs boson would be a sign of new physics, and it is therefore interesting to explore the size of such a component generated by the CP-violating phases in the complex MSSM. While we leave a detailed analysis of the current experimental sensitivity to this CP-odd component for future work, in this section we aim to give a qualitative picture of the size of this component and its dependence on the relevant parameters.

We calculate the CP-odd component using the tree-level mixing matrix \({\mathcal {M}}^2\) of  (10), inserting the parameter values at \(M_{H^+}\) obtained after the final iteration. First, we separate out the Goldstone field by rotating by \(D{\mathcal {M}}^2D^{\textrm{T}}\), where

$$\begin{aligned} D= \begin{pmatrix} I_2 &{}\quad 0 \\ 0 &{}\quad \begin{matrix} -\sin \beta &{}\quad \cos \beta \\ \cos \beta &{}\quad \sin \beta \end{matrix}\end{pmatrix} \end{aligned}$$
(69)

with \(I_2\) being the two-dimensional identity matrix. This rotates the mass matrix into the basis of the fields \(\phi _1, \phi _2, a, G\), where \(\phi _1,\phi _2\) are pure CP-even fields, \(a\equiv -a_1\sin \beta + a_2\cos \beta \) is pure CP-odd, and G is the the Goldstone boson.

In this basis, the mass matrix is block-diagonal, with a \(3\times 3\) block for the fields \(\phi _1,\phi _2,a\), and a \(1\times 1\) block for the Goldstone boson that does not mix with any of the other fields. The rotation matrix P that diagonalizes the \(3\times 3\) matrix relates the physical fields, h, \(H_2\), \(H_3\) to the interaction fields after electroweak symmetry breaking, \(\phi _1\), \(\phi _2\), a, by

$$\begin{aligned} \begin{pmatrix} h \\ H_2 \\ H_3 \end{pmatrix} = P\begin{pmatrix} \phi _1 \\ \phi _2 \\ a \end{pmatrix}. \end{aligned}$$
(70)

The third column of P gives the size of the CP-odd component of each of the physical fields. We choose to plot the CP-odd percentage, which is just \(P_{i3}^2*100\) for \(i = 1\ldots 3\), rather than the mixing-matrix component itself.

In Fig. 7, we show how the CP-odd component depends on \(\varphi _A\) and \(M_{H^+}\). The solid lines in Fig. 7a, b are the same and depict the result of the procedure just described. It is seen that this component is maximized for values of the phases from \(110^\circ \) to \(120^\circ \), justifying our claim from Sect. 6.1. In addition, the results confirm that the CP-odd component very quickly drops to vanishingly small values as \(M_{H^+}\) increases. This agrees with previous findings, see e.g. Ref. [98]. Since data disfavours a charged Higgs boson with small mass, it also disfavours a sizeable CP-odd component of lightest Higgs boson in the MSSM.

Figure 7a also shows the difference between the one-loop corrected mixing matrix evaluated at zero external momentum (\(p^2=0\)) using the mass matrix given in Eq. (60) and the tree-level case, both with parameters at the scale \(M_{H^+}\), which is evidently numerically small. The largest effect can be seen for \(M_{H^+} = 200\) GeV in the peak region where the inclusion of one-loop corrections leads to an increase of the squared CP-odd component of about 0.2 percentage points. If not otherwise stated, we use the tree-level mixing in the following.

Fig. 8
figure 8

The dependence of the size of the CP-odd component on \(M_{H^+}\) and \(\tan \beta \) is shown for a scenario \(M_s=30\)  TeV, \(\varphi _{M_3} = 0,\) \(\varphi _{A} = 2.1 \approx 120^\circ \), and \(|A| = |\mu | = 3M_s\). The color coding gives the percentage of the lightest Higgs boson that is CP-odd

In Fig. 7b, we compare the approach described above (solid line) with the approach where the mass matrix \(({\mathcal {M}}^2)^{\text {resum}}\) from Eq. (61) is used to calculated the CP-odd admixture (dashed line). The maximal difference is again about 0.2 percentage points and decreases with larger values of \(M_{H^+}\). Not shown are one-loop effects in the resummed approach as well as results when parameters are evaluated at the scale \(m_t\). The corresponding results would lie approximately between the dashed lines in Fig. 7a, b.

Figure 8 shows the percentage of the lightest Higgs boson that is CP-odd (shaded contours) and the mass of lightest Higgs boson (dashed line contours) in the \((M_{H^+},\tan \beta )\) plane. As before, the CP-odd component is calculated according to the approach described above using \(\mathcal M^2\) and parameters evaluated at the scale \(M_{H^+}\) without including loop effects. Again, the CP-odd component drops rapidly with increasing \(M_{H^+}\), and one sees that it is only weakly dependent on \(\tan \beta \). In light of the discussion in Sect. 6.1 and the mass of the observed Higgs boson being \({\sim }\,125\) GeV, the results of this section suggest that a large CP-odd component of the lightest Higgs boson is strongly disfavoured by experimental data. Already for \(M_{H^+} = 260\) GeV in this scenario the CP-odd component drops below 0.5%.

According to Ref. [131], an angle \(\varphi _\tau \) of about \(4^\circ \) can be reached with the high-luminosity run of the LHC. The angle \(\varphi _\tau \) is defined via the effective Lagrangian given in four-component Dirac notation

$$\begin{aligned} {\mathcal {L}}^{\tau -\text {Yukawa}}_{\text {eff}} = - \frac{m_\tau }{v} \kappa _\tau \left( \cos \varphi _\tau {\bar{\tau }} \tau +\sin \varphi _\tau {\bar{\tau }} \text {i} \gamma _5 \tau \right) h \end{aligned}$$
(71)

where \(\kappa _\tau \) denotes the change of the absolute strength of Yukawa coupling with respect to the SM while \(\varphi _\tau \) governs the amount of CP-violation. To connect to the notation of (9), we write this in two-component notation as

$$\begin{aligned} {\mathcal {L}}^{\tau -\text {Yukawa}}_{\text {eff}} = -\frac{m_\tau }{v}\kappa _\tau \left( \cos \varphi _\tau \tau \tau _c - \text {i}\sin \varphi _\tau \tau \tau _c \right) h + \mathrm {h.c.}\nonumber \\ \end{aligned}$$
(72)

In the case of the 2HDM discussed here, assuming that the coupling of the second Higgs doublet to the tau leptons is negligible, \(\tan \varphi _\tau = -s_\beta P_{13}/P_{11}\approx - P_{13}/P_{11}\) for scenarios with \(t_\beta \gtrapprox 5\). A CP-odd component of \(P_{13}^2\star 100 = 0.5\%\) as discussed above will lead to \(\varphi _\tau \gtrapprox 4^\circ \), since \(P_{11} \le 1\). Going back to the scenarios discussed in Fig. 7, the maximal value of \(\varphi _\tau \) is larger than \(20^\circ \) for \(M_{H^+} = 200\) GeV and larger than \(3^\circ \) but below \(4^\circ \) for \(M_{H^+} = 500\) GeV. Hence, these values suggest that if the MSSM is the final answer, then while it will certainly be difficult, it might be possible to determine a CP-odd component experimentally with further improvements of the measurements. However, we neither checked whether the Higgs signal rates are in the experimentally allowed regime – this is particularly relevant if \(\varphi _\tau \) is enhanced by a smaller \(P_{11}\) as is actually the case for the scenario in Fig. 7 for \(M_{H^+} = 500\) GeV – nor took into account any constraints from the electric dipole moments. Thus, to find out whether there is still a viable region of parameter space with a measurable CP-odd admixture of the SM-like Higgs boson in the MSSM with heavy superpartners, further investigations are required.

Fig. 9
figure 9

The dependence of the mass of the second lightest and the heaviest Higgs boson on the phase \(\varphi _A\) is shown for the parameter scenario \(M_s= 30\) TeV, \(M_{H^+}=500\) GeV, \(|A|=|\mu | = 3M_s\), \(\varphi _{M_3} = \varphi _\mu = 0\). The solid lines correspond to the mass of the second lightest Higgs boson, \(m_{H_2}\), and the dashed lines to the mass of the heaviest Higgs boson, \(m_{H_3}\)

Fig. 10
figure 10

The dependence of the CP-odd component in percent for the next-to lightest Higgs boson, \(P_{23}^2\star 100\), in the left plot and of the heaviest Higgs boson, \(P_{33}^2\star 100\), in the right plot on the phase \(\varphi _A\) for the scenario \(\tan \beta =5\), \(M_s=30\) TeV, \(\varphi _\mu = \varphi _{M_3} = 0\), and \(|A| = |\mu | = 3M_s.\) The solid curves are calculated using the tree-level mixing matrix, and the dashed curves are calculated using the one-loop mixing matrix at zero-external momentum (\(p^2 = 0\))

6.6 Masses and mixings of the heavy Higgs bosons

In this section, we explore the masses and mixings of the heavy Higgs bosons \(H_2\) and \(H_3\). The masses are calculated according to Sect. 5 where, in step 3 of the calculation procedure, option (d) is exploited where large logarithms \(\ln (M_{H^+}/m_t)\) are resummed. In Fig. 9, the dependence of the masses of the two heavy neutral Higgs bosons \(m_{H_2}\) and \(m_{H_3}\) on \(\varphi _A\) is shown with solid and dashed lines, respectively. Similar to the mass of the light Higgs boson, the phase dependence is larger for low \(\tan \beta \), with a change of the mass of the second-to lightest neutral Higgs boson of roughly 3.5 GeV and of the heaviest neutral Higgs boson of 1.5 GeV for \(\tan \beta = 5\). This phase dependence is gradually washed out for increasing \(\tan \beta \).

In Fig. 10, the CP-odd percentage is shown with \(P_{23}^2\) and \(P_{33}^2\) defined in Eq. (70) and evaluated at tree level with the 2HDM parameters at the scale \(M_{H^+}\). Figure 10 demonstrates that the mixing of the two heavy neutral Higgs bosons is very nearly independent of \(M_{H^+}\). Only for \(M_{H^+}= 200\) GeV can slight deviations from the other results be observed. For \(M_{H^+}= 200\) GeV, the CP-oddness is also shared with the lightest neutral Higgs boson, while for larger \(M_{H^+}\) the CP-odd component is mostly shared between the second lightest and heavy Higgs boson. One can also read off of Fig. 10 that for real values of |A| i.e. for \(\varphi _A = 0,180, 360^\circ \), the heavy Higgs boson \(H_3\) is CP-odd and \(H_2\) is CP-even. However, for \(\varphi _A \approx 90,\, 270^\circ \), the next-to lightest neutral Higgs boson is mostly CP-odd. That means both Higgs bosons oscillate between \(\sim 100\%\) CP-odd and CP-even with the phase \(\varphi _A\).

In addition, it is shown in Fig. 10 that the loop effects are negligible for the mixing of the heavy Higgs bosons, since the dashed lines, which include loop effects, mostly overlap with the solid lines, which depict the tree-level result. Not shown are differences between the different methods of evaluating the CP-mixing, i.e. whether the parameters are evaluated at \(M_{H^+}\) or at \(m_t\), or whether \({\mathcal {M}}^2\) or \((\mathcal M^2)^{\text {resum}}\) is used. The differences are on the same order as the difference shown between the tree-level and the one-loop results, in other words negligible.

7 Conclusion

In this paper, we have explored the Higgs sector of the CP-violating MSSM in a mass scenario with heavy SUSY particles and light Higgs bosons using effective field theory techniques. We matched the complex MSSM to a type-III 2HDM (where both Higgs doublets couple to the right-handed top and the right-handed bottom quarks), and calculated the complex threshold corrections to the 2HDM at one-loop level including contributions proportional to the top Yukawa, the bottom Yukawa and the strong coupling for a common soft-SUSY breaking mass scale for the superpartner particles and the corresponding RGEs for the 2HDM with complex parameters at two-loop level. Using these matching conditions and evolving the parameters down to a low scale with these RGEs, we resum the respective contributions at NLL order. We explored the effect of including the complex phases of the quartic couplings in the RGEs and found that in particular the absolute values of \(\lambda _5\) to \(\lambda _7\) change substantially compared to a scenario where only the absolute values (and signs) but not the phases are included in the RGEs.

In order to calculate the pole masses of the Higgs bosons, we exploited different methods:

  1. (a)

    Calculation of the pole mass with the 2HDM with parameters at the scale of the mass of the charged Higgs boson \(M_{H^+}\).

  2. (b)

    Calculation of the pole mass with the 2HDM with parameters at the scale of the running top quark mass.

  3. (c)

    Matching to the SM at the scale \(M_{H^+}\) and calculation of the pole mass of the Higgs boson at the scale of the running top quark mass within the SM.

  4. (d)

    Exploiting an approximation which resums the most important logarithms of the form \(\ln (\frac{M_{H^+}}{m_t})\) but still uses the full 2HDM for the calculation of the pole masses.

The approximation (d) agrees well with the pure 2HDM result for small values of \(M_{H^+}\) and approaches the result where the SM is used as the low-scale effective theory for larger \(M_{H^+}\). Therefore, it can be used as an interpolation between the two regimes.

We investigated several different scenarios to discuss the phase dependence of both the masses of the Higgs bosons and the CP-violating components of the Higgs bosons. All the masses of the Higgs bosons show a sizeable dependence on the common phase \(\varphi _A\) of \(A_t\) and \(A_b\), which can be of the order of several GeV for low \(\tan \beta \), in particular for the mass of the lightest Higgs boson. The heaviest Higgs boson shows the least sensitivity to \(\varphi _A\). The dependence of the Higgs masses on the gluino phase is much weaker and of the order of one GeV.

Additionally, we found that the size of the CP-odd component of the two heavy Higgs bosons shows a negligible dependence on the mass of the charged Higgs bosons, and that the heavy Higgs bosons interchange their CP-oddness when varying \(\varphi _A\). For the light Higgs boson, the size of the CP-odd component decreases quickly with larger values of \(M_{H^+}\) as has been discussed before in e.g. Ref. [98]. Even though we find that the CP-odd admixture in the considered scenario for \(M_{H^+} = 500\) GeV is just below the expected experimental reach according to the discussion in Ref. [131], it is likely that this particular scenario is excluded by the experimental results for the Higgs signal rates as well as of the measurement of electric dipole moments. In order to find out whether there is a viable scenario with heavy SUSY partners that leads to an observable size of a CP-odd component, further work is needed.