Phenomenology of on-shell Higgs production in the MSSM with complex parameters

Stefan Liebler¹,
Shruti Patel¹ &
Georg Weiglein¹

1196 Accesses
16 Citations
1 Altmetric
Explore all metrics

Abstract

A computation of inclusive cross sections for neutral Higgs boson production through gluon fusion and bottom-quark annihilation is presented in the MSSM with complex parameters. The predictions for the gluon-fusion process are based on an explicit calculation of the leading-order cross section for arbitrary complex parameters which is supplemented by higher-order corrections: massive top- and bottom-quark contributions at NLO QCD, in the heavy top-quark effective theory the top-quark contribution up to N$^{3}$LO QCD including a soft expansion for the $\mathcal {CP}$-even component of the light Higgs boson. For its $\mathcal {CP}$-odd component and the heavy Higgs bosons the contributions are incorporated up to NNLO QCD. Two-loop electroweak effects are also incorporated, and SUSY QCD corrections at NLO are interpolated from the MSSM with real parameters. Finite wave function normalisation factors ensuring correct on-shell properties of the external Higgs bosons are incorporated from the code FeynHiggs. For the typical case of a strong admixture of the two heavy Higgs bosons it is demonstrated that squark effects are strongly dependent on the phases of the complex parameters. The remaining theoretical uncertainties for cross sections are discussed. The results have been implemented into an extension of the numerical code SusHi called SusHiMi.

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

Article Open access 14 July 2018

Heavy Higgs decays into sfermions in the complex MSSM: a full one-loop analysis

Article Open access 07 May 2015

Neutral Higgs production at proton colliders in the CP-conserving NMSSM

Article Open access 14 May 2015

Use our pre-submission checklist

Avoid common mistakes on your manuscript.

1 Introduction

In 2012 the experimental collaborations ATLAS and CMS announced the discovery of a Higgs-like boson [1, 2] produced in collisions of protons at the Large Hadron Collider (LHC). Apart from the precise measurement of its production and decay properties in order to test whether there are deviations from the expectations for a Standard Model (SM) Higgs boson, an essential part of the programme of the LHC experiments in the upcoming years will be the search for additional Higgs bosons. The observed state can be easily accommodated in extended Higgs sectors like a Two-Higgs-Doublet Model (2HDM) or supersymmetric extensions, e.g. the Minimal Supersymmetric Standard Model (MSSM). For the search for additional Higgs bosons and the test of deviations from the SM expectations for the SM-like Higgs boson, the precise knowledge of production cross sections through gluon fusion and bottom-quark annihilation for these Higgs bosons is a key ingredient. Current efforts in this direction are summarised in the reports of the LHC Higgs Cross Section Working Group, see Refs. [3,4,5,6]. So far, the searches for additional Higgs bosons have been interpreted in various scenarios beyond the Standard Model, including several supersymmetric ones. However, those analyses do not yet cover the most general case where $\mathcal {CP}$ is violated and leads to mixing between $\mathcal {CP}$-even and -odd eigenstates. The reason that an analysis for the general case taking account the possibility of $\mathcal {CP}$-violation has not been possible so far has mainly been the lack of appropriate theoretical predictions for the Higgs production rates at the LHC for complex parameters in the MSSM and of a practical prescription for taking into account relevant interference effects in Higgs production and decay. A discussion of the latter has recently been given in Ref. [7]. Therefore it is our goal in the present paper to provide state-of-the-art cross-section predictions in the MSSM, taking into account $\mathcal {CP}$-violating effects, for the two main Higgs production channels at the LHC, which can be used as input for future experimental analyses in $\mathcal {CP}$-violating Higgs scenarios. We present in this paper precise predictions for neutral Higgs boson production through gluon fusion and bottom-quark annihilation in the MSSM with complex parameters, in which $\mathcal {CP}$-even and $\mathcal {CP}$-odd Higgs states form three admixed Higgs mass eigenstates ${h_a},a\in \lbrace 1,2,3\rbrace $. Complex parameters in the MSSM give rise to additional sources of $\mathcal {CP}$ violation beyond the one induced by the mixing of the quarks of the SM, described by the Cabibbo–Kobayashi–Maskawa (CKM) matrix [8, 9]. In order to explain the baryon asymmetry of the universe, such additional sources of $\mathcal {CP}$ violation beyond the CKM phase are actually needed, see e.g. Refs. [10,11,12] for reviews. It is thus of interest to investigate the MSSM with complex parameters. Its Higgs sector is influenced by the additional phases only beyond tree level. Still, these phases are of relevance in the Higgs boson collider phenomenology as they can induce a large mixing among the heavy Higgs bosons, and squark and gluino loop contributions also directly affect Higgs boson production and decay.

For a brief summary of higher-order corrections to the most important production processes – gluon fusion and bottom-quark annihilation – in the SM and the MSSM with real parameters we refer to Sect. 3 and focus here on studies performed for the Higgs sector of the MSSM with complex parameters. Early investigations of Higgs production through gluon fusion at hadron colliders in the MSSM with complex parameters were carried out in Refs. [13,14,15]. A thorough analysis taking different production channels into account was presented in Ref. [16], and results for Higgsstrahlung can be found in Ref. [17]. Large effects of stops on the cross section for a $\mathcal {CP}$-odd Higgs boson neglecting $\mathcal {CP}$-even and -odd Higgs mixing were discussed in Ref. [18]. References [19, 20] discuss the production of a light Higgs through gluon fusion including its decay into two photons in the MSSM with complex parameters. It should be noted that the mentioned references were published before the Higgs discovery in 2012 and mostly employ only the lowest order in perturbation theory for the production processes. It is therefore timely to improve these predictions by including up-to-date higher-order corrections and to investigate the compatibility with the experimental results obtained for the observed signal at 125 GeV. For this purpose we incorporate the prediction within the MSSM with complex parameters into the numerical code SusHi [21, 22], which calculates Higgs production through gluon fusion and heavy-quark annihilation [23] in the SM, the MSSM, the Two-Higgs-Doublet-Model (2HDM) and the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [24]. However until now, SusHi did not support complex parameters in the MSSM and thus did not provide predictions for $\mathcal {CP}$-admixed Higgs bosons.

For the calculation of the masses and the wave function normalisation factors ensuring the correct on-shell properties of external Higgs bosons, which involves the evaluation of Higgs boson self-energies and their renormalisation, we use the code FeynHiggs [25,26,27,28,29]. It employs a Feynman-diagrammatic approach and includes the full one-loop [28] and the dominant two-loop corrections of $\mathcal {O}(\alpha _t\alpha _s)$ [30] and $\mathcal {O}(\alpha _t^2)$ [31, 32] in the MSSM with complex parameters.^{Footnote 1} A detailed description of the prediction for the Higgs boson masses and the wave function normalisation factors as implemented in FeynHiggs can be found in Refs. [7, 28, 35,36,37,38,39]. Whereas the Higgs sector at tree level remains $\mathcal {CP}$-conserving, at higher orders an admixture of all three neutral Higgs bosons, i.e. the two $\mathcal {CP}$-even Higgs bosons h, H and the $\mathcal {CP}$-odd Higgs boson A, is induced. The case where the light Higgs boson describes the SM-like Higgs at ${\sim }125$ GeV is typically accompanied with a strong admixture of the two heavy Higgs bosons. For a proper prediction in such a case interference effects need to be taken into account in the full process involving production and decay of the Higgs bosons, which requires going beyond the usual narrow-width approximation (see also Refs. [40,41,42,43,44]). A convenient way to incorporate interference effects is a generalised narrow-width approximation for the production and decay of on-shell particles as described in Refs. [7, 45, 46], where in Ref. [46] only lowest-order contributions have been considered, while in Refs. [7, 45] also the inclusion of higher-order corrections has been addressed. The results for the cross sections for on-shell Higgs boson production obtained in the present paper are suitable for direct incorporation into the framework of a generalised narrow-width approximation.

Our paper is organised as follows: We start by outlining the relevant quantities in the Higgs, the gluino and the squark sector of the MSSM with complex parameters in Sect. 2. We move to the description of the gluon-fusion cross section in Sect. 3, where we discuss the calculation of the cross section at LO and the applicability of higher-order corrections. Next we introduce in Sect. 4 the code SusHi and its extension SusHiMi, which we use for our phenomenological studies carried out in Sect. 5. We discuss the remaining theoretical uncertainties in Sect. 6. Lastly, we conclude in Sect. 7 and list Higgs-(s)quark couplings in Appendix A.

2 The MSSM with complex parameters

In this section we discuss the relevant sectors of the MSSM with complex parameters, namely the gluino, the squark as well as the Higgs sector. While the discussion of the gluino and the squark sector at tree level is sufficient for our purposes, we will briefly describe the inclusion of higher-order corrections in the Higgs sector. The MSSM with complex parameters allows for 12 physical, independent phases of the complex parameters, once the phases of the wino soft-breaking parameter $M_2$ and the soft-breaking parameter $m_{12}^2$ are rotated away. Those independent phases are the ones of the soft-breaking gaugino masses $M_1$ and $M_3$, the Higgsino mass parameter $\mu $ and trilinear soft-breaking couplings $A_f, f\in \lbrace e,\mu ,\tau ,u,d,c,s,t,b\rbrace $. In the subsequent discussion we focus on these phases and their effect on the gluino, the squark and the Higgs sector as well as the Higgs boson cross sections.

2.1 Gluino and squark sector

The gluino $\tilde{g}$ does not mix with other fields and enters the Lagrangian in the form

$$\begin{aligned} {\mathcal {L}}\supset -\frac{1}{2}\overline{\tilde{g}}m_{\tilde{g}}\tilde{g}, \end{aligned}$$

(1)

where $m_{\tilde{g}}$ is the absolute value of the complex soft-breaking parameter $M_3=m_{\tilde{g}}e^{i\phi _{M_3}}$.^{Footnote 2} In the Feynman diagrams for the Higgs boson self-energies and the Higgs boson production via gluon fusion, the gluino only contributes beyond the one-loop level. However it affects the bottom-quark Yukawa coupling already at the one-loop level, where it enters the leading corrections to the relation between the bottom-quark mass and the bottom-quark Yukawa coupling which can be resummed to all orders, see below.

In the MSSM without flavour mixing in the squark sector, squarks $\tilde{q}_{L,R}$ of one generation mix into mass eigenstates $\tilde{q}_{1,2}$. The term of the Lagrangian containing the squark mass matrix of one generation is given by [37]

$$\begin{aligned} \mathcal {L}&\supset -(\tilde{q}_L^\dagger , \tilde{q}_R^\dagger )M_{\tilde{q}}^2\begin{pmatrix}\tilde{q}_L\\ \tilde{q}_R\end{pmatrix}\qquad \text {with}\nonumber \\&M_{\tilde{q},11}^2 = M^2_{\tilde{q}_L} + m_q^2 + M^2_Z \cos 2 \beta (I^3_q - Q_q s_W^2) \nonumber \\&M_{\tilde{q},12}^2 = m_q X^*_q \nonumber \\&M_{\tilde{q},21}^2 = m_q X_q \nonumber \\&M_{\tilde{q},22}^2 = M^2_{\tilde{q}_R} + m^2_q + M_Z^2 \cos 2 \beta Q_q s^2_W. \end{aligned}$$

(2)

Here $X_q := A_q - \mu ^* \cdot \lbrace { \cot \beta , \tan \beta \rbrace }$, where $\cot \beta $ and $\tan \beta $ apply to up- and down-type quarks, respectively. The soft-breaking masses $M^2_{\tilde{q}_L}$ and $M^2_{\tilde{q}_R}$, the third component of the weak isospin $I^3_q$, the electric charge $Q_q$ and the mass of the quark $m_q$ are real parameters. This also applies to the Z-boson mass $M_Z$ and the sine of the weak mixing angle $s_W\equiv \sin \theta _W$. Contrarily, in the $\mathcal {CP}$-violating MSSM the parameters $A_q = |A_q| e^{i \phi _{A_q}}$ and $\mu = |\mu | e^{i \phi _{\mu }}$, and hence $X_q$, can be complex. These complex parameters enter the Higgs sector via the Higgs-sfermion couplings, see Appendix A, which are thus also of direct relevance for Higgs boson production.

The mass matrix is diagonalised through the unitary matrix $U_{\tilde{q}}$ having real diagonal elements and complex off-diagonal elements

$$\begin{aligned} \begin{pmatrix} \tilde{q}_1 \\ \tilde{q}_2 \end{pmatrix} = U_{\tilde{q}} \begin{pmatrix} \tilde{q}_L \\ \tilde{q}_R \end{pmatrix}. \end{aligned}$$

(3)

The squark masses (using the convention $m_{\tilde{q}_1} \le m_{\tilde{q}_2}$) are calculated as the eigenvalues of Eq. (2). The fact that the left-handed soft-breaking parameter $M^2_{\tilde{q}_L}$ is the same for the fields in an SU(2) doublet gives rise to a tree-level relation between the stop and the sbottom masses. At the loop level, the corresponding relation between the physical squark masses receives a finite shift, see Ref. [47], which we have incorporated as a shift in the left-handed soft-breaking parameter $M^2_{\tilde{q}_L}$ in the sbottom sector, as obtained from FeynHiggs.

In the $b / \tilde{b}$ sector we take into account higher-order corrections to the relation between the bottom-quark mass and the bottom-Yukawa coupling [48,49,50,51,52,53]. The couplings of the lowest-order mass eigenstates $\phi $, where $\phi \in \lbrace h,H,A\rbrace $, see Sect. 2.2 below, of the Higgs bosons to bottom quarks are given by the effective Lagrangian

$$\begin{aligned} \mathcal {L}_{\text {eff}}= & {} \frac{m_be}{v_d \sqrt{2}} \sum _{\phi ^e \in \lbrace h,H \rbrace } \bar{b} \left[ g_{b_L}^{\phi ^e} P_L + (g_{b_L}^{\phi ^e})^* P_R \right] b\phi ^e \nonumber \\&+\,\frac{im_be}{v_d \sqrt{2}} \bar{b} \left[ g_{b_L}^{A} P_L - (g_{b_L}^{A})^* P_R \right] bA \end{aligned}$$

(4)

in terms of the left-handed and right-handed couplings $g_{b_L}^{\phi }$ and $g_{b_R}^{\phi }=(g_{b_L}^{\phi })^*$, where $P_{L/R}=\frac{1}{2}(1\mp \gamma _5)$ are the left- and right-handed projection operators, respectively. The explicit form of the couplings is given by

$$\begin{aligned} g_{b_L}^{h}= & {} \frac{f_{\alpha \beta }^h}{1+ \Delta _b} \left[ 1 - \frac{\cot \alpha }{\tan \beta } \Delta _b \right] ,\nonumber \\ g_{b_L}^{H}= & {} \frac{f_{\alpha \beta }^H}{1+ \Delta _b} \left[ 1 + \frac{\tan \alpha }{\tan \beta } \Delta _b \right] ,\nonumber \\ g_{b_L}^{A}= & {} \frac{f_{\alpha \beta }^A}{1+ \Delta _b} \left[ 1- \frac{\Delta _b}{\tan ^2 \beta } \right] , \end{aligned}$$

(5)

with $f_{\alpha \beta }^h=\sin \alpha /\cos \beta $, $f_{\alpha \beta }^H=\cos \alpha /\cos \beta $ and $f_{\alpha \beta }^A=\tan \beta $ (see also Refs. [37, 54]). The effective Lagrangian provides a resummation of leading $\tan \beta $-enhanced contributions entering via the quantity $\Delta _b$. The leading QCD contribution to $\Delta _b$ has the form

$$\begin{aligned} \Delta _b = \frac{2}{3}\frac{\alpha _s(\mu _d)}{\pi }M_3^*\mu ^*\tan \beta \, I\left( m_{\tilde{b}_1}^2,m_{\tilde{b}_2}^2,m_{\tilde{g}}^2\right) , \end{aligned}$$

(6)

where $\alpha _s$ is typically evaluated at an averaged SUSY scale $\mu _d=(m_{\tilde{b}_1}+m_{\tilde{b}_2}+m_{\tilde{g}})/3$, and the function I(a, b, c) is given by $I(a,b,c)=(ab\log (\tfrac{a}{b})+bc\log (\tfrac{b}{c})+ca\log (\tfrac{c}{a} ))/((a-b)(b-c)(a-c)))$. As one can see from Eq. (6), the leading contribution to $\Delta _b$ has an explicit dependence on the complex parameters $M_3$ and $\mu $. In our numerical analysis below we use the value for $\Delta _b$ as obtained from FeynHiggs (see Ref. [55]), which includes additional QCD and electroweak contributions [56,57,58,59]. In our implementation in the program SusHiMi, see Sect. 4 below, both the $\Delta _b$ value from FeynHiggs and the leading contribution from Eq. (6) can be selected.

We will use the expression for the bottom-quark Yukawa coupling according to the effective Lagrangian of Eqs. (4) and (5) in our leading-order expressions for the (loop-induced) gluon-fusion process. For bottom-quark annihilation and the implementation of higher-order corrections to the gluon-fusion process, see Sect. 3, we use as a simplified version [37]

$$\begin{aligned} g_{b}^\phi \equiv g_{b_L}^\phi =g_{b_R}^\phi =\frac{1}{|1+\Delta _b|}f^\phi _{\alpha \beta }, \end{aligned}$$

(7)

in which the left- and right-handed couplings to bottom quarks are identical to each other. We will compare the numerical impact of the two implementations at LO in Sect. 5. The effective Yukawa coupling in Eq. (4) is complex. The phase of this coupling could be rotated away by an appropriate redefinition of the (s)quark fields, see e.g. Ref. [56]. We prefer to use the general expression for a complex Yukawa coupling. In our phenomenological discussion in Sect. 5 below we will compare the effect of the complex Yukawa coupling of Eq. (4) with the simplified real coupling of Eq. (7) (which are not equivalent to each other) and we will show that the numerical differences are small.

2.2 Higgs sector

The MSSM contains two Higgs doublets with opposite hypercharges $Y_{{\mathcal {H}}_{1,2}} = \pm 1$ in order to introduce masses for both the up- and down-type fermions. The neutral fields of the two Higgs doublets can be decomposed in $\mathcal {CP}$-even ($\phi _1^0, \phi _2^0$) and $\mathcal {CP}$-odd ($\chi _1^0, \chi _2^0$) components as follows^{Footnote 3}

$$\begin{aligned} {\mathcal {H}}_1= & {} \begin{pmatrix} h_d^0 \\ h_d^- \end{pmatrix} = \begin{pmatrix} v_d + \frac{1}{\sqrt{2}}(\phi _1^0 + i \chi _1^0)\\ \phi ^-_1 \end{pmatrix} \end{aligned}$$

(8)

$$\begin{aligned} {\mathcal {H}}_2= & {} \begin{pmatrix} h_u^+ \\ h_u^0 \end{pmatrix} = e^{i \xi } \begin{pmatrix} \phi ^+_2\\ v_u + \frac{1}{\sqrt{2}}(\phi _2^0 + i \chi _2^0) \end{pmatrix} , \end{aligned}$$

(9)

such that the Higgs potential $V_H$ in terms of the neutral Higgs states is given by

$$\begin{aligned} V_H^0= & {} \left( |\mu |^2 + m^2_{{\mathcal {H}}_2}\right) |h^0_u|^2 + \left( |\mu |^2 + m^2_{{\mathcal {H}}_1}\right) |h^0_d|^2 \nonumber \\&-\,[m_{12}^2 h_u^0h_d^0 +\hbox {h.c.}] + \frac{g_1^2+ g_2^2}{8}[|h_u^0|^2-|h_d^0|^2]^2.\nonumber \\ \end{aligned}$$

(10)

The quadratic terms of $V_H$ contain the SUSY parameter $|\mu |^2$ and the soft terms $m_{{\mathcal {H}}_1}$, $m_{{\mathcal {H}}_2}$. The bilinear terms have the soft coefficient $m^2_{12}$, which is a complex parameter in general but whose phase can be absorbed through a Peccei–Quinn transformation [60, 61]. The relative phase $\xi $ between the Higgs doublets vanishes when the Higgs potential is minimised, making the Higgs sector of the MSSM $\mathcal {CP}$-invariant at lowest order.

The tree-level neutral mass eigenstates $\lbrace h, H, A, G\rbrace $ are related to the tree-level neutral fields $\lbrace \phi _1^0, \phi _2^0, \chi _1^0, \chi _2^0\rbrace $ through a unitary matrix as follows

$$\begin{aligned} \left( \begin{array}{c} h\\ H \\ A \\ GV \end{array}\right) = \left( \begin{array}{c@{\quad }c@{\quad }c@{\quad }c} -s_{\alpha } &{} c_{\alpha } &{} 0 &{} 0 \\ c_{\alpha } &{} s_{\alpha } &{} 0 &{} 0 \\ 0 &{} 0 &{} s_{\beta _n} &{} c_{\beta _n}\\ 0 &{} 0 &{} -c_{\beta _n} &{} s_{\beta _n} \end{array}\right) \left( \begin{array}{c} \phi _1^0 \\ \phi _2^0 \\ \chi _1^0 \\ \chi _2^0 \end{array}\right) . \end{aligned}$$

(11)

Similarly, for the charged Higgs states one obtains

$$\begin{aligned} \left( \begin{array}{c} H^{\pm } \\ G^{\pm } \end{array}\right) = \left( \begin{array}{cc} s_{\beta _c} &{} c_{\beta _c} \\ -c_{\beta _c} &{} s_{\beta _c} \end{array}\right) \left( \begin{array}{c} \phi _1^{\pm } \\ \phi _2^{\pm } \end{array}\right) , \end{aligned}$$

(12)

where $s_x \equiv \sin x, c_x \equiv \cos x$. $\alpha $, $\beta _n$ and $\beta _c$ are the mixing angles for the $\mathcal {CP}$-even Higgs bosons (h, H), the neutral $\mathcal {CP}$-odd states (A, G), and the charged states ($H^{\pm }$, $G^{\pm }$), respectively. Minimising the Higgs potential leads to $\beta :=\beta _n = \beta _c$ at tree level. The masses of the charged Higgs bosons and the neutral $\mathcal {CP}$-odd Higgs boson at tree level are given by

$$\begin{aligned} m^2_{H^{\pm }} = m_A^2 + M_W^2,\quad m_A^2= \frac{2 m^2_{12}}{\sin (2 \beta )}. \end{aligned}$$

(13)

The Higgs sector of the MSSM at lowest order is fully determined (besides the gauge couplings) by two parameters, which are usually chosen as $m_{H^{\pm }}$ ($m_A$) and $\tan \beta := \frac{v_u}{v_d}$ for the case of the MSSM with complex (real) parameters.

2.3 Higgs mixing at higher orders

$\mathcal {CP}$-violating mixing between the neutral Higgs bosons $\lbrace h,H,A\rbrace $ arises as a consequence of radiative corrections and results in the neutral mass eigenstates $\lbrace h_1, h_2, h_3\rbrace $, where by convention $m_{h_1} \le m_{h_2} \le m_{h_3}$. The full mixing in higher orders takes place not just between $\lbrace h, H, A\rbrace $, but also with the Goldstone boson and the electroweak gauge bosons. In general, $(6 \times 6)$-mixing contributions involving the fields $\left\{ h,H,A,G,Z,\gamma \right\} $ need to be taken into account. For the calculation of the Higgs boson masses and wave function normalisation factors at the considered order it is sufficient to restrict to a $(3 \times 3)$-mixing matrix among $\left\{ h,H,A \right\} $, since mixing effects with $\lbrace G,Z,\gamma \rbrace $ only appear at the sub-leading two-loop level and beyond. In processes with external Higgs bosons, on the other hand, mixing contributions with G and Z already enter at the one-loop level, but the numerical effect of these contributions has been found to be very small, see e.g. Refs. [35,36,37, 62]. In our numerical analysis of the Higgs production through gluon fusion and bottom-quark annihilation below we will neglect these kinds of (electroweak) mixing contributions of the external Higgs bosons with Goldstone and gauge bosons. Concerning electroweak corrections, we only incorporate the potentially numerically large contributions to the Higgs boson masses and wave function normalisation factors as well as the electroweak contribution to the correction affecting the relation between the bottom-Yukawa coupling and the bottom quark mass (see above), while all other contributions considered here like e.g. electroweak corrections to gluon fusion involve at least one power of the strong coupling. For the contribution of the Z boson and the Goldstone boson to the gluon-fusion process via $gg\rightarrow \lbrace Z^*,G^*\rbrace \rightarrow h_i$ (the photon only enters at higher orders) it should be noted that contributions from mass-degenerate quark weak-isodoublets vanish and only top- and bottom-quark contributions proportional to their masses are of relevance, see the discussion of the Higgsstrahlung process in Refs. [63, 64]. This is a consequence of the fact that only the axial component of the quark-quark-Z boson coupling contributes to the loop-induced coupling of the Z boson to two gluons. Similarly, squark contributions in $gg\rightarrow \lbrace Z^*, G^*\rbrace $ are completely absent at the one-loop level, even in case of $\mathcal {CP}$ violation in the squark sector. The one-loop contributions to $gg\rightarrow \lbrace Z^*, G^*\rbrace $ therefore have no dependence on the phases of complex parameters.

Thus, we focus our discussion on the contributions to the $(3\times 3)$-mass matrix ${\mathbf {M}}$, which contains the tree level masses $m_i^2$ on the diagonal and has non-zero (off-)diagonal self-energies involving the Higgs states. It enters the Lagrangian, with $\Phi =\left( h,H,A\right) $, as follows

$$\begin{aligned}&{\mathcal {L}}\supset -\frac{1}{2}\Phi {\mathbf {M}}\Phi ^T \quad \text {with} \nonumber \\&\quad {\mathbf {M}}=\left( \begin{array}{c@{\quad }c@{\quad }c} m_h^2 - {\hat{\Sigma }}_{hh}(p^2) &{} -{\hat{\Sigma }}_{hH}(p^2) &{} - {\hat{\Sigma }}_{hA}(p^2) \\ - {\hat{\Sigma }}_{Hh}(p^2) &{} m_H^2 - {\hat{\Sigma }}_{HH}(p^2) &{} - {\hat{\Sigma }}_{HA}(p^2) \\ - {\hat{\Sigma }}_{Ah}(p^2) &{} - {\hat{\Sigma }}_{AH}(p^2) &{} m_A^2 - {\hat{\Sigma }}_{AA}(p^2) \end{array}\right) .\nonumber \\ \end{aligned}$$

(14)

The propagator matrix is then given by

$$\begin{aligned}&[-\varvec{\Delta }_{hHA} (p^2)]^{-1} ={\hat{\varvec{\Gamma }}}_{hHA}(p^2) \quad \text {with}\nonumber \\&\quad [{\hat{\varvec{\Gamma }}}_{hHA}(p^2)]_{ij}=\hat{\Gamma }_{ij} = i[(p^2 m_i^2) \delta _{ij} + {\hat{\Sigma }}_{ij}(p^2)], \end{aligned}$$

(15)

and the roots of the determinant of this matrix yield the loop-corrected Higgs boson masses. The non-diagonal and diagonal propagators can be written as follows

$$\begin{aligned}&\Delta _{ij} = \frac{\hat{\Gamma }_{ij}\hat{\Gamma }_{kk}-\hat{\Gamma }_{jk}\hat{\Gamma }_{ki}}{\hat{\Gamma }_{ii}\hat{\Gamma }_{jj}\hat{\Gamma }_{kk} +2\hat{\Gamma }_{ij}\hat{\Gamma }_{jk}\hat{\Gamma }_{ki}-\hat{\Gamma }_{ii}\hat{\Gamma }_{jk}^2-\hat{\Gamma }_{jj}\hat{\Gamma }_{ki}^2-\hat{\Gamma }_{kk}\hat{\Gamma }_{ij}^2}, \nonumber \\ \end{aligned}$$

(16)

$$\begin{aligned}&\Delta _{ii} = \frac{i}{p^2 - m_i^2 + {\hat{\Sigma }}_{ii}^{\text {eff}}}, \end{aligned}$$

(17)

with effective self-energies that contain mixed terms

$$\begin{aligned} {\hat{\Sigma }}_{ii}^{\text {eff}}= & {} {\hat{\Sigma }}_{ii} - i \frac{2\hat{\Gamma }_{ij}\hat{\Gamma }_{jk} \hat{\Gamma }_{ki}-\hat{\Gamma }_{jj}\hat{\Gamma }_{ki}^2-\hat{\Gamma }_{kk}\hat{\Gamma }_{ij}^2}{\hat{\Gamma }_{jj}\hat{\Gamma }_{kk}-\hat{\Gamma }_{jk}^2}\nonumber \\= & {} {\hat{\Sigma }}_{ii} + \frac{\Delta _{ij}}{\Delta _{ii}} {\hat{\Sigma }}_{ij} + \frac{\Delta _{ik}}{\Delta _{ii}} {\hat{\Sigma }}_{ik}, \end{aligned}$$

(18)

where we suppressed the $p^2$ arguments of all terms and $i\ne j\ne k$.

2.4 Wave function normalisation factors for external Higgs bosons

For Higgs bosons that appear as external particles in a process appropriate on-shell properties are required for a correct normalisation of the S-matrix. Unless the field renormalisation constants have been chosen such that all mixing contributions between the mass eigenstates $\lbrace h_1, h_2, h_3\rbrace $ vanish on-shell and the propagators of the external particles have unit residue, the correct on-shell properties need to be ensured via the introduction of finite wave function normalisation factors, see e.g. Refs. [65,66,67], as a consequence of the LSZ formalism [68]. The matrix of those so-called ${\hat{\mathbf Z}}$ factors contains the correction factors for the external Higgs bosons $\lbrace h_1, h_2, h_3\rbrace $ relative to the lowest-order mass eigenstates $\lbrace h, H, A\rbrace $. The matrix elements $\hat{{\mathbf {Z}}}_{aj}$ [39] (see also Refs. [35, 37]) are composed of the root of the external wave function normalisation factor

$$\begin{aligned} \hat{Z}_i^a := \text {Res}_{\mathcal {M}_a^2} \lbrace \Delta _{ii} (p^2) \rbrace \end{aligned}$$

(19)

and the on-shell transition ratio

$$\begin{aligned} \hat{Z}^a_{ij} = \frac{\Delta _{ij}(p^2)}{\Delta _{jj}(p^2)} \Bigg |_{p^2 = \mathcal {M}_a^2}, \end{aligned}$$

(20)

which are evaluated at the complex pole $\mathcal {M}_a^2$. Here the indices $\lbrace a, b, c \rbrace $ refer to the loop-corrected mass eigenstates, while $\lbrace i, j, k \rbrace $ label the lowest-order mass eigenstates. With an appropriate assignment of the indices of the two types of states (see Ref. [39]) the matrix elements can be written as

$$\begin{aligned} \hat{{\mathbf {Z}}}_{aj} = \sqrt{\hat{Z}}_a \hat{Z}_{aj}, \end{aligned}$$

(21)

corresponding to the (non-unitary) matrix

$$\begin{aligned} {\hat{\mathbf Z}}= \left( \begin{array}{c@{\quad }c@{\quad }c} \sqrt{\hat{Z}}_1 \hat{Z}_{1h} &{} \sqrt{\hat{Z}}_1 \hat{Z}_{1H} &{} \sqrt{\hat{Z}}_1 \hat{Z}_{1A} \\ \sqrt{\hat{Z}}_2 \hat{Z}_{2h} &{} \sqrt{\hat{Z}}_2 \hat{Z}_{2H} &{} \sqrt{\hat{Z}}_2 \hat{Z}_{2A} \\ \sqrt{\hat{Z}}_3 \hat{Z}_{3h} &{} \sqrt{\hat{Z}}_3 \hat{Z}_{3H} &{} \sqrt{\hat{Z}}_3 \hat{Z}_{3A} \end{array}\right) . \end{aligned}$$

(22)

As explained above, these ${\hat{\mathbf Z}}$ factors provide the correct normalisation of a matrix element with an external on-shell Higgs boson $h_a$, $a\in \lbrace 1,2,3\rbrace $, at $p^2 = \mathcal {M}_a^2$. The application of the ${\hat{\mathbf Z}}$ factors yields an expression of the amplitude ${\mathcal {A}}_{h_a}$ for an external on-shell Higgs boson $h_a$ in terms of a linear combination of the amplitudes resulting from the one-particle irreducible diagrams for each of the lowest-order mass eigenstates $\lbrace h, H, A\rbrace $ according to

$$\begin{aligned} {\mathcal {A}}_{h_a}= & {} {\hat{\mathbf Z}}_{ah}{\mathcal {A}}_{h} + {\hat{\mathbf Z}}_{aH} {\mathcal {A}}_{H} + {\hat{\mathbf Z}}_{aA} {\mathcal {A}}_{A}+\cdots \nonumber \\= & {} \sqrt{\hat{Z}_a} ( \hat{Z}_{ah} {\mathcal {A}}_h + \hat{Z}_{aH}{\mathcal {A}}_H + \hat{Z}_{aA}{\mathcal {A}}_A)+\cdots . \end{aligned}$$

(23)

The ellipsis indicate additional mixing effects with Goldstone bosons and gauge bosons, which we neglect in our numerical analysis, see the discussion in Sect. 2.3.

3 The gluon-fusion cross section

In this section we discuss the calculation of the gluon-fusion cross section with particular emphasis on the effects of complex parameters. We first focus on individual ingredients and then combine them in Sect. 4. For this purpose our notation closely follows Ref. [21]. At leading order (LO) the gluon-fusion cross section is known since a long time [69]. In addition to the quark-induced contributions, the squark-induced contributions to the gluon-fusion process are also of relevance in supersymmetric extensions of the SM, even though they are suppressed by inverse powers of the supersymmetric particle masses if those masses are heavy. Subsequently we present our calculation of the LO cross section for the case of the MSSM with complex parameters for the three physical Higgs bosons $h_a, a\in \lbrace 1,2,3\rbrace $. Differences with respect to the calculation in the MSSM with real parameters are induced through^{Footnote 4}

${\hat{\mathbf Z}}$ factors, which relate the amplitude for an external on-shell Higgs $h_a$ (in the mass eigenstate basis) to the amplitudes of both the $\mathcal {CP}$-even lowest-order states h and H and the $\mathcal {CP}$-odd state A, see Sect. 2.4.
Non-vanishing couplings of squarks $g^A_{\tilde{f}ii}$ to the pseudoscalar component A.
Different left- and right-handed quark couplings $g^\phi _{q_L}$ and $g^\phi _{q_R}$ with $\phi \in \lbrace h,H,A\rbrace $, see Sect. 2.1.

3.1 Lowest-order cross section

The LO production cross section of the mass eigenstates $h_a$ can be written as follows

$$\begin{aligned}&\sigma _{{\mathrm{LO}{}}}(pp\rightarrow h_a)=\sigma _0^{h_a}\tau _{h_a}{\mathcal {L}}^{gg}(\tau _{h_a})\quad \text {with}\nonumber \\&\quad {\mathcal {L}}^{gg}(\tau )=\int _{\tau }^1\frac{\mathrm{d}x}{x}g(x)g(\tau /x), \end{aligned}$$

(24)

where $\tau _{h_a}=m_{h_a}^2/s$. The hadronic squared centre-of-mass energy is denoted by s, and the gluon-gluon luminosity by ${\mathcal {L}}^{gg}$. Therein, the partonic LO cross section for $gg\rightarrow h_a$ is given by

$$\begin{aligned} \sigma _0^{h_a}= & {} \frac{G_F \alpha _s^2 (\mu _R)}{288 \sqrt{\pi }}\left[ \left| \mathcal {A}^{h_a,{\mathrm {e}}} \right| ^2 + \left| \mathcal {A}^{h_a,{\mathrm {o}}} \right| ^2\right] \nonumber \\&\quad \text {with }\mathcal {A}^{h_a,{\mathrm {e}}}={\hat{\mathbf Z}}_{ah}\mathcal {A}_+^{h} + {\hat{\mathbf Z}}_{aH} \mathcal {A}_+^{H} + {\hat{\mathbf Z}}_{aA} \mathcal {A}_-^{A}\nonumber \\&\quad \text {and }\mathcal {A}^{h_a,{\mathrm {o}}}={\hat{\mathbf Z}}_{ah}\mathcal {A}_-^{h} + {\hat{\mathbf Z}}_{aH} \mathcal {A}_-^{H} + {\hat{\mathbf Z}}_{aA} \mathcal {A}_+^{A},\nonumber \\ \end{aligned}$$

(25)

where $G_F$ denotes Fermi’s constant, and ${\hat{\mathbf Z}}_{a\phi }$ are the elements of the ${\hat{\mathbf Z}}$ factor matrix. $\mu _\text {R}$ is the renormalisation scale, which at LO only enters through the scale dependence of the strong coupling constant $\alpha _s$. We denote the cross section $pp\rightarrow h_a$, which involves one-loop diagrams in the production process $pp\rightarrow \phi $, as “LO cross section” despite the fact that it contains higher-order effects through the application of the ${\hat{\mathbf Z}}$ factors (Fig. 1). We note that in the effective field theory approach of heavy quark and SUSY masses, where the gluon-gluon-Higgs interaction is condensed into a single vertex, the amplitudes of the first term in Eq. (25) can be identified with a contribution that stems from ${\mathcal {L}}\supset G^{\mu \nu }G_{\mu \nu }\phi $ with the gluon field strength $G^{\mu \nu }$. The amplitudes of the second term stem from ${\mathcal {L}}\supset \tilde{G}^{\mu \nu }G_{\mu \nu }\phi $, which involves the dual of the gluon field strength tensor $\tilde{G}^{\mu \nu }$, resulting in the cross section being expressible as the sum of two non-interfering squared amplitudes. This explains the naming of the first and the second term with $\mathcal {A}^{h_a,{\mathrm {e}}}$ and $\mathcal {A}^{h_a,{\mathrm {o}}}$, respectively. Similarly, we can split $\sigma _{{\mathrm{LO}{}}}$ into $\sigma _{{\mathrm{LO}{}}}^{{\mathrm {e}}}$ and $\sigma _{{\mathrm{LO}{}}}^{{\mathrm {o}}}$.

For the two $\mathcal {CP}$-even lowest-order mass eigenstates $\phi ^e\in \lbrace h,H\rbrace $ we obtain the amplitudes

$$\begin{aligned} \mathcal {A}_+^{\phi ^e} = \sum _{q \in \lbrace t,b \rbrace } \left( a^{\phi ^e}_{q,+} + \tilde{a}_{q}^{\phi ^e}\right) ,\quad \mathcal {A}_-^{\phi ^e} = \sum _{q \in \lbrace t,b \rbrace } a^{\phi ^e}_{q,-} \end{aligned}$$

(26)

with

$$\begin{aligned} a^{\phi ^e}_{q,+}= & {} \frac{1}{2} \left( g_{q_L}^{\phi ^e} + g_{q_R}^{\phi ^e} \right) \frac{3}{2} \tau _q^{h_a} \left[ 1 + \left( 1- \tau _q^{h_a}\right) f(\tau _q^{h_a}) \right] ,\nonumber \\ a^{\phi ^e}_{q,-}= & {} \frac{i}{2}\left( g_{q_R}^{\phi ^e} - g_{q_L}^{\phi ^e} \right) \frac{3}{2} \tau _q^{h_a}f\left( \tau _q^{h_a}\right) \nonumber \\ \tilde{a}_{q}^{\phi ^e}= & {} -\frac{3}{8} \tau _q^{h_a} \sum _{i=1}^2 g^{\phi ^e}_{\tilde{q}ii} \left[ 1- \tau _{\tilde{q}i}^{h_a} f\left( \tau _{\tilde{q}i}^{h_a}\right) \right] , \end{aligned}$$

(27)

where $g_{q_L}^{\phi ^e}$ and $g_{q_R}^{\phi ^e}$ are the couplings of the Higgs boson $\phi ^e$ to the left- and right-handed quarks, respectively. They are normalised to the SM Higgs-quark couplings. $g^{\phi ^e}_{\tilde{q}ij}$ are the couplings of the Higgs boson $\phi ^e$ to squarks $\tilde{q}_i$ and $\tilde{q}_j$. The explicit expressions for the Higgs-squark and relative Higgs-quark couplings are listed in Appendix A. Similarly, for the $\mathcal {CP}$-odd Higgs boson A we have

$$\begin{aligned} \mathcal {A}_-^{A} = \sum _{q \in \lbrace t,b \rbrace } \left( a^{A}_{q,-} + \tilde{a}_{q}^{A}\right) ,\quad \mathcal {A}_+^{A} = \sum _{q \in \lbrace t,b \rbrace } a^{A}_{q,+} \end{aligned}$$

(28)

with

$$\begin{aligned} a^{A}_{q,+}= & {} \frac{1}{2} \left( g_{q_{L}}^A + g_{q_R}^A \right) \frac{3}{2} \tau _q^{h_a} f(\tau _q^{h_a}),\nonumber \\ a^{A}_{q,-}= & {} \frac{i}{2} \left( g_{q_L}^A - g_{q_R}^A \right) \frac{3}{2} \tau _q^{h_a} \left[ 1 + \left( 1- \tau _q^{h_a}\right) f(\tau _q^{h_a}) \right] \nonumber \\ \tilde{a}_{q}^{A}= & {} -\frac{3}{8}\tau _{q}^A \sum _{i=1}^2 g_{\tilde{q}ii}^A \left[ 1- \tau _{\tilde{q}_i}^{h_a} f\left( \tau _{\tilde{q}_i}^{h_a}\right) \right] . \end{aligned}$$

(29)

Within the previous formulas we use the notation

$$\begin{aligned} \tau _q^{h_a} = \frac{4m^2_q}{m_{h_a}^2}, \tau _{\tilde{q}_i}^{h_a} = \frac{4m_{\tilde{q}_i}^2}{m_{h_a}^2}, \end{aligned}$$

(30)

and $f(\tau )$ is given by

$$\begin{aligned} f(\tau ) = \left\{ \begin{array}{ll} \text {arcsin}^2 \frac{1}{\sqrt{\tau }} &{}\quad {\mathrm {for}}\ \tau \ge 1\\ -\frac{1}{4} \left( \log \frac{1+ \sqrt{1-\tau }}{1-\sqrt{1-\tau }} - i\pi \right) ^2 &{} \quad {\mathrm {for}}\ \tau < 1 \end{array}. \right. \end{aligned}$$

(31)

Our result is consistent with Ref. [13], which however assumes $g_{q_L}^\phi =g_{q_R}^\phi $ (see our discussion of this issue in Sect. 2.1) and does not take into account the mixing among the tree-level mass eigenstates $\phi \in \lbrace h,H,A\rbrace $. All squark contributions, i.e. $\tilde{a}_q^{\phi ^e}$ and $\tilde{a}_q^{A}$, enter the first term, $\mathcal {A}^{h_a,{\mathrm {e}}}$, in Eq. (25). Quark contributions to $\mathcal {A}^{h_a,{\mathrm {e}}}$ which couple to the $\mathcal {CP}$-odd lowest-order mass eigenstate A are proportional to the difference between the left- and right-handed quark Yukawa couplings. The same holds for the contributions to the second term $\mathcal {A}^{h_a,{\mathrm {o}}}$ in Eq. (25) which couple to the $\mathcal {CP}$-even lowest-order mass eigenstates $\phi ^e$. All these terms are therefore denoted with the subscript ${\mathcal {A}}_-$. It should be noted that the amplitudes ${\mathcal {A}}_-^{\phi _e, A}$ only arise due to the complex nature of the Yukawa couplings, which is a consequence of the incorporation of higher-order contributions entering via $\Delta _b$, see Eq. (5), and our choice of working with a complex Yukawa coupling. Accordingly, the amplitudes ${\mathcal {A}}_-^{\phi _e, A}$ are zero in the case of the MSSM with real parameters.

3.2 Higher-order contributions

Gluon fusion receives sizeable corrections at higher orders in QCD. The NLO corrections for the SM quark contributions are known for arbitrary quark masses [70,71,72,73,74,75]. NNLO (SM-) QCD contributions were calculated in the limit of a heavy top-quark mass [76,77,78], similar to the recently published N$^{3}$LO contributions for a $\mathcal {CP}$-even Higgs boson in an expansion around the threshold of Higgs production [79,80,81,82,83].^{Footnote 5} Finite top-quark mass effects at NNLO are known in an expansion of inverse powers of the top-quark mass [86,87,88,89,90,91,92,93]. All of the previously mentioned corrections are implemented in SusHi [21, 22] and can be added in all supported models. We will later discuss in more detail for which Higgs mass ranges these corrections are applicable, which also explains why the above mentioned N$^{3}$LO contributions are only employed for the $\mathcal {CP}$-even component of the light Higgs boson.

As explained above a complex Yukawa coupling is only induced for the bottom quark through the incorporation of $\Delta _b$ contributions. According to this approach, for the top-quark Yukawa coupling $g_t^\phi $ left- and right handed components are identical also in the MSSM with complex parameters. Therefore we can directly adapt the known higher-order QCD corrections to the top-quark loop contribution for the MSSM with complex parameters. They are incorporated in the extension SusHiMi, see Sect. 4. For the incorporation of the bottom-quark contribution at NLO (SM-) QCD, on the other hand, we have to rely on the simplified version of the $\Delta _b$ corrections to the bottom-Yukawa coupling as specified in Eq. (7). Electroweak two-loop corrections as discussed in Refs. [94,95,96] can be added as well. We take into account the contributions mediated by light quarks, which can be reweighted to the MSSM with complex parameters. We follow Ref. [97] and define the correction factor

$$\begin{aligned} \delta _{\text {EW}}^\text {lf}=\frac{\alpha _{\text {EM}}}{\pi } \frac{2\text {Re}\left( {\mathcal {A}}^{h_a,{\mathrm {e}}} {\mathcal {A}}^{h_a,\text {EW}*}\right) }{|{\mathcal {A}}^{{h_a},{\mathrm {e}}}|^2}, \end{aligned}$$

(32)

where ${\mathcal {A}}^{h_a,{\mathrm {e}}}$, which has been given in Eq. (25), denotes the $\mathcal {CP}$-even part of the LO amplitude including quark and squark contributions. Accordingly, this electroweak correction factor is only applied to the $\mathcal {CP}$-even component of the LO and NLO cross section, see Sect. 4. The electroweak amplitude is given by [96]

(33)

with the abbreviation

$$\begin{aligned} x_V = \frac{1}{m_{h_a}^2}\left( M_V-i\frac{\Gamma _V}{2}\right) ^2, \quad V\in \{W,Z\}. \end{aligned}$$

(34)

In Eq. (32) $\alpha _{\text {EM}}$ denotes the electro-magnetic coupling, and $s_W\equiv \sin \theta _W=(1-c_W^2)^{1/2}=(1-M_W^2/M_Z^2)^{1/2}$ is the sine of the weak mixing angle. $M_V$ and $\Gamma _V$ are the mass and the width of the heavy gauge bosons $V\in \{W,Z\}$, and the function $A_1[x]$ can be found in Ref. [96].

In the MSSM with real parameters analytical NLO virtual contributions involving squarks, quarks and gluinos are either known in the limit of a vanishing Higgs mass [98,99,100,101] or in an expansion of heavy SUSY masses [102,103,104].^{Footnote 6} Even NNLO corrections of stop-induced contributions to gluon fusion are known [108, 109]; SusHi can approximate these NNLO stop effects [110] in the $\mathcal {CP}$-conserving MSSM. We neglect those contributions in our analysis for the MSSM with complex parameters.

At NLO in the MSSM with complex parameters, supersymmetric contributions are present both in virtual and real corrections. The real corrections show a similar behaviour as observed for the LO cross section, i.e. the squark induced contributions of $\mathcal {CP}$-odd components proportional to $g_{\tilde{q}ii}^A$ are added as a complex component to the $\mathcal {CP}$-even couplings. Since beyond LO we employ the simplified version of the $\Delta _b$ resummation according to Eq. (7), the higher-order quark contributions, both real and virtual, are of the same structure as in the $\mathcal {CP}$-conserving MSSM. The NLO virtual contributions as described in the previous paragraph are however not easily adjustable to the MSSM with complex parameters. We therefore interpolate the NLO virtual contributions between phases 0 and $\pi $ of the various MSSM parameters using a cosine interpolation, see Refs. [111, 112]. This interpolation makes use of on-shell stop- and sbottom-quark masses defined at phases 0 and $\pi $. Thus, within the interpolated result we have to ensure the correct subtraction of the NLO contributions that have already been taken into account through $\Delta _b$ effects in the bottom-quark Yukawa coupling. This is done by expanding the $\Delta _b$ correction to next-to-leading order in the subtraction term. For a certain value of the phase $\phi _z$ of a complex parameter z, the virtual NLO amplitude $\mathcal {A}_{\mathrm{NLO}{}}^{\phi } (\phi _z)$ can be approximated using

$$\begin{aligned} \mathcal {A}_{\mathrm{NLO}{}}^{\phi } (\phi _z) = \frac{1 + \cos \phi _z}{2}\mathcal {A}_{\mathrm{NLO}{}}^{\phi } (0) + \frac{1 - \cos \phi _z}{2}\mathcal {A}_{\mathrm{NLO}{}}^{\phi } ( \pi )\nonumber \\ \end{aligned}$$

(35)

for each of the lowest-order mass eigenstates $\phi \in \lbrace h,H,A\rbrace $. Here $\mathcal {A}_{\mathrm{NLO}{}}^{\phi }(0)$ is the analytical result for the MSSM with real parameters, and $\mathcal {A}_{\mathrm{NLO}{}}^{\phi }(\pi )$ is the analytical result with $z \rightarrow -z$. Using the factors $\cos \phi _z$ ensures a smooth interpolation such that the known results for a vanishing phase are recovered. Whereas a dependence on the phases of $A_q$ and $\mu $ is already apparent in the lowest-order diagrams of $gg\rightarrow \phi $, the phase of $M_3$ only enters through the NLO virtual corrections. Besides the $\Delta _b$ contributions, where the full phase dependence is incorporated, the treatment of the phase of $M_3$ therefore relies on the performed interpolation. While the implemented routines for the MSSM with real parameters are expressed in terms of the gluino mass, they can also be used for a negative soft-breaking parameter $M_3$, such that we can obtain interpolated results for a complex-valued parameter $M_3$. We note that the NLO virtual amplitudes with a negative $M_3$ are identical to the virtual amplitudes for positive $M_3$ with opposite signs of the parameters $A_t, A_b$ and $\mu $. This can be understood from the structure of the NLO diagrams involving the squark–quark–gluino couplings. It should however be noted in this context that due to the generation of Higgs-squark couplings $g_{\tilde{q}ii}^A$ for non-vanishing phases a new class of NLO virtual diagrams arises which is not present in the MSSM with real parameters. Since the interpolation is based on the result for the MSSM with real parameters as input for the predictions at the phases 0 and $\pi $, the additional set of diagrams may not be adequately approximated in this way.

Despite this fact, we expect that the interpolation of the virtual two-loop contributions involving squarks and gluinos to the gluon-fusion amplitude provides a reasonable approximation, for the following reasons (we discuss the theoretical uncertainty associated with the interpolation in Sect. 6 and assign a conservative estimate of the uncertainty in our numerical analysis). We focus here on the gluon fusion amplitude without ${\hat{\mathbf Z}}$ factors, since in the ${\hat{\mathbf Z}}$ factors the full phase dependence is incorporated without approximations. Gluino contributions are generally suppressed for gluino masses that are sufficiently heavy to be in accordance with the present bounds from LHC searches, while gluon-exchange contributions do not add an additional phase dependence compared to the dependence on the phases of $A_q$ and $\mu $ in the LO cross section, which is fully taken into account. The dependence of the NLO amplitude on the phases of $A_q$ and $\mu $ is therefore expected to follow a similar pattern as the LO amplitude, which is also what we find in the application of the interpolation method.

One can also compare the higher-order corrections to the gluon-fusion process with the ones to the Higgs boson masses and ${\hat{\mathbf Z}}$ factors. In fact, a similar interpolation was probed in the prediction for Higgs boson masses in the MSSM with complex parameters, see e.g. Refs. [28, 30, 112, 113], where the phase dependence of sub-leading two-loop contributions beyond $\mathcal{O}(\alpha _t \alpha _s)$ were approximated with an interpolation before the full phase dependence of the corresponding two-loop corrections at $\mathcal{O}(\alpha _t^2)$ was calculated [31, 32]. Generally good agreement was found between the full result and the approximation [31, 32]. In order to investigate the interpolation of the phase of $M_3$ that is associated with the gluino we performed a similar check concerning the phase dependence of two-loop squark and gluino loop contributions. We numerically compared the full result for the Higgs mass prediction at this order from FeynHiggs with an approximation where the phases at the two-loop level are interpolated. Despite the fact that also for the Higgs mass calculation new diagrams proportional to $g_{\tilde{q}ii}^A$ arise away from phases 0 and $\pi $, the phase dependence of the interpolated results generically follows the behaviour of the full results very well.

Based on the NLO amplitude that has been obtained as described above, we can construct the NLO cross sections $\sigma _{\mathrm{NLO}{}}^{{\mathrm {e}}}$ and $\sigma _{\mathrm{NLO}{}}^{{\mathrm {o}}}$ individually, following Ref. [21], by defining the NLO correction factors $C^{{\mathrm {e}}}$ and $C^{{\mathrm {o}}}$:

$$\begin{aligned} C^{{\mathrm {e/o}}}=2\text {Re}\left[ \frac{\mathcal {A}_{{\mathrm{NLO}}}^{h_a,{\mathrm {e/o}}}}{\mathcal {A}^{h_a,{\mathrm {e/o}}}}\right] +\pi ^2+\beta _0\log \left( \frac{\mu _\text {R}^2}{\mu _\text {F}^2}\right) . \end{aligned}$$

(36)

The amplitudes are given by $\mathcal {A}_{{\mathrm{NLO}}}^{h_a,{\mathrm {e}}}={\hat{\mathbf Z}}_{ah}\mathcal {A}_{{\mathrm{NLO}}}^h+{\hat{\mathbf Z}}_{aH}\mathcal {A}_{{\mathrm{NLO}}}^H$ and $\mathcal {A}_{{\mathrm{NLO}}}^{h_a,{\mathrm {o}}}={\hat{\mathbf Z}}_{aA}\mathcal {A}_{{\mathrm{NLO}}}^A$, $\mu _\text {F}$ denotes the factorisation scale, and $\beta _0=11/2-n_f/3$ with $n_f=5$. Note that the LO amplitudes $\mathcal {A}^{h_a,{\mathrm {e/o}}}$ entering Eq. (36) are taken in the limit of large stop and sbottom masses, see Ref. [21]. The correction factors enter the NLO cross section as follows

$$\begin{aligned}&\sigma _{{\mathrm{NLO}{}}}^{{\mathrm {e/o}}}(pp\rightarrow h_a+X)=\sigma _0^{h_a,{\mathrm {e/o}}}\tau _{h_a}{\mathcal {L}}^{gg}(\tau _{h_a})\left[ 1+C^{{\mathrm {e/o}}}\frac{\alpha _s}{\pi }\right] \nonumber \\&\quad +\,\Delta \sigma _{gg}^{{\mathrm {e/o}}}+\Delta \sigma _{gq}^{{\mathrm {e/o}}}+\Delta \sigma _{q\bar{q}}^{{\mathrm {e/o}}}. \end{aligned}$$

(37)

The terms $\Delta \sigma $ denote the real corrections, which are not fully displayed here. We emphasise again that at NLO we work with Eq. (7), such that in the real corrections the only new ingredients are Higgs-squark couplings $g_{\tilde{q}ii}^A$, which are added to the $\mathcal {CP}$-even components $\Delta \sigma ^{{\mathrm {e}}}$. The real corrections can be split in $\Delta \sigma ^{{\mathrm {e}}}$ and $\Delta \sigma ^{{\mathrm {o}}}$ since no interference terms arise.

4 The program SusHi and the extension SusHiMi

SusHi is a numerical FORTRAN code [21, 22] which combines analytical results for the calculation of Higgs boson cross sections through gluon fusion and heavy-quark annihilation in models beyond the Standard Model up to the highest known orders in perturbation theory. However, the current release does not allow for $\mathcal {CP}$ violation in the Higgs sector. Following our discussion in Sect. 3 we present the calculation of Higgs boson production in the context of the MSSM with complex parameters, which we included in an extension of SusHi named SusHiMi.^{Footnote 7} For this purpose we proceed along the lines of Fig. 2 and calculate the Higgs boson production cross section through gluon fusion as follows: SusHiMi calls SusHi twice and in these two calls performs a “$\mathcal {CP}$-even” calculation for $\sigma _{{\mathrm{NLO}{}}}^{{\mathrm {e}}}$ and a “$\mathcal {CP}$-odd” calculation for $\sigma _{{\mathrm{NLO}{}}}^{{\mathrm {o}}}$ according to Eq. (37). Thus, the total gluon-fusion cross section is the sum of the two parts

$$\begin{aligned}&\sigma _{\mathrm{N}^{k}\mathrm{LO}}(pp\rightarrow h_a+X)=\sigma _{\mathrm{N}^{k}\mathrm{LO}}^{{\mathrm {e}}}(pp\rightarrow h_a+X)\nonumber \\&\quad +\,\sigma _{\mathrm{N}^{k}\mathrm{LO}}^{{\mathrm {o}}}(pp\rightarrow h_a+X). \end{aligned}$$

(38)

We obtain the result beyond LO QCD through^{Footnote 8}

$$\begin{aligned}&\sigma ^{{\mathrm {e}}}_{{\mathrm{N}^{k}\mathrm{LO}}} = \sigma _{{\mathrm{NLO}{}}}^{{\mathrm {e}}}(1+\delta _{\mathrm {EW}}^{\mathrm {lf}}) +\left( \sigma _{{\mathrm{N}^{k}\mathrm{LO}}, {\mathrm{EFT}{}}}^{t,{\mathrm {e}}} - \sigma _{{\mathrm{NLO}{}, {\mathrm{EFT}{}}}}^{t,{\mathrm {e}}}\right) \end{aligned}$$

(39)

$$\begin{aligned}&\sigma ^{{\mathrm {o}}}_{{\mathrm{N}^{k}\mathrm{LO}}} = \sigma _{{\mathrm{NLO}{}}}^{{\mathrm {o}}} +\left( \sigma _{{{\mathrm{N}^{k}\mathrm{LO}},{\mathrm{EFT}{}}}}^{t,{\mathrm {o}}} - \sigma _{{\mathrm{NLO}{}, {\mathrm{EFT}{}}}}^{t,{\mathrm {o}}}\right) , \end{aligned}$$

(40)

whereas $\sigma ^{{\mathrm {e/o}}}_{{\mathrm{LO}{}}}$ was specified in Sect. 3, and $k\in \lbrace 1,2,3\rbrace $. The $\mathcal {CP}$-odd component $\sigma _{{{\mathrm{N}^{k}\mathrm{LO}}, {\mathrm{EFT}{}}}}^{t,{\mathrm {o}}}$ is only implemented up to $k=2$ (see below). In the previous formulas $\sigma _{{\mathrm{NLO}{}}}^{{\mathrm {e/o}}}$ are the NLO cross sections including real contributions and the interpolated NLO virtual corrections as discussed in Sect. 3. They employ the simplified $\Delta _b$ resummation according to Eq. (7), named $\Delta _{b1}$. $\sigma _{{{\mathrm{N}^{k}\mathrm{LO}}, {\mathrm{EFT}{}}}}^{t,{\mathrm {e}}}$ and $\sigma _{{{\mathrm{N}^{k}\mathrm{LO}}, {\mathrm{EFT}{}}}}^{t,{\mathrm {o}}}$ are cross sections including the top-quark contribution only. They are based on a K-factor calculated in the EFT approach of an infinitely heavy top-quark obtained for a SM Higgs boson H and a pseudoscalar A (in a 2HDM with $\tan \beta =1$) with mass $m_{h_a}$, respectively. This K-factor is subsequently reweighted with the exact LO cross section. For this purpose the employed LO cross sections $\sigma _{{\mathrm{LO}}}^{t,{\mathrm {e}}}$ and $\sigma _{{\mathrm{LO}}}^{t,{\mathrm {o}}}$ are again evaluated as discussed in Sect. 3 with full ${\hat{\mathbf Z}}$ factors, but include only the top-quark contribution. They are multiplied with the K-factors in $\sigma _{{{\mathrm{N}^{k}\mathrm{LO}}, {\mathrm{EFT}{}}}}^{t,{\mathrm {e}}}$ and $\sigma _{{{\mathrm{N}^{k}\mathrm{LO}}, {\mathrm{EFT}{}}}}^{t,{\mathrm {o}}}$, respectively. Due to their small numerical impact in $\sigma _{{{\mathrm{N}^{k}\mathrm{LO}}, {\mathrm{EFT}{}}}}^{t,{\mathrm {e/o}}}$ we do not take into account top-quark mass effects beyond NLO even though they are implemented in SusHi. An alternative approach, which is not discussed in this paper but can be implemented in SusHiMi, is to include the relative couplings $g_t^\phi $ and the ${\hat{\mathbf Z}}$ factors into the complex-valued Wilson coefficients of the EFT directly.

As already mentioned N$^{3}$LO QCD corrections are only taken into account for the $\mathcal {CP}$-even component of the light Higgs boson, which allows us to match the precision of the light Higgs boson cross section in the SM employed in up-to-date predictions. This is motivated by the fact that the light Higgs boson that is identified with the observed signal at $125~\text {GeV}$ is usually assumed to have a dominant $\mathcal {CP}$-even component, which is also the case in the scenarios which are considered in our numerical discussion. For the $\mathcal {CP}$-odd component of the light Higgs and the heavy Higgs bosons we employ the NNLO corrections for the top-quark induced contributions to gluon fusion in the effective theory of a heavy top-quark, i.e. we do not take into account top-quark mass effects beyond NLO, but only factor out the LO QCD cross sections $\sigma _{{\mathrm{LO}}}^{t,{\mathrm {e}}}$ and $\sigma _{{\mathrm{LO}}}^{t,{\mathrm {o}}}$. The strategy to employ the EFT result at NNLO beyond the top-quark mass threshold can be justified from the comparison of NLO corrections, which are known in the EFT approach and exactly with full quark-mass dependence and agree also beyond the top-quark mass threshold. On the other hand, the N$^{3}$LO QCD corrections that were obtained for the top-quark contribution are only known in the EFT approach and for an expansion around the threshold of Higgs production at $x=m_{h_a}^2/s\rightarrow 1$, which we can take into account up to $\mathcal {O}(1-x)^{16}$. Since the combination of the EFT approach and the threshold expansion becomes questionable above the top-quark mass threshold, we apply N$^{3}$LO QCD corrections only for the $\mathcal {CP}$-even component of the light Higgs boson and thus match the precision of the SM prediction. The electroweak correction factor $\delta _{\mathrm {EW}}^{\mathrm {lf}}$ multiplied in the “$\mathcal {CP}$-even” run is obtained from Eq. (32).

As shown in Fig. 2 we call SusHiMi three times in order to take into account the different possibilities of the resummation of $\tan \beta $ enhanced sbottom effects in the LO QCD contributions. We add the results as follows

$$\begin{aligned} \sigma (pp\rightarrow h_a+X)=\sigma _{{{\mathrm{N}^{k}\mathrm{LO}}}}^{\Delta _{b1}} +\sigma _{{\mathrm{LO}{}}}^{\Delta _{b2}} -\sigma _{{\mathrm{LO}{}}}^{\Delta _{b1}}, \end{aligned}$$

(41)

where in the N$^{k}$LO QCD cross section following Eq. (38) the simplified resummation according to Eq. (7) is employed, indicated through the index $\Delta _{b1}$. We add and subtract the LO QCD cross section using the full resummation according to Eq. (5), named $\Delta _{b2}$, and the simplified resummation, respectively. As we will demonstrate the differences between the two versions of resummation are small, which can partially be understood from a possible rephasing of complex Yukawa couplings by a redefinition of all (s)quark fields (see the discussion in Sect. 2.1).

SusHi also allows one to obtain differential cross sections as a function of the transverse momentum or the (pseudo-)rapidity of the Higgs boson. These effects can be studied also in the MSSM with complex parameters. In the case of non-vanishing transverse momentum, which is only possible through additional radiation, i.e. real corrections, the precision for massive quark contributions in extended Higgs sectors is currently limited to the LO prediction [114, 115]. The predictions of the $p_T$ distributions in SusHiMi have been obtained from the LO contributions with arbitrary complex parameters, and in contrast to the total cross sections are therefore not affected by additional interpolation uncertainties from higher orders in comparison to the case of the MSSM with real parameters.

Higgs production through bottom-quark annihilation is calculated in SusHi for a SM Higgs boson at NNLO QCD. In the employed five-flavour scheme, where the bottom quarks are understood as partons, the result equals the cross section of a pseudoscalar A (in a 2HDM with $\tan \beta =1$). For the production of the Higgs boson $h_a$ in the MSSM with complex parameters, as implemented in SusHiMi, the results for the SM Higgs boson are reweighted to the MSSM with $|{\hat{\mathbf Z}}_{ah}g_b^h+{\hat{\mathbf Z}}_{aH}g_b^H|^2+|{\hat{\mathbf Z}}_{aA}g_b^A|^2$, which includes $\tan \beta $-enhanced squark effects through $\Delta _b$ according to Eq. (7). This procedure equals the application of a K-factor on the full LO cross section including ${\hat{\mathbf Z}}$ factors. In case of non-equal left- and right-handed couplings $g_{bL}$ and $g_{bR}$ due to the application of the full resummation in Eq. (5), the SM cross section has to be multiplied with

$$\begin{aligned}&|{\hat{\mathbf Z}}_{ah}(g_{bL}^h+g_{bR}^h)+{\hat{\mathbf Z}}_{aH}(g_{bL}^H+g_{bR}^H)+i{\hat{\mathbf Z}}_{aA}(g_{bL}^A-g_{bR}^A)|^2\nonumber \\&\quad +|i{\hat{\mathbf Z}}_{ah}(g_{bR}^h-g_{bL}^h)+i{\hat{\mathbf Z}}_{aH}(g_{bR}^H-g_{bL}^H)\nonumber \\&\quad +{\hat{\mathbf Z}}_{aA}(g_{bL}^A+g_{bR}^A)|^2. \end{aligned}$$

(42)

Though, due to the similarity of both approaches we only discuss Higgs production through bottom-quark annihilation with simplified $\Delta _b$ resummation.

5 Numerical results

For our numerical analysis we slightly modify two standard MSSM scenarios introduced in Ref. [116], namely the $m_h^{\mathrm {mod}+}$ and the light-stop scenario. The scenarios have been chosen for illustration, featuring relatively large squark and gluino contributions to the gluon fusion process. The corresponding effects will be relevant in our discussion of the associated theoretical uncertainties.

The light-stop inspired scenario that we use for our numerical analysis is defined as follows

$$\begin{aligned}&M_1=340~\text {GeV},\quad M_2= \mu =400~\text {GeV},\quad M_3=1.5~\text {TeV}\nonumber \\&X_t=X_b=X_{\tau }=1.0~\text {TeV},\quad A_q=A_l=0 \nonumber \\&{\tilde{m}}_{Q_2}={\tilde{m}}_L=1~\text {TeV}, \quad {\tilde{m}}_{Q_3}=0.5~\text {TeV}, \end{aligned}$$

(43)

where the modified values of $M_1$ and $M_2$ have been chosen to avoid direct bounds from stop searches obtained in LHC Run I (assuming R-parity conservation).^{Footnote 9} For the $m_h^{\mathrm {mod}+}$-inspired scenario we choose for vanishing phases of the complex parameters:

$$\begin{aligned}&M_1=250~\text {GeV},\quad M_2=500~\text {GeV},\quad M_3=1.5~\text {TeV}\nonumber \\&X_t=X_b=X_{\tau }=1.5~\text {TeV},\quad A_q=A_l=0 \nonumber \\&\mu ={\tilde{m}}_{Q}={\tilde{m}}_L=1~\text {TeV}. \end{aligned}$$

(44)

We use for the SM parameters the values $m_t^{\text {OS}} = 173.20$ GeV, $m_b^{\overline{\text {MS}}}(m_b)=4.16$ GeV, $m_b^{\text {OS}}(m_b)=4.75$ GeV and $\alpha _s (M_Z) = 0.119$. The depicted on-shell bottom-quark mass is used as internal mass for propagators and for the bottom-quark Yukawa coupling in the gluon-fusion process. The depicted value of $\alpha _s$ is only used for the evaluations of FeynHiggs, for the cross sections the value of $\alpha _s$ associated with the employed PDF set is taken. We employ the MMHT2014 PDF sets at LO, NLO and NNLO QCD [118]. The central choice for the renormalisation and factorisation scales $\mu _\text {R}^0$ and $\mu _\text {F}^0$, respectively, is $(\mu _\text {R}^0,\mu _\text {F}^0)=(m_{h_a}/2,m_{h_a}/2)$ for gluon fusion and $(\mu _\text {R}^0,\mu _\text {F}^0)=(m_{h_a},m_{h_a}/4)$ for bottom-quark annihilation. More details are described in Sect. 6.

Whereas for the $m_h^{\mathrm {mod}+}$-inspired scenario we pick heavy Higgs bosons through $m_{H^\pm }=900$ GeV with $\tan \beta =10$ and 40 for the study of $\Delta _b$ effects, we choose $m_{H^\pm }=500$ GeV with $\tan \beta =16$ for the light-stop inspired scenario. A detailed discussion of squark effects for the Higgs boson cross sections in the light-stop scenario can also be found in Ref. [119]. For the chosen parameter point the squark effects are sizeable, both for the light Higgs boson and in particular also for the heavy $\mathcal {CP}$-even Higgs boson, where they reduce the gluon-fusion cross section by about ${\sim }90$%. The Higgs boson masses and the ${\hat{\mathbf Z}}$ factors are obtained from FeynHiggs 2.11.2. The cross sections are evaluated with SusHiMi, which is based on the latest release of SusHi, version 1.6.1. We will mostly focus on the gluon-fusion cross section and present the bottom-quark annihilation cross section only for the scenario with $\tan \beta =40$.

For the parameter points associated with the mentioned scenarios in the MSSM with real parameters we vary the phases of $A_t=|A_t|e^{i\phi _{A_t}}$ and $M_3=m_{\tilde{g}}e^{i\phi _{M_3}}$ leaving the absolute values constant in order to address various aspects in the phenomenology of Higgs boson production. The phases of $A_b$ and $\mu $ do not introduce new phenomenological features, and we do not display results for the variation of those phases. A variation of the phase of $X_t$ leads to very similar cross sections for all Higgs bosons as observed for the variation of the phase of $A_t$. This can be understood from the fact that we choose not too large values of $\mu $ and $\tan \beta \ge 10$, and so $X_t \approx A_t$. Note that the stop masses are constant as a function of the phase of $X_t$, if the absolute value of $X_t$ is fixed. Before we proceed we want to briefly discuss experimental constraints on the phases: The most restrictive constraints on the phases arise from bounds on the electric dipole moments (EDMs) of the electron and the neutron, see Refs. [120,121,122] and references therein. EDMs from heavy quarks [123, 124] and the deuteron [125] also have an impact. MSSM contributions to these EDMs already contribute at the one-loop level and primarily involve the first two generations of sleptons and squarks. Thus, EDMs lead to severe constraints on the phases of $A_{q}$ for $q\in \{ u,d,s,c\}$ and $A_l$ for $l\in \{e,\mu \}$. Using the convention that the phase of the wino soft-breaking mass $M_2$ is rotated away, one finds tight constraints on the phase of $\mu $ [126]. On the other hand constraints on the phases of the third-generation trilinear couplings are significantly weaker. We refer the reader to Ref. [127] for a review. While recent constraints from EDMs [128] taking into account two-loop contributions [129] have the potential to rule out the largest values of the phase of $A_t$, there is still significant room for variation of the phases of $A_t$ and $M_3$. We therefore display the full range of the phases of $A_t$ and $M_3$ in our considered scenarios without explicitly imposing EDM constraints, following the common approach in benchmark scenarios for Higgs phenomenology (see e.g. Ref. [130] for a recent discussion). It should be noted in this context that in particular the variation of $A_t$ affects the value of the stop masses. Additionally, the Higgs boson masses are a function of the phases of the complex parameters. The impact is particularly pronounced for the mass of the light Higgs boson. In order to factor out the impact of phase space effects, we normalise the prediction for the cross section of the light Higgs boson in the MSSM to the cross section of a SM Higgs boson with identical mass as the light Higgs mass eigenstate $m_{h_1}$. In case of the heavy Higgs bosons for which the phase space effects are much less severe, we stick to the inclusive cross sections without such a normalisation. The predicted value for the Higgs boson mass $m_{h_1}$ deviates from 125 GeV by up to a few GeV in our illustrative studies. Deviations from the experimental value in this ballpark are still commensurate with the remaining theoretical uncertainties from unknown higher-order corrections of current state-of-the-art calculations of the light Higgs boson mass in the MSSM [6].

Subsequently we discuss three aspects: We start with a discussion of squark effects for the Higgs boson production cross sections. They are of relevance both for the heavy Higgs bosons and the light Higgs boson. Secondly, we focus on the admixture of the two heavy Higgs bosons (described through ${\hat{\mathbf Z}}$ factors) and its effect on production cross sections. Lastly we discuss $\Delta _b$ corrections in the context of the $m_h^{\mathrm {mod}+}$-inspired scenario with large $\tan \beta $, for which the bottom-quark annihilation process for the heavy Higgs bosons is relevant as well.

Note that given the large admixture of the two heavy Higgs bosons in the MSSM with complex parameters, interference effects in the full processes of production and decay can be large. However, we restrict our discussion in the present paper to Higgs boson production. The results for the cross sections obtained in our paper can be employed in a generalised narrow-width approximation as described in Ref. [7] in order to incorporate interference effects. We will address this issue elsewhere.

The prediction for Higgs boson cross sections is affected by various theoretical uncertainties, which we discuss in detail in Sect. 6. In order to demonstrate the improvement in precision through the inclusion of higher-order corrections, all subsequent figures which show the LO cross section and our best prediction cross section according to Eq. (38) include renormalisation and factorisation scale uncertainties. The procedure for obtaining these scale uncertainties is outlined in Sect. 6.

5.1 Squark contributions in the light-stop inspired scenario

We start with a discussion of squark effects to the Higgs boson cross section $\sigma (gg\rightarrow h_i)$ for all three Higgs bosons $h_i$ in the context of the light-stop inspired scenario with $m_{H^\pm }=500$ GeV and $\tan \beta =16$. The variation of the light Higgs boson mass $m_{h_1}$ as well as the stop masses is depicted in Fig. 3a as a function of $\phi _{A_t}$. The light Higgs mass in this scenario may appear to be too light to be compatible with the signal observed at the LHC, however we regard it as sufficiently close in view of the facts that our discussion should demonstrate phenomenological effects only and that on the other hand there are still sizeable theoretical uncertainties in the MSSM prediction for the light Higgs boson mass. The lightest stop mass has its minimum at around 293 GeV. For the heavy Higgs bosons with masses between 492 and 494 GeV, which are not shown in the figures, the NLO squark-gluino contributions [102,103,104] which assume heavy squarks and gluinos are thus well applicable. The variation of the heavy Higgs boson masses as a function of the phases of $\phi _{A_t}$ and $\phi _{M_3}$ turns out to be small, namely within 0.6 GeV. Due to the strong admixture of the left- and right-handed stops through a large value of $A_t$, also a phase dependence of the stop masses is observed. We checked that if we instead choose a phase for $X_t$ keeping $|X_t|$ constant, we obtain constant stop masses if the phase of $X_t$ is varied.

Figure 3b shows the production cross section through gluon fusion for the light Higgs boson $h_1$. The black, dot-dashed curve depicts the cross section with top-quark and bottom-quark contributions and electroweak corrections in the production amplitudes only, i.e. in the formulas of Sect. 4 we omit all squark contributions which enter either directly or through $\Delta _b$. Note however that squark contributions are always part of the ${\hat{\mathbf Z}}$ factors. Due to the decoupling with large values of $m_{H^\pm }$ in our scenarios the light Higgs $h_1$ has mostly SM-like couplings to quarks and gauge bosons. Thus, thanks to the inclusion of N$^{3}$LO QCD contributions for the top-quark induced contribution, our prediction of the gluon fusion cross section omitting the squark contributions (black, dot-dashed curve in Fig. 3b) is very close to the one for the SM Higgs boson with the same mass as provided by the LHC Higgs Cross Section Working Group [6, 131]. The inclusion of squark contributions explicitly and through $\Delta _b$ resummation lowers the gluon-fusion cross section by about 20%, as can be inferred from the blue, solid curve, which is shown together with its renormalisation and factorisation scale uncertainty, see Sect. 6. For completeness we also show the LO cross section calculated according to Eq. (24) including squark effects as the red curve. It is apparent that the scale uncertainties are significantly reduced from LO QCD to our best prediction cross section calculated according to Eq. (38). Figure 3b also includes the cross section for $h_1$ normalised to the cross section of a SM Higgs boson with the same mass. Here the ${\sim } 20$% reduction due to squark effects is apparent once again, whereas the quark-induced cross section shows the well-known decoupling behaviour. Not shown in the figures are the following effects, which we state here for completeness: The variation of $\phi _{M_3}$ leads to a very similar picture, even though the light Higgs mass variation is not as pronounced (2 GeV) and the stop masses are unaffected. Moreover, in the comparison of the simplified and the full resummation of $\Delta _b$ contributions in the LO gluon-fusion cross section of $h_1$ we observe a well-known behaviour, namely the simplified resummation of Eq. (7) does not yield a decoupled bottom-quark Yukawa coupling, whereas the full resummation of Eq. (5) does.

In Fig. 4a, b we show the gluon-fusion cross sections of the heavy Higgs bosons $h_2$ and $h_3$, respectively, as a function of $\phi _{A_t}$. The colour coding is identical to Fig. 3 except for the fact that we show the K-factor of our best prediction for the cross section with respect to the LO cross section, $\sigma /\sigma _{\mathrm{LO}{}}$, rather than a cross section normalised to the SM Higgs boson cross section. In fact, the heavy Higgs masses change only slightly as a function of the phase $\phi _{A_t}$, and therefore the associated phase space effect is small. For vanishing phase $\phi _{A_t}=0$ it is known that squark effects are huge and reduce the cross section by ${\sim } 89$% ($h_2$) and ${\sim } 22$% ($h_3$) [119]. These squark effects are strongly dependent on the phase $\phi _{A_t}$ and induce a large positive correction at phase $\phi _{A_t}=\pi $ in case of $h_2$. For $h_3$ the effects are not as pronounced, but still sizeable. The K-factor for both processes $gg\rightarrow h_2$ and $gg\rightarrow h_3$ remains within [1, 1.6], i.e. higher-order corrections mainly follow the phase dependence of the LO cross section. The dependence of the K-factor on $\phi _{A_t}$ follows the black, dot-dashed curve, which shows the cross section with quark contributions only. The significant dependence of the cross section where only quark contributions are included on the phase $\phi _{A_t}$ is induced by the admixture of the two Higgs bosons through ${\hat{\mathbf Z}}$ factors. We will discuss this feature in detail for the $m_h^{\mathrm {mod}+}$-inspired scenario in Sect. 5.2.

The phase dependence on $\phi _{M_3}$ is less pronounced. We show the corresponding cross sections for the two heavy Higgs bosons $h_2$ and $h_3$ in Fig. 5. As in previous figures we observe a significant reduction in the scale dependence from LO QCD to our best prediction for the cross section. The inclusion of squark and gluino contributions through the ${\hat{\mathbf Z}}$ factors and through $\Delta _b$ induces a dependence on the gluino phase already for the LO cross section. The almost flat black dot-dashed curves show the cross section with quark contributions only, and any variation with $\phi _{M_3}$ is an effect of the ${\hat{\mathbf Z}}$ factors, which in this case is negligible since $\phi _{M_3}$ only enters at the two-loop level. The K-factor, which takes into account our interpolated NLO virtual corrections, only shows a relatively mild dependence on the phase. We will discuss the interpolation uncertainty for this scenario in Sect. 6, since we obtain the largest relative interpolation uncertainty in the cross section variation with phases for the interpolation of the gluino phase $\phi _{M_3}$.

5.2 Admixture of Higgs bosons in the $m_h^{\mathrm {mod}+}$-inspired scenario

In this subsection we discuss the $m_h^{\mathrm {mod}+}$-inspired scenario with $\tan \beta =10$ and $m_{H^{\pm }}=900$ GeV. Since the squark masses are at the TeV level in this scenario, the numerical effect of the squark loops in the gluon fusion vertex contributions is rather small for the production cross section of the light Higgs boson $h_1$. We do not discuss the results for $h_1$ in this section. The results for the two heavy Higgs bosons are displayed in Fig. 6. The effects from squark loops are at the level of about ${\pm } 20$% in this case. The considered scenario is typical for the decoupling region of supersymmetric theories, where a light SM-like Higgs boson (that is interpreted as the signal observed at about 125 GeV) is accompanied by additional heavy Higgs bosons that are nearly mass-degenerate. In the general case where the possibility of $\mathcal {CP}$-violating interactions is taken into account, there can be a large mixing between the $\mathcal {CP}$-even and $\mathcal {CP}$-odd neutral Higgs states. This feature is clearly visible in Fig. 6. The dependence on the phase $\phi _{A_t}$ is seen to be closely correlated to the mixing character of the two neutral heavy Higgs bosons.

Figure 6a depicts the masses of the two heavy Higgs bosons $h_2$ and $h_3$ as a function of $\phi _{A_t}$ together with the $\mathcal {CP}$-odd character of $h_2$ and $h_3$, being defined as $|{\hat{\mathbf Z}}_{aA}|^2$. For illustration here and in the following we call the mass eigenstates $h_2$ and $h_3$ either $h_e$ or $h_o$, depending on their mixing character: if $|{\hat{\mathbf Z}}_{aA}|^2\gtrsim 1/2$ the mass eigenstate $h_a$ is called $h_o$, otherwise it is called $h_e$. It can be seen in Fig. 6b, c that the behaviour of the cross sections as a function of $\phi _{A_t}$ closely follows the variation in the $\mathcal {CP}$-even and $\mathcal {CP}$-odd character of the Higgs states. A similar effect was already apparent in the top- and bottom-quark induced cross sections depicted in the light-stop inspired scenario, see Fig. 4, however there the effects of squark contributions are dominant. Also in this case our best prediction for the cross section is significantly reduced in comparison with the prediction in LO QCD. The variation of the K-factors between about 1.2 and 1.5 with the phase $\phi _{A_t}$ also follows the modification of the mixing character of the two neutral heavy Higgs bosons.

Since the two heavy Higgs bosons are nearly mass degenerate, it may not be possible in such a case to experimentally resolve the two Higgs bosons as separate signals. Rather than the individual cross sections times their respective branching ratios, the experimentally measurable quantity then consists of the sum of the cross sections of the two Higgs states times their respective branching ratios together with the interference contribution involving the two Higgs states. The latter can be particularly important if the mass difference between the two Higgs states is smaller than the sum of their total widths [7]. While we defer the incorporation of such interference effects into the prediction for the production and decay process to a forthcoming publication, one can already infer from the plots of Fig. 6b, c that in the overall contribution there will be sizeable cancellations between the phase dependencies of the separate contributions.

5.3 $\Delta _b$ corrections in the $m_h^{\mathrm {mod}+}$-inspired scenario

We finally discuss the impact of $\Delta _b$ effects, which we investigate for the two heavy Higgs bosons in the $m_h^{\mathrm {mod}+}$-inspired scenario with $\tan \beta =40$. In this scenario the admixture between the two heavy Higgs bosons is again sizeable both as a function of $\phi _{A_t}$ and as a function of $\phi _{M_3}$. This even leads to mass crossings as seen in Fig. 7. It is therefore convenient to discuss the results in terms of the predominantly $\mathcal {CP}$-even mass eigenstate $h_e$ and the predominantly $\mathcal {CP}$-odd mass eigenstate $h_o$, as defined in Sect. 5.2, as for those states a smooth behaviour of the cross section as function of the phases is obtained. The masses of the two heavy Higgs bosons and their $\mathcal {CP}$-character (defining $h_o$ and $h_e$) are shown in Fig. 7 as a function of $\phi _{A_t}$ and $\phi _{M_3}$. One can see that the states $h_2$ and $h_3$ drastically change their $\mathcal {CP}$ character upon variation of the phases $\phi _{A_t}$ and $\phi _{M_3}$, while on the other hand the state $h_e$ is almost purely $\mathcal {CP}$-even and $h_o$ is almost purely $\mathcal {CP}$-odd for the whole range of phase values. It should be kept in mind in this context that $|{\hat{\mathbf Z}}_{aA}|^2$ arises from a non-unitary matrix and can therefore have values above 1. For vanishing phases the mass eigenstate $h_2$ corresponds to $h_e$ and $h_3$ to $h_o$.

In the following we show results for the predominantly $\mathcal {CP}$-even mass eigenstate $h_e$. The observations for $h_o$ are very similar and are not shown here, we will only add comments where appropriate. In Fig. 8 we show the gluon-fusion cross section as a function of the phases $\phi _{A_t}$ and $\phi _{M_3}$. In both cases the behaviour for the full prediction, including the squark contributions, is dominated by $\Delta _b$ corrections. For vanishing phases those corrections significantly reduce the cross sections compared to the case where only quark contributions are taken into account. For phase values around $\pi $, however, the $\Delta _b$ corrections can also give rise to a significant enhancement of the cross section. In particular, for $\phi _{M_3}$ the quantity $\Delta _b$ changes sign between $\phi _{M_3}=0$ and $\phi _{M_3}=\pi $, such that the bottom-Yukawa coupling is suppressed for small values of $\phi _{M_3}$ and enhanced for $\phi _{M_3}$ values close to $\pi $ as a consequence of the resummation of the $\Delta _b$ corrections. The reduction of the scale uncertainties from LO QCD to our best prediction for the cross section is similar as in the previous plots. The K-factors in the lower panel show that the dependence of the NLO cross sections on the phases $\phi _{A_t}$ and $\phi _{M_3}$ follows a similar trend as the LO cross section. In the plot on the right, the asymmetric K-factor dependence on $\phi _{M_3}$ is related to the direct dependence of $\Delta _b$ on the phase $\phi _{M_3}$.

In Fig. 9 we separately analyse the squark contributions for the LO cross section, i.e. the prediction omitting the squark loop contributions (black dot-dashed curves) is compared with the ones where first the pure LO squark contributions are added (depicted in cyan), and then the resummation of the $\Delta _b$ contributions to the bottom-quark Yukawa coupling is taken into account. For the latter both the results for the full ($\Delta _{b2}$, blue) and the simplified ($\Delta _{b1}$, red) resummation are shown. While the pure LO squark contributions are seen to have a moderate effect, it can be seen that the incorporation of the resummation of the $\Delta _b$ contribution leads to a significant enhancement of the squark loop effects. We furthermore confirm that for the heavy neutral Higgs bosons considered here the simplified resummation approximates the full resummation of the $\Delta _b$ contribution very well. The curves corresponding to $\Delta _{b2}$ and $\Delta _{b1}$ hardly differ from each other both for the variation of $\phi _{A_t}$ and $\phi _{M_3}$. As before all curves include the same ${\hat{\mathbf Z}}$ factors obtained from FeynHiggs. The results for $h_o$, which are not shown here, are qualitatively very similar. The LO squark contributions are less relevant for the $h_o$ cross section, since those contributions are absent in the MSSM with real parameters. We also note that the curves for $h_o$ follow a similar behaviour as the ones for $h_e$, which implies that there are no large cancellations expected in the sum of the cross sections for the two heavy Higgs bosons times their respective branching ratios. Thus, the phases entering $\Delta _b$ could potentially lead to observable effects in the production of the two heavy Higgs bosons even if the two states cannot be experimentally resolved as separate signals.

Having discussed the three different sources for $\mathcal {CP}$-violating effects relevant for Higgs boson production through gluon fusion in the MSSM – squark loop contributions, admixtures through ${\hat{\mathbf Z}}$ factors and resummation of $\Delta _b$ contributions – for completeness we also briefly discuss the bottom-quark annihilation cross section for the $m_h^{\mathrm {mod}+}$-inspired scenario with $\tan \beta =40$. The corresponding cross section is shown in Fig. 10 as a function of $\phi _{A_t}$ and $\phi _{M_3}$. For such a large value of $\tan \beta $ this cross section exceeds the gluon-fusion cross section by far. It shows a very significant dependence on the phases $\phi _{A_t}$ and $\phi _{M_3}$, which is mainly induced by the $\Delta _b$ contribution.

6 Remaining theoretical uncertainties

In the previous section we analysed our cross section predictions regarding $\mathcal {CP}$-violating effects entering via squark loop contributions, ${\hat{\mathbf Z}}$ factors and $\Delta _b$ contributions. Therein, we included renormalisation and factorisation scale uncertainties, which as expected are reduced upon inclusion of higher-order corrections. However, the cross section predictions are also affected by other relevant theoretical uncertainties, which we want to discuss in detail in this section.

Some of the theoretical uncertainties of cross sections in the MSSM with complex parameters are very similar to the ones in the MSSM with real parameters as discussed in Ref. [119]. Therefore, we can directly transfer the discussion of PDF$+\alpha _s$ uncertainties as well as the uncertainty associated with the renormalisation prescription for the bottom-quark Yukawa coupling from the case of the MSSM with real parameters:

PDF$+\alpha _s$ uncertainties: The fitted parton distribution functions (PDF) and the associated value of $\alpha _s$ induce an uncertainty in the prediction of the gluon-fusion cross section and, in particular, also the bottom-quark annihilation cross section. In our calculation we employ the MMHT2014 PDF sets at LO, NLO and NNLO [118], which can be used for both gluon fusion and bottom-quark annihilation. In Refs. [24, 119] it was observed that despite the effects of squarks in supersymmetric models, the PDF$+\alpha _s$ uncertainties are mostly a function of the Higgs boson mass $m_{h_a}$. We will therefore not discuss them in more detail, since – similar to the prescription for MSSM Higgs boson cross sections by the LHC Higgs Cross Section Working Group [6] – relative uncertainties can be taken over from tabulated relative uncertainties obtained for the SM Higgs boson or a pseudoscalar (in a 2HDM with $\tan \beta =1$) as a function of its mass. For Higgs masses in the range between 50 GeV and 1 TeV the typical size of PDF$+\alpha _s$ uncertainties for gluon fusion is ${\pm } (3{-}5)$% following the prescription of Ref. [132]. They increase up to ${\pm } 10$% for Higgs masses up to 2 TeV. For bottom-quark annihilation they are in the range $\pm (3{-}8)$% for Higgs masses between 50 GeV and 1 TeV and up to ${\pm } 16$% for Higgs masses below 2 TeV.
Renormalisation of the bottom-quark mass and definition of the bottom Yukawa coupling: In our calculation the bottom-quark mass is renormalised on-shell, and the bottom-Yukawa coupling is obtained from the bottom-quark mass as described in Sect. 2.1. The renormalisation of the bottom-quark mass and the freedom in the definition of the bottom-Yukawa coupling are known to have a sizeable numerical impact on the cross section predictions. This is in particular the case for large values of $\tan \beta $ where the bottom-Yukawa coupling of the heavy Higgs bosons is significantly enhanced and the top-quark Yukawa coupling is suppressed. On the other hand, in these regions of parameter space bottom-quark annihilation is the dominant process, for which there is less ambiguity regarding an appropriate choice for the renormalisation scale. The described uncertainties in the MSSM with complex parameters are analogous to the case of real parameters. We therefore refer to the discussion in Ref. [119] and references therein for further details.

We neglect approximate NNLO stop-quark contributions and accordingly the uncertainty associated with the approximation of the involved Wilson coefficients, which was discussed in Ref. [119]. The impact of the NNLO stop-quark contributions for the case of the MSSM with real parameters can be compared with our estimate for the renormalisation and factorisation scale uncertainty of our calculation. As an example, the NNLO stop-quark contributions lower the inclusive cross section for the light Higgs boson by about 2 pb for zero phases in the light-stop inspired scenario, which is at the lower edge of the scale uncertainty depicted in Fig. 3b. Other uncertainties discussed in Ref. [119] are renormalisation and factorisation scale uncertainties and an uncertainty related to higher-order contributions to $\Delta _b$. Moreover, we add another uncertainty related to the performed interpolation of supersymmetric NLO QCD contributions. We discuss in the following our estimates for the three previously mentioned uncertainties:

We obtain the renormalisation and factorisation scale uncertainty as follows: The central scale choice is $(\mu _\text {R}^0,\mu _\text {F}^0)=(m_{h_a}/2,m_{h_a}/2)$ for gluon fusion and $(\mu _\text {R}^0,\mu _\text {F}^0)=(m_{h_a},m_{h_a}/4)$ for bottom-quark annihilation. We obtain the scale uncertainty by taking the maximal deviation from the central scale choice $\Delta \sigma $ obtained from the additional scale choices $(\mu _\text {R},\mu _\text {F})\in \{(2\mu _\text {R}^0,2\mu _\text {F}^0),(2\mu _\text {R}^0,\mu _\text {F}^0),(\mu _\text {R}^0,2\mu _\text {F}^0), (\mu _\text {R}^0,\mu _\text {F}^0/2),(\mu _\text {R}^0/2,\mu _\text {F}^0),(\mu _\text {R}^0/2,\mu _\text {F}^0/2)\}$. We perform this procedure individually for all three cross sections in Eq. (41) and then obtain the overall absolute uncertainty through
$$\begin{aligned} \Delta \sigma ^{\mathrm {scale}}=\sqrt{\left( \Delta \sigma _{{{\mathrm{N}^{k}\mathrm{LO}}}}^{\Delta _{b1}}\right) ^2 +\left( \Delta \sigma _{{\mathrm{LO}{}}}^{\Delta _{b2}} -\Delta \sigma _{{\mathrm{LO}{}}}^{\Delta _{b1}}\right) ^2}, \end{aligned}$$
(45)
where we assume the two LO cross sections to be fully correlated. The uncertainty bands that we have displayed in the plots shown above correspond to the cross section range covered by $\sigma \pm \Delta \sigma ^{\mathrm {scale}}$.
In order to display the propagation of an uncertainty arising from higher-order contributions to $\Delta _b$ to our cross section calculation, we vary the value of $\Delta _b$ obtained from FeynHiggs by ${\pm } 10\%$. This variation by ${\pm } 10\%$ roughly corresponds to the effect of a variation of the renormalisation scales, see the discussion in Ref. [119]. We label the obtained uncertainty as $\Delta \sigma ^{\mathrm {resum}}$ and assign an uncertainty band of $\sigma \pm \Delta \sigma ^{\mathrm {resum}}$.
The employed interpolation for the two-loop virtual squark-gluino contributions following Eq. (35) leads to a further uncertainty. A conservative estimate for it can be obtained as follows: We determine the cross section $\sigma (\phi _z)$ following Eq. (38) not only for the correct phase $\phi _z$ in Eq. (35), but also leave the phase within Eq. (35) constant, i.e. fixed to 0 and $\pi $. We call the obtained cross sections $\sigma (0)$ and $\sigma (\pi )$. For each value of $\phi _z$ we take the difference $\Delta \sigma ^{\mathrm {int}}=\sin ^2(\phi _z)|\sigma (0)-\sigma (\pi )|/2$. It is reweighted with $\sin ^2(\phi _z)$, since we know that our result is correct at phases 0 and $\pi $. The obtained uncertainty band is given by $\sigma \pm \Delta \sigma ^{\mathrm {int}}$.

In the following we display the effects of the estimated uncertainties for certain scenarios, where we choose the displayed scenarios and the displayed cross sections such that the effect of the uncertainties is largest. While the scale uncertainties were included in all previous figures for the LO prediction as well as for our best prediction already, we will discuss the interpolation uncertainty for the light-stop inspired scenario with $\tan \beta =16$ and the resummation uncertainty for the $m_h^{\mathrm {mod}+}$-inspired scenario with $\tan \beta =40$.

Figure 11 shows the renormalisation and factorisation scale uncertainties $\Delta \sigma ^{\mathrm {scale}}$ as before and in addition the above described interpolation uncertainty $\Delta \sigma ^{\mathrm {int}}$, which in case of the variation of $\phi _{M_3}$ can be substantial. As can be seen in Fig. 11, the interpolation uncertainty obtained from our conservative estimate can in this scenario even exceed the scale uncertainty for the gluon-fusion cross section of $h_2$. It should be noted that this is an extreme case, while the interpolation uncertainty, which is an NLO effect related to the squark and gluino loop contributions, remains small for the other previously described scenarios (which we do not show here explicitly). This is simply a consequence of the fact that the relative impact of the squark and gluino contributions in the other scenarios is much smaller than in the light-stop inspired scenario. The interpolation uncertainty for the gluon-fusion cross section of $h_3$ in Fig. 11 is much less pronounced than for $h_2$, since as discussed above the squark loop corrections are significantly smaller in this case and would vanish if $h_3$ were a pure $\mathcal {CP}$-odd state. The behaviour in the lower panels of Fig. 11 displays the fact that by construction the assigned interpolation uncertainty vanishes for the phases 0 and $\pi $, where the interpolated result in the MSSM with complex parameters merges the known result of the MSSM with real parameters. For the variation of $\phi _{A_t}$ the LO cross section incorporating squark contributions already includes the dominant effect on the cross section, such that the uncertainty due to the interpolated NLO contributions is also less pronounced than in case of the variation of $\phi _{M_3}$.

The described $\Delta _b$ uncertainties are depicted in Fig. 12. Since $\Delta _b$ crosses 0 as a function of $\phi _{M_3}$ twice, the uncertainty that we have associated to it according to the prescription discussed above also vanishes there, as can be seen in the lower panel of Fig. 12b. Even for the large value of $\tan \beta $ chosen here the assigned $\Delta _b$ uncertainty of $\pm 10\%$ is much smaller than the scale uncertainty of the displayed cross sections. Despite the different behaviour with the phases $\phi _{A_t}$ and $\phi _{M_3}$ displayed in the lower panel of Fig. 12 the qualitative effect of the resummation uncertainties on the Higgs boson production cross sections is nevertheless rather similar. The latter is also true for the bottom-quark annihilation cross section, which is not depicted here. The resummation uncertainties are of most relevance for large values of $\tan \beta $, where the cross section of bottom-quark annihilation exceeds the gluon-fusion cross section.

7 Conclusions

In this paper we have presented theoretical predictions for inclusive cross sections for neutral Higgs boson production via gluon fusion and bottom-quark annihilation in the MSSM with complex parameters, and demonstrated the relevance of the $\mathcal {CP}$-violating phases on these cross sections.

The cross section predictions for the gluon-fusion process at leading-order are based on an explicit calculation taking into account the dependence on all complex parameters in the MSSM, and the complete form of the analytical formulae for the general $\mathcal {CP}$-violating case including Higgs mixing has been presented in the literature for the first time. The wave function normalisation factors arising from the $(3\times 3)$-mixing of the lowest-order mass eigenstates of the Higgs bosons $\lbrace h, H, A\rbrace $ into the loop-corrected mass eigenstates $\lbrace h_1, h_2, h_3\rbrace $ have been described with full propagator corrections using the self-energies of the neutral Higgs bosons as provided by FeynHiggs. Furthermore, the LO predictions for the gluon-fusion process in the MSSM with complex parameters deviate from those of the MSSM with real parameters due to non-zero couplings of the squarks to the pseudoscalar A and potentially different left- and right-handed bottom-Yukawa couplings arising from the resummation of $\tan \beta $-enhanced sbottom contributions in $\Delta _b$. We have supplemented the LO computation of the cross section by higher-order contributions: using for the treatment of the higher-order corrections a simplified version of the $\Delta _b$ resummation we have included the full massive top- and bottom quark contributions at NLO QCD and have interpolated the NLO SUSY QCD corrections from the amplitudes in the MSSM with real parameters. We have thoroughly discussed the uncertainties involved in using such an interpolation. The interpolation uncertainty at NLO, which is most relevant in scenarios where the squarks and the gluino are relatively light in view of the present limits from the LHC searches, could be avoided if an explicit result for the squark-gluino contributions at NLO QCD in the MSSM becomes available for the general case of complex parameters. For the top-quark contribution in the effective theory of a heavy top-quark we have added NNLO QCD contributions for all Higgs bosons, and N$^{3}$LO QCD contributions in an expansion around the threshold of Higgs production for the $\mathcal {CP}$-even component of the light Higgs boson $h_1$ to match the precision of the predictions for the SM Higgs boson. Electroweak effects, which include two-loop contributions with couplings of the heavy gauge bosons to the $\mathcal {CP}$-even component of the Higgs bosons mediated by light quarks, have been added to the $\mathcal {CP}$-even component of the gluon-fusion cross section.

The results presented in this paper are currently the state of the art for neutral Higgs production in the MSSM with complex parameters. Our calculations have been implemented in an extension of the code SusHi called SusHiMi, which is linked to FeynHiggs. SusHiMi is available upon request. Using SusHiMi, we have investigated the phenomenological effects of $\mathcal {CP}$-violating phases on the production of Higgs bosons in the MSSM with complex parameters in two slightly modified benchmark scenarios, light-stop and $m_h^{\text {mod+}}$. We have found in our analysis of Higgs boson production through gluon fusion that a proper description of squark and gluino loop contributions is essential. This refers both to the loop contributions to the gluon–gluon–Higgs vertex and to the corrections entering through $\Delta _b$. Squark and gluino loop contributions furthermore enter the wave function normalisation factors that are necessary to ensure the correct on-shell properties of the produced Higgs boson. Where squark and gluino contributions are sizeable the production cross sections show a significant dependence on the $\mathcal {CP}$-violating phases. We have discussed the remaining theoretical uncertainties in the cross section predictions taking into account renormalisation and factorisation scale uncertainties, a resummation uncertainty for $\Delta _b$ and an uncertainty due to the performed interpolation of NLO SUSY QCD corrections. We have furthermore briefly commented on other uncertainties that can directly be taken over from the case of the MSSM with real parameters.

A further important feature that occurs in the production processes for the two heavy states $h_2$ and $h_3$ in the general case where $\mathcal {CP}$-violating interactions are taken into account is the fact that there can be a large mixing between these often nearly mass-degenerate states. Their mixing effects are incorporated in the wave function normalisation factors for the external Higgs bosons. For a proper interpretation of experimental exclusion limits arising from MSSM Higgs searches, which so far have only been analysed in the framework of the $\mathcal {CP}$-conserving MSSM, it will be important to take into account interference effects in the full process of Higgs production and decay. Our results for the cross sections for on-shell Higgs bosons can be directly used in the context of a generalised narrow-width approximation to incorporate these interference effects. This topic will be addressed in a forthcoming publication.

Notes

Another approach to calculate the Higgs boson sector in the MSSM with complex parameters is based on the renormalisation group improved effective potential approach and implemented in e.g. the code CPsuperH [33, 34].
The soft-breaking parameter $M_1$ associated with the bino can also be complex, but has a minor impact on the Higgs sector, and we neglect its phase dependence in the following.
We note that the convention differs from the convention employed by FeynHiggs by a different sign of $\chi _1^0$ and $\phi _1^-$, which induces different signs in the corresponding elements of the matrices in Eqs. (11) and (12) and the $\chi _1^0$ couplings to (s)quarks displayed in the Appendix.
In the MSSM with real parameters only couplings involving $\tilde{f}_i-\tilde{f}_j-A$ with $i\ne j$ are non-vanishing, and left- and right-handed quark Yukawa couplings are identical, $g^\phi _q\equiv g^\phi _{q_L}=g^\phi _{q_R}$.
Most recently also N$^{3}$LO QCD corrections for $\mathcal {CP}$-odd Higgs bosons became available [84, 85]. We neglect those corrections in our analysis.
Exact numerical and for certain contributions analytical results for NLO virtual contributions were presented in Refs. [74, 75, 105,106,107].
The name is inspired by the mixing of the Higgs bosons. SusHiMi can be obtained upon request.
These formulas equal the master formulas employed in previous SusHi releases [21, 22].
Indirect bounds from the effects of stops on the measured Higgs rates are much weaker, see e.g. Ref. [117].

References

G. Aad et al. (ATLAS Collaboration), Phys. Lett. B 716, 1 (2012). arXiv:1207.7214
S. Chatrchyan et al. (CMS Collaboration), Phys. Lett. B 716, 30 (2012). arXiv:1207.7235
S. Dittmaier et al., Handbook of LHC Higgs cross sections: 1. Inclusive observables. arXiv:1101.0593
S. Dittmaier et al., Handbook of LHC Higgs cross sections: 2. Differential distributions. arXiv:1201.3084
S. Heinemeyer et al., Handbook of LHC Higgs cross sections: 3. Higgs properties. arXiv:1307.1347
D. de Florian et al., Handbook of LHC Higgs cross sections: 4. Deciphering the nature of the Higgs sector. arXiv:1610.07922
E. Fuchs, S. Thewes, G. Weiglein, Eur. Phys. J. C 75, 254 (2015). arXiv:1411.4652
ADS Google Scholar
N. Cabibbo, Phys. Rev. Lett. 10, 531 (1963)
ADS Google Scholar
M. Kobayashi, T. Maskawa, Prog. Theor. Phys. 49, 652 (1973)
ADS Google Scholar
A. Riotto, arXiv:hep-ph/9807454
J.M. Cline, arXiv:hep-ph/0609145
S. Davidson, E. Nardi, Y. Nir, Phys. Rep. 466, 105 (2008). arXiv:0802.2962
ADS Google Scholar
A. Dedes, S. Moretti, Nucl. Phys. B 576, 29 (2000). arXiv:hep-ph/9909418
ADS Google Scholar
S.Y. Choi, J.S. Lee, Phys. Rev. D 61, 115002 (2000). arXiv:hep-ph/9910557
ADS Google Scholar
S.Y. Choi, K. Hagiwara, J.S. Lee, Phys. Lett. B 529, 212 (2002). arXiv:hep-ph/0110138
ADS Google Scholar
M. Carena, J.R. Ellis, S. Mrenna, A. Pilaftsis, C.E.M. Wagner, Nucl. Phys. B 659, 145 (2003). arXiv:hep-ph/0211467
ADS Google Scholar
A. Arhrib, D.K. Ghosh, O.C.W. Kong, Phys. Lett. B 537, 217 (2002). arXiv:hep-ph/0112039
ADS Google Scholar
Q.H. Cao, D. Nomura, K. Tobe, C.-P. Yuan, Phys. Lett. B 632, 688 (2006). arXiv:hep-ph/0508311
ADS Google Scholar
S. Hesselbach, S. Moretti, S. Munir, P. Poulose, Eur. Phys. J. C 54, 129 (2008). arXiv:0706.4269
ADS Google Scholar
S. Hesselbach, S. Moretti, S. Munir, P. Poulose, Phys. Rev. D 82, 074004 (2010). arXiv:0903.0747
ADS Google Scholar
R.V. Harlander, S. Liebler, H. Mantler, Comput. Phys. Commun. 184, 1605 (2013). arXiv:1212.3249
ADS Google Scholar
R.V. Harlander, S. Liebler and H. Mantler, Comput. Phys. Commun. 212, 239 (2017). doi:10.1016/j.cpc.2016.10.015. arXiv:1605.03190 [hep-ph]
R.V. Harlander, Eur. Phys. J. C 76(5), 252 (2016). arXiv:1512.04901
ADS Google Scholar
S. Liebler, Eur. Phys. J. C 75(5), 210 (2015). arXiv:1502.07972
ADS Google Scholar
S. Heinemeyer, W. Hollik, G. Weiglein, Comput. Phys. Commun. 124, 76 (2000). arXiv:hep-ph/9812320
ADS Google Scholar
S. Heinemeyer, W. Hollik, G. Weiglein, Eur. Phys. J. C 9, 343 (1999). arXiv:hep-ph/9812472
ADS Google Scholar
G. Degrassi, S. Heinemeyer, W. Hollik, P. Slavich, G. Weiglein, Eur. Phys. J. C 28, 133 (2003). arXiv:hep-ph/0212020
ADS Google Scholar
M. Frank, T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, JHEP 0702, 047 (2007). arXiv:hep-ph/0611326
ADS Google Scholar
T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Phys. Rev. Lett. 112(14), 141801 (2014). arXiv:1312.4937
ADS Google Scholar
S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Phys. Lett. B 652, 300 (2007). arXiv:0705.0746
ADS Google Scholar
W. Hollik, S. Paßehr, Phys. Lett. B 733, 144 (2014). arXiv:1401.8275
ADS Google Scholar
W. Hollik, S. Paßehr, JHEP 1410, 171 (2014). arXiv:1409.1687
ADS Google Scholar
J.S. Lee, A. Pilaftsis, M. Carena, S.Y. Choi, M. Drees, J.R. Ellis, C.E.M. Wagner, Comput. Phys. Commun. 156, 283 (2004). arXiv:hep-ph/0307377
ADS Google Scholar
J.S. Lee, M. Carena, J. Ellis, A. Pilaftsis, C.E.M. Wagner, Comput. Phys. Commun. 180, 312 (2009). arXiv:0712.2360
ADS Google Scholar
K.E. Williams, G. Weiglein, Phys. Lett. B 660, 217 (2008). arXiv:0710.5320
ADS Google Scholar
A.C. Fowler, G. Weiglein, JHEP 1001, 108 (2010). arXiv:0909.5165
ADS Google Scholar
K.E. Williams, H. Rzehak, G. Weiglein, Eur. Phys. J. C 71, 1669 (2011). arXiv:1103.1335
ADS Google Scholar
E. Fuchs, DESY-THESIS-2015-037, Ph.D. thesis (2015)
E. Fuchs, G. Weiglein, arXiv:1610.06193
N. Kauer, Phys. Lett. B 649, 413 (2007). arXiv:hep-ph/0703077
ADS Google Scholar
D. Berdine, N. Kauer, D. Rainwater, Phys. Rev. Lett. 99, 111601 (2007). arXiv:hep-ph/0703058
ADS Google Scholar
C.F. Uhlemann, N. Kauer, Nucl. Phys. B 814, 195 (2009). arXiv:0807.4112
ADS Google Scholar
G. Cacciapaglia, A. Deandrea, S. De Curtis, Phys. Lett. B 682, 43 (2009). arXiv:0906.3417
ADS Google Scholar
J. Reuter, eConf C 0705302, SUS14 (2007). arXiv:0709.0068
A.C. Fowler, Ph.D. thesis. http://etheses.dur.ac.uk/449/1/Thesis_AFowler.pdf (2010)
D. Barducci, A. Belyaev, J. Blamey, S. Moretti, L. Panizzi, H. Prager, JHEP 1407, 142 (2014). arXiv:1311.3977
ADS Google Scholar
A. Djouadi, P. Gambino, S. Heinemeyer, W. Hollik, C. Junger, G. Weiglein, Phys. Rev. D 57, 4179 (1998). arXiv:hep-ph/9710438
ADS Google Scholar
T. Banks, Nucl. Phys. B 303, 172 (1988)
ADS Google Scholar
L.J. Hall, R. Rattazzi, U. Sarid, Phys. Rev. D 50, 7048 (1994). arXiv:hep-ph/9306309
ADS Google Scholar
R. Hempfling, Phys. Rev. D 49, 6168 (1994)
ADS Google Scholar
M.S. Carena, M. Olechowski, S. Pokorski, C.E.M. Wagner, Nucl. Phys. B 426, 269 (1994). arXiv:hep-ph/9402253
ADS Google Scholar
M.S. Carena, D. Garcia, U. Nierste, C.E.M. Wagner, Nucl. Phys. B 577, 88 (2000). arXiv:hep-ph/9912516
ADS Google Scholar
M. Carena, D. Garcia, U. Nierste, C.E.M. Wagner, Phys. Lett. B 499, 141 (2001). arXiv:hep-ph/0010003
ADS Google Scholar
J. Baglio, R. Gröber, M. Mühlleitner, D.T. Nhung, H. Rzehak, M. Spira, J. Streicher, K. Walz, Comput. Phys. Commun. 185(12), 3372 (2014). arXiv:1312.4788
ADS Google Scholar
M. Frank, L. Galeta, T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Phys. Rev. D 88(5), 055013 (2013). arXiv:1306.1156
ADS Google Scholar
L. Hofer, U. Nierste, D. Scherer, JHEP 0910, 081 (2009). arXiv:0907.5408
ADS Google Scholar
D. Noth, M. Spira, JHEP 1106, 084 (2011). arXiv:1001.1935
ADS Google Scholar
D. Noth, M. Spira, Phys. Rev. Lett. 101, 181801 (2008). arXiv:0808.0087
ADS Google Scholar
A. Bauer, L. Mihaila, J. Salomon, JHEP 0902, 037 (2009). arXiv:0810.5101
ADS Google Scholar
R.D. Peccei, H.R. Quinn, Phys. Rev. D 16, 1791 (1977)
ADS Google Scholar
R.D. Peccei, H.R. Quinn, Phys. Rev. Lett. 38, 1440 (1977)
ADS Google Scholar
A. Bharucha, A. Fowler, G. Moortgat-Pick, G. Weiglein, JHEP 1305, 053 (2013). arXiv:1211.3134
ADS Google Scholar
O. Brein, A. Djouadi, R. Harlander, Phys. Lett. B 579, 149 (2004). arXiv:hep-ph/0307206
ADS Google Scholar
R.V. Harlander, S. Liebler, T. Zirke, JHEP 1402, 023 (2014). arXiv:1307.8122
ADS Google Scholar
P.H. Chankowski, S. Pokorski, J. Rosiek, Nucl. Phys. B 423, 437 (1994). arXiv:hep-ph/9303309
ADS Google Scholar
A. Dabelstein, Nucl. Phys. B 456, 25 (1995). arXiv:hep-ph/9503443
ADS Google Scholar
S. Heinemeyer, W. Hollik, J. Rosiek, G. Weiglein, Eur. Phys. J. C 19, 535 (2001). arXiv:hep-ph/0102081
ADS Google Scholar
H. Lehmann, K. Symanzik, W. Zimmermann, Nuovo Cim. 1, 205 (1955)
Google Scholar
H.M. Georgi, S.L. Glashow, M.E. Machacek, D.V. Nanopoulos, Phys. Rev. Lett. 40, 692 (1978)
ADS Google Scholar
A. Djouadi, M. Spira, P.M. Zerwas, Phys. Lett. B 264, 440 (1991)
ADS Google Scholar
S. Dawson, Nucl. Phys. B 359, 283 (1991)
ADS Google Scholar
M. Spira, A. Djouadi, D. Graudenz, P.M. Zerwas, Nucl. Phys. B 453, 17 (1995). arXiv:hep-ph/9504378
ADS Google Scholar
R. Harlander, P. Kant, JHEP 0512, 015 (2005). arXiv:hep-ph/0509189
ADS Google Scholar
C. Anastasiou, S. Beerli, S. Bucherer, A. Daleo, Z. Kunszt, JHEP 0701, 082 (2007). arXiv:hep-ph/0611236
ADS Google Scholar
U. Aglietti, R. Bonciani, G. Degrassi, A. Vicini, JHEP 0701, 021 (2007). arXiv:hep-ph/0611266
ADS Google Scholar
R.V. Harlander, W.B. Kilgore, Phys. Rev. Lett. 88, 201801 (2002). arXiv:hep-ph/0201206
ADS Google Scholar
C. Anastasiou, K. Melnikov, Nucl. Phys. B 646, 220 (2002). arXiv:hep-ph/0207004
ADS Google Scholar
V. Ravindran, J. Smith, W.L. van Neerven, Nucl. Phys. B 665, 325 (2003). arXiv:hep-ph/0302135
ADS Google Scholar
C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog, B. Mistlberger, Phys. Lett. B 737, 325 (2014). arXiv:1403.4616
ADS Google Scholar
Y. Li, A. von Manteuffel, R.M. Schabinger, H.X. Zhu, Phys. Rev. D 91, 036008 (2015). arXiv:1412.2771
ADS Google Scholar
C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog, B. Mistlberger, JHEP 1503, 091 (2015). arXiv:1411.3584
Google Scholar
C. Anastasiou, C. Duhr, F. Dulat, F. Herzog, B. Mistlberger, Phys. Rev. Lett. 114, 212001 (2015). arXiv:1503.06056
ADS Google Scholar
C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog, A. Lazopoulos, B. Mistlberger, JHEP 1605, 058 (2016). arXiv:1602.00695
ADS Google Scholar
T. Ahmed, M.C. Kumar, P. Mathews, N. Rana, V. Ravindran, Eur. Phys. J. C 76, 355 (2016). arXiv:1510.02235
ADS Google Scholar
T. Ahmed, M. Bonvini, M.C. Kumar, P. Mathews, N. Rana, V. Ravindran, L. Rottoli, Eur. Phys. J. C 76(12), 663 (2016). doi:10.1140/epjc/s10052-016-4510-1. arXiv:1606.00837 [hep-ph]
R.V. Harlander, H. Mantler, S. Marzani, K.J. Ozeren, Eur. Phys. J. C 66, 359 (2010). arXiv:0912.2104
ADS Google Scholar
R.V. Harlander, K.J. Ozeren, JHEP 0911, 088 (2009). arXiv:0909.3420
ADS Google Scholar
R.V. Harlander, K.J. Ozeren, Phys. Lett. B 679, 467 (2009). arXiv:0907.2997
ADS Google Scholar
A. Pak, M. Rogal, M. Steinhauser, JHEP 1002, 025 (2010). arXiv:0911.4662
ADS Google Scholar
A. Pak, M. Rogal, M. Steinhauser, Phys. Lett. B 679, 473 (2009). arXiv:0907.2998
ADS Google Scholar
A. Pak, M. Rogal, M. Steinhauser, JHEP 1109, 088 (2011). arXiv:1107.3391
ADS Google Scholar
S. Marzani, R.D. Ball, V. Del Duca, S. Forte, A. Vicini, Nucl. Phys. B 800, 127 (2008). arXiv:0801.2544
ADS Google Scholar
T. Neumann, C. Williams, Phys. Rev. D 95(1), 014004 (2017). doi:10.1103/PhysRevD.95.014004. arXiv:1609.00367 [hep-ph]
S. Actis, G. Passarino, C. Sturm, S. Uccirati, Phys. Lett. B 670, 12 (2008). arXiv:0809.1301
ADS Google Scholar
U. Aglietti, R. Bonciani, G. Degrassi, A. Vicini, Phys. Lett. B 595, 432 (2004). arXiv:hep-ph/0404071
ADS Google Scholar
R. Bonciani, G. Degrassi, A. Vicini, Comput. Phys. Commun. 182, 1253 (2011). arXiv:1007.1891
ADS Google Scholar
E. Bagnaschi, G. Degrassi, P. Slavich, A. Vicini, JHEP 1202, 088 (2012). arXiv:1111.2854
ADS Google Scholar
R.V. Harlander, M. Steinhauser, Phys. Lett. B 574, 258 (2003). arXiv:hep-ph/0307346
ADS Google Scholar
R.V. Harlander, M. Steinhauser, JHEP 0409, 066 (2004). arXiv:hep-ph/0409010
ADS Google Scholar
R.V. Harlander, F. Hofmann, JHEP 0603, 050 (2006). arXiv:hep-ph/0507041
ADS Google Scholar
G. Degrassi, P. Slavich, Nucl. Phys. B 805, 267 (2008). arXiv:0806.1495
ADS Google Scholar
G. Degrassi, P. Slavich, JHEP 1011, 044 (2010). arXiv:1007.3465
ADS Google Scholar
G. Degrassi, S. Di Vita, P. Slavich, JHEP 1108, 128 (2011). arXiv:1107.0914
ADS Google Scholar
G. Degrassi, S. Di Vita, P. Slavich, Eur. Phys. J. C 72, 2032 (2012). arXiv:1204.1016
ADS Google Scholar
C. Anastasiou, S. Beerli, A. Daleo, Phys. Rev. Lett. 100, 241806 (2008). arXiv:0803.3065
ADS Google Scholar
M. Mühlleitner, H. Rzehak, M. Spira, PoS RADCOR 2009, 043 (2010). arXiv:1001.3214
Google Scholar
M. Mühlleitner, M. Spira, Nucl. Phys. B 790, 1 (2008). arXiv:hep-ph/0612254
ADS Google Scholar
A. Pak, M. Steinhauser, N. Zerf, Eur. Phys. J. C 71, 1602 (2011). arXiv:1012.0639 [Eur. Phys. J. C 72, 2182 (2012)]
A. Pak, M. Steinhauser, N. Zerf, JHEP 1209, 118 (2012). arXiv:1208.1588
ADS Google Scholar
R. Harlander, M. Steinhauser, Phys. Rev. D 68, 111701 (2003). arXiv:hep-ph/0308210
ADS Google Scholar
S. Heinemeyer, W. Hollik, D. Stockinger, A.M. Weber, G. Weiglein, JHEP 0608, 052 (2006). arXiv:hep-ph/0604147
ADS Google Scholar
T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, arXiv:0710.4891
T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Comput. Phys. Commun. 180, 1426 (2009)
ADS Google Scholar
E. Bagnaschi, R.V. Harlander, H. Mantler, A. Vicini, M. Wiesemann, JHEP 1601, 090 (2016). arXiv:1510.08850
ADS Google Scholar
S. Liebler, H. Mantler, M. Wiesemann, arXiv:1608.02949
M. Carena, S. Heinemeyer, O. Stl, C.E.M. Wagner, G. Weiglein, Eur. Phys. J. C 73(9), 2552 (2013). arXiv:1302.7033
ADS Google Scholar
S. Liebler, S. Profumo, T. Stefaniak, JHEP 1604, 143 (2016). arXiv:1512.09172
ADS Google Scholar
L.A. Harland-Lang, A.D. Martin, P. Motylinski, R.S. Thorne, Eur. Phys. J. C 75(5), 204 (2015). arXiv:1412.3989
ADS Google Scholar
E. Bagnaschi, R.V. Harlander, S. Liebler, H. Mantler, P. Slavich, A. Vicini, JHEP 1406, 167 (2014). arXiv:1404.0327
ADS Google Scholar
D.A. Demir, O. Lebedev, K.A. Olive, M. Pospelov, A. Ritz, Nucl. Phys. B 680, 339 (2004). arXiv:hep-ph/0311314
ADS Google Scholar
D. Chang, W.Y. Keung, A. Pilaftsis, Phys. Rev. Lett. 82, 900 (1999). arXiv:hep-ph/9811202 [Erratum: Phys. Rev. Lett. 83, 3972 (1999)]
A. Pilaftsis, Phys. Lett. B 471, 174 (1999). arXiv:hep-ph/9909485
ADS Google Scholar
W. Hollik, J.I. Illana, S. Rigolin, D. Stöckinger, Phys. Lett. B 416, 345 (1998). arXiv:hep-ph/9707437
ADS Google Scholar
W. Hollik, J.I. Illana, S. Rigolin, D. Stöckinger, Phys. Lett. B 425, 322 (1998). arXiv:hep-ph/9711322
ADS Google Scholar
O. Lebedev, K.A. Olive, M. Pospelov, A. Ritz, Phys. Rev. D 70, 016003 (2004). arXiv:hep-ph/0402023
ADS Google Scholar
V.D. Barger, T. Falk, T. Han, J. Jiang, T. Li, T. Plehn, Phys. Rev. D 64, 056007 (2001). arXiv:hep-ph/0101106
ADS Google Scholar
Y. Li, S. Profumo, M. Ramsey-Musolf, JHEP 1008, 062 (2010). arXiv:1006.1440
ADS Google Scholar
Y. Nakai, M. Reece. arXiv:1612.08090
G.F. Giudice, A. Romanino, Phys. Lett. B 634, 307 (2006). arXiv:hep-ph/0510197
ADS Google Scholar
M. Carena, J. Ellis, J.S. Lee, A. Pilaftsis, C.E.M. Wagner, JHEP 1602, 123 (2016). arXiv:1512.00437
ADS Google Scholar
S. Forte et al., LHCHXSWG-DRAFT-INT-2016-011
J. Butterworth et al., J. Phys. G 43, 023001 (2016). arXiv:1510.03865
ADS Google Scholar
S. Liebler, W. Porod, Nucl. Phys. B 849, 213 (2011). arXiv:1011.6163 [Erratum: Nucl. Phys. B 856, 125 (2012)]
S. Liebler, Ph.D. thesis (2011)

Download references

Acknowledgements

We thank Sebastian Paßehr and Pietro Slavich for discussions and Elina Fuchs for discussions and comments on the manuscript. The authors acknowledge support by Deutsche Forschungsgemeinschaft through the SFB 676 “Particles, Strings and the Early Universe” and by the European Commission through the “HiggsTools” Initial Training Network PITN-GA-2012-316704.

Author information

Authors and Affiliations

DESY, Notkestraße 85, 22607, Hamburg, Germany
Stefan Liebler, Shruti Patel & Georg Weiglein

Authors

Stefan Liebler
View author publications
You can also search for this author in PubMed Google Scholar
Shruti Patel
View author publications
You can also search for this author in PubMed Google Scholar
Georg Weiglein
View author publications
You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Shruti Patel.

Appendix A: Formulas: Higgs-quark and Higgs-squark couplings

In SusHiMi the Higgs-(s)quark couplings are expressed in terms of the $\mathcal {CP}$-even and $\mathcal {CP}$-odd neutral gauge eigenstates $\phi _g\in \lbrace \phi _1^0,\phi _2^0\rbrace $ and $\chi _g\in \lbrace \chi _1^0,\chi _2^0\rbrace $, respectively. In order to obtain the couplings of the squarks with the lowest-order mass eigenstates $\phi \in \lbrace h,H,A,G\rbrace $ the gauge eigenstates are rotated using the tree-level mixing matrix ${\mathcal {R}}$ as depicted in the following Feynman diagrams

(46)

with $v = 2m_W/g = 1/ \sqrt{\sqrt{2}G_F}$ and the tree-level mixing matrix ${\mathcal {R}}$ given in Eq. (11). At the amplitude level the results will then also be multiplied with the corresponding ${\hat{\mathbf Z}}$ factor.

The couplings between the gauge eigenstates and the third generation quarks are $g_{q_L}=1/\cos \beta $ for $\phi _1^0$ and $\chi _1^0$ and $g_{q_L}=1/\sin \beta $ for $\phi _2^0$ and $\chi _2^0$. For $\Delta _b$ corrections we refer to Eqs. (5) and (7).

The couplings between the gauge eigenstates and the third generation squarks contain terms from the squark mass diagonalisation matrix $U_{\tilde{q}}$ (see Eq. (3)) which is a unitary matrix with real diagonal elements and complex off-diagonal elements, i.e. it can be written as follows

$$\begin{aligned} U_{\tilde{q}} = \left( \begin{array}{c@{\quad }c} U_{\tilde{q}11} &{} U_{\tilde{q}12} \\ -U_{\tilde{q}12}^* &{} U_{\tilde{q}22} \end{array}\right) . \end{aligned}$$

(47)

They are obtained with MaCoR [133, 134]. For the $\mathcal {CP}$-even state $\phi _1^0$ we have the stop couplings (using $s_{\beta } \equiv \sin {\beta }$, $c_{\beta } \equiv \cos {\beta }$, $t_{\beta } \equiv \tan {\beta }, s_W \equiv \sin \theta _{W}$ and $c_W \equiv \cos \theta _W$):

$$\begin{aligned} g_{\tilde{t},11}^{\phi _1^0}= & {} U_{\tilde{t}12}^* \left[ -\frac{U_{\tilde{t}11}m_t \mu }{s_\beta } + \frac{4}{3} U_{\tilde{t}12} c_\beta m_Z^2 s^2_W \right] \nonumber \\&-U_{\tilde{t}11}^* \left[ \frac{U_{\tilde{t}12}m_t \mu ^*}{s_\beta } + c_\beta m_Z^2 U_{\tilde{t}11} \left( \frac{1}{3} s^2_W - c_W^2 \right) \right] \nonumber \\ g_{\tilde{t},12}^{\phi _1^0}= & {} U_{\tilde{t}12}^* \left[ -\frac{U_{\tilde{t}21}m_t \mu }{s_\beta } + \frac{4}{3} U_{\tilde{t}22} c_\beta m_Z^2 s^2_W\right] \nonumber \\&-U_{\tilde{t}11}^* \left[ \frac{U_{\tilde{t}22}m_t \mu ^*}{s_\beta } + c_\beta m_Z^2 U_{\tilde{t}21} \left( \frac{1}{3} s^2_W - c_W^2 \right) \right] \nonumber \\ g_{\tilde{t},21}^{\phi _1^0}= & {} U_{\tilde{t}22}^* \left[ -\frac{U_{\tilde{t}11}m_t \mu }{s_\beta } + \frac{4}{3} U_{\tilde{t}12} c_\beta m_Z^2 s^2_W \right] \nonumber \\&- U_{\tilde{t}21}^* \left[ \frac{U_{\tilde{t}12}m_t \mu ^*}{s_\beta } + c_\beta m_Z^2 U_{\tilde{t}11} \left( \frac{1}{3} s^2_W - c_W^2 \right) \right] \nonumber \\ g_{\tilde{t},22}^{\phi _1^0}= & {} U_{\tilde{t}22}^* \left[ -\frac{U_{\tilde{t}21}m_t \mu }{s_\beta } + \frac{4}{3} U_{\tilde{t}22} c_\beta m_Z^2 s^2_W \right] \nonumber \\&- U_{\tilde{t}21}^* \left[ \frac{U_{\tilde{t}22}m_t \mu ^*}{s_\beta } + c_\beta m_Z^2 U_{\tilde{t}21} \left( \frac{1}{3} s^2_W - c_W^2 \right) \right] .\nonumber \\ \end{aligned}$$

(48)

For the $\mathcal {CP}$-even state $\phi _2^0$, the couplings are:

$$\begin{aligned} g_{\tilde{t},11}^{\phi _2^0}= & {} U_{\tilde{t}12}^* \left[ \frac{U_{\tilde{t}11}m_tA_t^*}{s_\beta } + U_{\tilde{t}12} \left( \frac{2 m^2_t}{s_\beta } -\frac{4}{3} s_\beta m_z^2s^2_W \right) \right] \nonumber \\&+\, U_{\tilde{t}11}^* \left[ \frac{U_{\tilde{t}12}m_t A_t}{s_\beta } + \frac{2U_{\tilde{t}11}m^2_t}{s_\beta }\right. \nonumber \\&\left. +\, s_\beta m_Z^2 U_{\tilde{t}11} \left( \frac{1}{3} s^2_W - c^2_W \right) \right] \nonumber \\ g_{\tilde{t},12}^{\phi _2^0}= & {} U_{\tilde{t}12}^* \left[ \frac{U_{\tilde{t}21}m_tA_t^*}{s_\beta } + U_{\tilde{t}22} \left( \frac{2 m^2_t}{s_\beta } -\frac{4}{3} s_\beta m_z^2s^2_W \right) \right] \nonumber \\&+\, U_{\tilde{t}11}^* \left[ \frac{U_{\tilde{t}22}m_t A_t}{s_\beta } + \frac{2U_{\tilde{t}21}m^2_t}{s_\beta }\right. \nonumber \\&\left. +\, s_\beta m_Z^2 U_{\tilde{t}21} \left( \frac{1}{3} s^2_W - c^2_W \right) \right] \nonumber \\ g_{\tilde{t},21}^{\phi _2^0}= & {} U_{\tilde{t}22}^* \left[ \frac{U_{\tilde{t}11}m_tA_t^*}{s_\beta } + U_{\tilde{t}12} \left( \frac{2 m^2_t}{s_\beta } -\frac{4}{3} s_\beta m_z^2s^2_W \right) \right] \nonumber \\&+\, U_{\tilde{t}21}^* \left[ \frac{U_{\tilde{t}12}m_t A_t}{s_\beta } + \frac{2U_{\tilde{t}11}m^2_t}{s_\beta }\right. \nonumber \\&\left. +\, s_\beta m_Z^2 U_{\tilde{t}11} \left( \frac{1}{3} s^2_W - c^2_W \right) \right] \nonumber \\ g_{\tilde{t},22}^{\phi _2^0}= & {} U_{\tilde{t}22}^* \left[ \frac{U_{\tilde{t}21}m_tA_t^*}{s_\beta } + U_{\tilde{t}22} \left( \frac{2 m^2_t}{s_\beta } -\frac{4}{3} s_\beta m_z^2s^2_W \right) \right] \nonumber \\&+\, U_{\tilde{t}21}^* \left[ \frac{U_{\tilde{t}22}m_t A_t}{s_\beta } + \frac{2U_{\tilde{t}21}m^2_t}{s_\beta }\right. \nonumber \\&\left. +\, s_\beta m_Z^2 U_{\tilde{t}21} \left( \frac{1}{3} s^2_W - c^2_W \right) \right] . \end{aligned}$$

(49)

Similarly for the $\mathcal {CP}$-odd states $\chi _1^0$ and $\chi _2^0$ the stop couplings are given as:

$$\begin{aligned} g_{\tilde{t},11}^{\chi _1^0}= & {} i\frac{m_t}{s_\beta } \left[ -\mu U_{\tilde{t}12}^* U_{\tilde{t}11} + \mu ^* U_{\tilde{t}11}^* U_{\tilde{t}12} \right] \nonumber \\ g_{\tilde{t},11}^{\chi _2^0}= & {} i\frac{m_t}{s_\beta } \left[ -A_t^* U_{\tilde{t}12}^* U_{\tilde{t}11} + A_t U_{\tilde{t}11}^* U_{\tilde{t}12} \right] \nonumber \\ g_{\tilde{t},12}^{\chi _1^0}= & {} i\frac{m_t}{s_\beta } \left[ -\mu U_{\tilde{t}12}^* U_{\tilde{t}21} + \mu ^* U_{\tilde{t}11}^* U_{\tilde{t}22} \right] \nonumber \\ g_{\tilde{t},12}^{\chi _2^0}= & {} i\frac{m_t}{s_\beta } \left[ -A_t^* U_{\tilde{t}12}^* U_{\tilde{t}21} + A_t U_{\tilde{t}11}^* U_{\tilde{t}22} \right] \nonumber \\ g_{\tilde{t},21}^{\chi _1^0}= & {} i\frac{m_t}{s_\beta } \left[ -\mu U_{\tilde{t}22}^* U_{\tilde{t}11} + \mu ^* U_{\tilde{t}21}^* U_{\tilde{t}12} \right] \nonumber \\ g_{\tilde{t},21}^{\chi _2^0}= & {} i\frac{m_t}{s_\beta } \left[ -A_t^* U_{\tilde{t}22}^* U_{\tilde{t}11} + A_t U_{\tilde{t}21}^* U_{\tilde{t}12} \right] \nonumber \\ g_{\tilde{t},22}^{\chi _1^0}= & {} i\frac{m_t}{s_\beta } \left[ -\mu U_{\tilde{t}22}^* U_{\tilde{t}21} + \mu ^* U_{\tilde{t}21}^* U_{\tilde{t}22}\right] \nonumber \\ g_{\tilde{t},22}^{\chi _2^0}= & {} i\frac{m_t}{s_\beta } \left[ -A_t^* U_{\tilde{t}22}^* U_{\tilde{t}21} + A_t U_{\tilde{t}21}^* U_{\tilde{t}22} \right] . \end{aligned}$$

(50)

Analogously, the Higgs-sbottom couplings for the $\mathcal {CP}$-even state $\phi _1^0$ are:

$$\begin{aligned} g_{\tilde{b},11}^{\phi _1^0}= & {} U_{\tilde{b}12}^* \left[ \frac{U_{\tilde{b}11}A_b^* m_b}{c_\beta } + U_{\tilde{b}12} \left( \frac{2m_b^2}{c_\beta } - \frac{2}{3} c_\beta m_Z^2 s_W^2 \right) \right] \nonumber \\&+\, U_{\tilde{b}11}^* \left[ \frac{U_{\tilde{b}12}A_b m_b}{c_\beta } + \frac{2U_{\tilde{b}11}m_b^2}{c_\beta }\right. \nonumber \\&\left. -\, c_\beta m_Z^2U_{\tilde{b}11} \left( \frac{1}{3}s^2_W+ c_W^2 \right) \right] \nonumber \\ g_{\tilde{b},12}^{\phi _1^0}= & {} U_{\tilde{b}12}^* \left[ \frac{U_{\tilde{b}21}A_b^* m_b}{c_\beta } + U_{\tilde{b}22} \left( \frac{2m_b^2}{c_\beta } - \frac{2}{3} c_\beta m_Z^2 s_W^2 \right) \right] \nonumber \\&+\, U_{\tilde{b}11}^* \left[ \frac{U_{\tilde{b}22}A_b m_b}{c_\beta } + \frac{2U_{\tilde{b}21}m_b^2}{c_\beta }\right. \nonumber \\&\left. -\, c_\beta m_Z^2U_{\tilde{b}21} \left( \frac{1}{3}s^2_W + c_W^2 \right) \right] ,\nonumber \\ g_{\tilde{b},21}^{\phi _1^0}= & {} U_{\tilde{b}22}^* \left[ \frac{U_{\tilde{b}11}A_b^* m_b}{c_\beta } + U_{\tilde{b}12} \left( \frac{2m_b^2}{c_\beta } - \frac{2}{3} c_\beta m_Z^2 s_W^2 \right) \right] \nonumber \\&+\, U_{\tilde{b}21}^* \left[ \frac{U_{\tilde{b}12}A_b m_b}{c_\beta } + \frac{2U_{\tilde{b}11}m_b^2}{c_\beta }\right. \nonumber \\&\left. -\, c_\beta m_Z^2U_{\tilde{b}11} \left( \frac{1}{3}s^2_W + c_W^2 \right) \right] \nonumber \\ g_{\tilde{b},22}^{\phi _1^0}= & {} U_{\tilde{b}22}^* \left[ \frac{U_{\tilde{b}21}A_b^* m_b}{c_\beta } + U_{\tilde{b}22} \left( \frac{2m_b^2}{c_\beta } - \frac{2}{3} c_\beta m_Z^2 s_W^2 \right) \right] \nonumber \\&+\, U_{\tilde{b}21}^* \left[ \frac{U_{\tilde{b}22}A_b m_b}{c_\beta } + \frac{2U_{\tilde{b}21}m_b^2}{c_\beta }\right. \nonumber \\&\left. -\, c_\beta m_Z^2U_{\tilde{b}21} \left( \frac{1}{3}s^2_W + c_W^2 \right) \right] . \end{aligned}$$

(51)

For the $\mathcal {CP}$-even state $\phi _0^2$ they are given as:

$$\begin{aligned} g_{\tilde{b},11}^{\phi _2^0}= & {} U_{\tilde{b}12}^* \left[ - \frac{U_{\tilde{b}11} m_b \mu }{c_\beta } + \frac{2}{3} U_{\tilde{b}12} s_\beta m_z^2s_W^2 \right] \nonumber \\&+\, U_{\tilde{b}11}^* \left[ -\frac{U_{\tilde{b}12}m_b \mu ^*}{c_\beta } + s_\beta m_Z^2 U_{\tilde{b}11} \left( \frac{1}{3}s_W^2 + c_W^2 \right) \right] \nonumber \\ g_{\tilde{b},12}^{\phi _2^0}= & {} U_{\tilde{b}12}^* \left[ - \frac{U_{\tilde{b}21} m_b \mu }{c_\beta } + \frac{2}{3} U_{\tilde{b}22} s_\beta m_z^2s_W^2 \right] \nonumber \\&+ U_{\tilde{b}11}^* \left[ -\frac{U_{\tilde{b}22}m_b \mu ^*}{c_\beta } + s_\beta m_Z^2 U_{\tilde{b}21} \left( \frac{1}{3}s_W^2 + c_W^2 \right) \right] \nonumber \\ g_{\tilde{b},21}^{\phi _2^0}= & {} U_{\tilde{b}22}^* \left[ - \frac{U_{\tilde{b}11} m_b \mu }{c_\beta } + \frac{2}{3} U_{\tilde{b}12} s_\beta m_z^2s_W^2 \right] \nonumber \\&+\, U_{\tilde{b}21}^* \left[ -\frac{U_{\tilde{b}12}m_b \mu ^*}{c_\beta } + s_\beta m_Z^2 U_{\tilde{b}11} \left( \frac{1}{3}s_W^2 + c_W^2 \right) \right] \nonumber \\ g_{\tilde{b},22}^{\phi _2^0}= & {} U_{\tilde{b}22}^* \left[ - \frac{U_{\tilde{b}21} m_b \mu }{c_\beta } + \frac{2}{3} U_{\tilde{b}22} s_\beta m_z^2s_W^2 \right] \nonumber \\&+\, U_{\tilde{b}21}^* \left[ -\frac{U_{\tilde{b}22}m_b \mu ^*}{c_\beta } + s_\beta m_Z^2 U_{\tilde{b}21} \left( \frac{1}{3}s_W^2 + c_W^2 \right) \right] .\nonumber \\ \end{aligned}$$

(52)

Finally, for the $\mathcal {CP}$-odd states $\chi _1^0$ and $\chi _2^0$ the Higgs-sbottom couplings are:

$$\begin{aligned} g_{\tilde{b},11}^{\chi _1^0}= & {} i\frac{m_b}{c_\beta } \left[ -A_b^* U_{\tilde{b}12}^* U_{\tilde{b}11} + A_b U_{\tilde{b}11}^* U_{\tilde{b}12} \right] \nonumber \\ g_{\tilde{b},11}^{\chi _2^0}= & {} i\frac{m_b}{c_\beta } \left[ -\mu U_{\tilde{b}12}^* U_{\tilde{b}11} + \mu ^* U_{\tilde{b}11}^* U_{\tilde{b}12} \right] \nonumber \\ g_{\tilde{b},12}^{\chi _1^0}= & {} i\frac{m_b}{c_\beta } \left[ -A_b^* U_{\tilde{b}12}^* U_{\tilde{b}21} + A_b U_{\tilde{b}11}^* U_{\tilde{b}22} \right] \nonumber \\ g_{\tilde{b},12}^{\chi _2^0}= & {} i\frac{m_b}{c_\beta } \left[ -\mu U_{\tilde{b}12}^* U_{\tilde{b}21} + \mu ^* U_{\tilde{b}11}^* U_{\tilde{b}22} \right] \nonumber \\ g_{\tilde{b},21}^{\chi _1^0}= & {} i\frac{m_b}{c_\beta } \left[ -A_b^* U_{\tilde{b}22}^* U_{\tilde{b}11} + A_b U_{\tilde{b}21}^* U_{\tilde{b}12} \right] \nonumber \\ g_{\tilde{b},21}^{\chi _2^0}= & {} i\frac{m_b}{c_\beta } \left[ -\mu U_{\tilde{b}22}^* U_{\tilde{b}11} + \mu ^* U_{\tilde{b}21}^* U_{\tilde{b}12} \right] \nonumber \\ g_{\tilde{b},22}^{\chi _1^0}= & {} i\frac{m_b}{c_\beta } \left[ -A_b^* U_{\tilde{b}22}^* U_{\tilde{b}21} + A_b U_{\tilde{b}21}^* U_{\tilde{b}22} \right] \nonumber \\ g_{\tilde{b},22}^{\chi _2^0}= & {} i\frac{m_b}{c_\beta } \left[ -\mu U_{\tilde{b}22}^* U_{\tilde{b}21} + \mu ^* U_{\tilde{b}21}^* U_{\tilde{b}22} \right] . \end{aligned}$$

(53)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Funded by SCOAP³.

Reprints and permissions

About this article

Cite this article

Liebler, S., Patel, S. & Weiglein, G. Phenomenology of on-shell Higgs production in the MSSM with complex parameters. Eur. Phys. J. C 77, 305 (2017). https://doi.org/10.1140/epjc/s10052-017-4849-y

Download citation

Received: 27 March 2017
Accepted: 20 April 2017
Published: 12 May 2017
DOI: https://doi.org/10.1140/epjc/s10052-017-4849-y

Phenomenology of on-shell Higgs production in the MSSM with complex parameters

Abstract

Similar content being viewed by others

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

Heavy Higgs decays into sfermions in the complex MSSM: a full one-loop analysis

Neutral Higgs production at proton colliders in the CP-conserving NMSSM

1 Introduction

2 The MSSM with complex parameters

2.1 Gluino and squark sector

2.2 Higgs sector

2.3 Higgs mixing at higher orders

2.4 Wave function normalisation factors for external Higgs bosons

3 The gluon-fusion cross section

3.1 Lowest-order cross section

3.2 Higher-order contributions

4 The program SusHi and the extension SusHiMi

5 Numerical results

5.1 Squark contributions in the light-stop inspired scenario

5.2 Admixture of Higgs bosons in the \(m_h^{\mathrm {mod}+}\)-inspired scenario

5.3 \(\Delta _b\) corrections in the \(m_h^{\mathrm {mod}+}\)-inspired scenario

6 Remaining theoretical uncertainties

7 Conclusions

Notes

References

Acknowledgements

Author information

Authors and Affiliations

Corresponding author

Appendix A: Formulas: Higgs-quark and Higgs-squark couplings

Rights and permissions

About this article

Cite this article

Navigation

Phenomenology of on-shell Higgs production in the MSSM with complex parameters

Abstract

Similar content being viewed by others

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

Heavy Higgs decays into sfermions in the complex MSSM: a full one-loop analysis

Neutral Higgs production at proton colliders in the CP-conserving NMSSM

1 Introduction

2 The MSSM with complex parameters

2.1 Gluino and squark sector

2.2 Higgs sector

2.3 Higgs mixing at higher orders

2.4 Wave function normalisation factors for external Higgs bosons

3 The gluon-fusion cross section

3.1 Lowest-order cross section

3.2 Higher-order contributions

4 The program SusHi and the extension SusHiMi

5 Numerical results

5.1 Squark contributions in the light-stop inspired scenario

5.2 Admixture of Higgs bosons in the \(m_h^{\mathrm {mod}+}\)-inspired scenario

5.3 \(\Delta _b\) corrections in the \(m_h^{\mathrm {mod}+}\)-inspired scenario

6 Remaining theoretical uncertainties

7 Conclusions

Notes

References

Acknowledgements

Author information

Authors and Affiliations

Corresponding author

Appendix A: Formulas: Higgs-quark and Higgs-squark couplings

Appendix A: Formulas: Higgs-quark and Higgs-squark couplings

Rights and permissions

About this article

Cite this article

Share this article

Search

Navigation