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

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 $\texttt{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 $\texttt{SusHi}$ called $\texttt{SusHiMi}$.


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 CP is violated and leads to mixing between CP-even and -odd eigenstates. The reason that an analysis for the general case taking account the possibility of 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 CP-violating effects, for the two main Higgs production channels at the LHC, which can be used as input for future experimental analyses in 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 CP-even and CP-odd Higgs states form three admixed Higgs mass eigenstates h a , a ∈ {1, 2, 3}. Complex parameters in the MSSM give rise to additional sources of 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 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 processesgluon fusion and bottom-quark annihilation -in the SM and the MSSM with real parameters we refer to Section 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 CP-odd Higgs boson neglecting CP-even and -odd Higgs mixing were discussed in Ref. [18]. Refs. [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 Nextto-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 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 O(α t α s ) [30] and O(α 2 t ) [31,32] in the MSSM with complex parameters. 1 A detailed 1 Another approach to calculate the Higgs boson sector in the MSSM with complex parameters is based on the 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 CP-conserving, at higher orders an admixture of all three neutral Higgs bosons, i.e. the two CP-even Higgs bosons h, H and the CP-odd Higgs boson A, is induced. The case where the light Higgs boson describes the SM-like Higgs at ∼ 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 Section 2. We move to the description of the gluon-fusion cross section in Section 3, where we discuss the calculation of the cross section at LO and the applicability of higher-order corrections. Next we introduce in Section 4 the code SusHi and its extension SusHiMi, which we use for our phenomenological studies carried out in Section 5. We discuss the remaining theoretical uncertainties in Section 6. Lastly, we conclude in Section 7 and list Higgs-(s)quark couplings in Appendix A.

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 2 12 are rotated away. Those independent phases are the ones of the soft-breaking gaugino masses M 1 and M 3 , the Higgsino mass parameter µ and trilinear soft-breaking couplings A f , f ∈ {e, µ, τ, u, d, c, s, t, b}. 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.

Gluino and squark sector
The gluinog does not mix with other fields and enters the Lagrangian in the form renormalisation group improved effective potential approach and implemented in e.g. the code CPsuperH [33,34].
where mg is the absolute value of the complex soft-breaking parameter M 3 = mge iφ M 3 . 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 bottomquark 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, squarksq L,R of one generation mix into mass eigenstatesq 1,2 . The term of the Lagrangian containing the squark mass matrix of one generation is given by [37] Here X q := A q − µ * · {cot β, tan β}, where cot β and tan β apply to up-and down-type quarks, respectively. The soft-breaking masses M 2 q L and M 2 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 ≡ sin θ W . Contrarily, in the CP-violating MSSM the parameters A q = |A q |e iφ Aq and µ = |µ|e iφµ , 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 Uq having real diagonal elements and complex off-diagonal elements The squark masses (using the convention mq 1 ≤ mq 2 ) are calculated as the eigenvalues of Eq. (2). The fact that the left-handed soft-breaking parameter M 2 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 Refs. [47], which we have incorporated as a shift in the left-handed soft-breaking parameter M 2 q L in the sbottom sector, as obtained from FeynHiggs.
In the b/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 φ, where φ ∈ {h, H, A}, see Section 2.2 below, of the Higgs bosons to bottom quarks are given by the effective Lagrangian in terms of the left-handed and right-handed couplings g φ b L and g are the left-and right-handed projection operators, respectively. The explicit form of 2 The soft-breaking parameter M1 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. the couplings is given by with f h αβ = sin α/ cos β, f H αβ = cos α/ cos β and f A αβ = tan β (see also Refs. [37,54]). The effective Lagrangian provides a resummation of leading tan β-enhanced contributions entering via the quantity ∆ b . The leading QCD contribution to ∆ b has the form where α s is typically evaluated at an averaged SUSY scale µ d = (mb 1 +mb 2 +mg)/3, and the function I(a, b, c) is given by ). As one can see from Eq. (6), the leading contribution to ∆ b has an explicit dependence on the complex parameters M 3 and µ. In our numerical analysis below we use the value for ∆ 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 Section 4 below, both the ∆ 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 Eq. (4) and Eq. (5) in our leading-order expressions for the (loop-induced) gluonfusion process. For bottom-quark annihilation and the implementation of higher-order corrections to the gluon-fusion process, see Section 3, we use as a simplified version [37] 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 Section 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 Section 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.

Higgs sector
The MSSM contains two Higgs doublets with opposite hypercharges Y H 1,2 = ±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 CP-even (φ 0 1 , φ 0 2 ) and CP-odd (χ 0 1 , χ 0 2 ) components as follows 3 3 We note that the convention differs from the convention employed by FeynHiggs by a different sign of χ 0 1 and φ − 1 , which induces different signs in the corresponding elements of the matrices in Eq. (11) and Eq. (12) and the χ 0 1 couplings to (s)quarks displayed in the Appendix.
such that the Higgs potential V H in terms of the neutral Higgs states is given by The quadratic terms of V H contain the SUSY parameter |µ| 2 and the soft terms m H 1 , m 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 ξ between the Higgs doublets vanishes when the Higgs potential is minimised, making the Higgs sector of the MSSM CP-invariant at lowest order.
The tree-level neutral mass eigenstates {h, H, A, G} are related to the tree-level neutral fields {φ 0 1 , φ 0 2 , χ 0 1 , χ 0 2 } through a unitary matrix as follows Similarly, for the charged Higgs states one obtains where s x ≡ sin x, c x ≡ cos x. α, β n and β c are the mixing angles for the CP-even Higgs bosons (h, H), the neutral CP-odd states (A, G), and the charged states (H ± , G ± ), respectively. Minimising the Higgs potential leads to β := β n = β c at tree level. The masses of the charged Higgs bosons and the neutral CP-odd Higgs boson at tree level are given by 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 ± (m A ) and tan β := vu v d for the case of the MSSM with complex (real) parameters.

Higgs mixing at higher orders
The full mixing in higher orders takes place not just between {h, H, A}, but also with the Goldstone boson and the electroweak gauge bosons. In general, (6 × 6)-mixing contributions involving the fields {h, H, A, G, Z, γ} 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 × 3)-mixing matrix among {h, H, A}, since mixing effects with {G, Z, γ} 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 → {Z * , G * } → 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 → {Z * , G * } are completely absent at the one-loop level, even in case of CP violation in the squark sector. The one-loop contributions to gg → {Z * , G * } therefore have no dependence on the phases of complex parameters.
Thus, we focus our discussion on the contributions to the (3×3)-mass matrix M, which contains the tree level masses m 2 i on the diagonal and has non-zero (off-)diagonal self-energies involving the Higgs states. It enters the Lagrangian, with Φ = (h, H, A), as follows The propagator matrix is then given by (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 with effective self-energies that contain mixed termŝ where we suppressed the p 2 arguments of all terms and i = j = k.

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 {h 1 , h 2 , h 3 } vanish on-shell and the propagators of the external particles have unit residue, the correct onshell 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Ẑ factors contains the correction factors for the external Higgs bosons {h 1 , h 2 , h 3 } relative to the lowest-order mass eigenstates {h, H, A}. The matrix elementsẐ aj [39] (see also Refs. [35,37]) are composed of the root of the external wave function normalisation factor and the on-shell transition ratioẐ which are evaluated at the complex pole M 2 a . Here the indices {a, b, c} refer to the loop-corrected mass eigenstates, while {i, j, k} 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Ẑ corresponding to the (non-unitary) matrix As explained above, theseẐ factors provide the correct normalisation of a matrix element with an external on-shell Higgs boson h a , a ∈ {1, 2, 3}, at p 2 = M 2 a . The application of theẐ factors yields an expression of the amplitude A ha 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 {h, H, A} according to The ellipsis indicate additional mixing effects with Goldstone bosons and gauge bosons, which we neglect in our numerical analysis, see the discussion in Section 2.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 Section 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 ∈ {1, 2, 3}. Differences with respect to the calculation in the MSSM with real parameters are induced through 4 •Ẑ factors, which relate the amplitude for an external on-shell Higgs h a (in the mass eigenstate basis) to the amplitudes of both the CP-even lowest-order states h and H and the CP-odd state A, see Section 2.4.
• Non-vanishing couplings of squarks g Ã f ii to the pseudoscalar component A.
• Different left-and right-handed quark couplings g φ q L and g φ q R with φ ∈ {h, H, A}, see Section 2.1.

Lowest-order cross section
The LO production cross section of the mass eigenstates h a can be written as follows where τ ha = m 2 ha /s. The hadronic squared centre-of-mass energy is denoted by s, and the gluon-gluon luminosity by L gg . Therein, the partonic LO cross section for gg → h a is given by where G F denotes Fermi's constant, andẐ aφ are the elements of theẐ factor matrix. µ R is the renormalisation scale, which at LO only enters through the scale dependence of the strong coupling constant α s . We denote the cross section pp → h a , which involves one-loop diagrams 4 In the MSSM with real parameters only couplings involvingfi −fj − A with i = j are non-vanishing, and leftand right-handed quark Yukawa couplings are identical, in the production process pp → φ, as "LO cross section" despite the fact that it contains higherorder effects through the application of theẐ 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 L ⊃ G µν G µν φ with the gluon field strength G µν . The amplitudes of the second term stem from L ⊃G µν G µν φ, which involves the dual of the gluon field strength tensorG µν , 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 A ha,e and A ha,o , respectively. Similarly, we can split σ LO into σ e LO and σ o LO . For the two CP-even lowest-order mass eigenstates φ e ∈ {h, H} we obtain the amplitudes with where g φ e q L and g φ e q R are the couplings of the Higgs boson φ e to the left-and right-handed quarks, respectively. They are normalised to the SM Higgs-quark couplings. g φ ẽ qij are the couplings of the Higgs boson φ e to squarksq i andq j . The explicit expressions for the Higgs-squark and relative Higgs-quark couplings are listed in Appendix A. Similarly, for the CP-odd Higgs boson A we have with Within the previous formulas we use the notation and f (τ ) is given by Our result is consistent with Ref. [13], which however assumes g φ q L = g φ q R (see our discussion of this issue in Section 2.1) and does not take into account the mixing among the tree-level mass eigenstates φ ∈ {h, H, A}. All squark contributions, i.e.ã φ e q andã A q , enter the first term, A ha,e , in Eq. (25). Quark contributions to A ha,e which couple to the 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 A ha,o in Eq. (25) which couple to the CP-even lowest-order mass eigenstates φ e . All these terms are therefore denoted with the subscript A − . It should be noted that the amplitudes A φe,A − only arise due to the complex nature of the Yukawa couplings, which is a consequence of the incorporation of higherorder contributions entering via ∆ b , see Eq. (5), and our choice of working with a complex Yukawa coupling. Accordingly, the amplitudes A φe,A − are zero in the case of the MSSM with real parameters.

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 CP-even Higgs boson in an expansion around the threshold of Higgs production [79][80][81][82][83]. 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 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 ∆ b contributions. According to this approach, for the top-quark Yukawa coupling g φ t 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 Section 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 ∆ 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 where A ha,e , which has been given in Eq. (25), denotes the CP-even part of the LO amplitude including quark and squark contributions. Accordingly, this electroweak correction factor is only applied to the CP-even component of the LO and NLO cross section, see Section 4. The electroweak amplitude is given by [96] A ha,EW = − 3 8 · −Ẑ ah sin α cos β +Ẑ aH cos α sin β , with the abbreviation In Eq. (32) α EM denotes the electro-magnetic coupling, and is the sine of the weak mixing angle. M V and Γ V are the mass and the width of the heavy gauge bosons V ∈ {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] 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 CPconserving 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 CP-odd components proportional to g Ã qii are added as a complex component to the CP-even couplings. Since beyond LO we employ the simplified version of the ∆ b resummation according to Eq. (7), the higher-order quark contributions, both real and virtual, are of the same structure as in the 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 π of the various MSSM parameters using a cosine interpolation, see Refs. [111,112]. This interpolation makes use of on-shell stop-and sbottomquark masses defined at phases 0 and π. Thus, within the interpolated result we have to ensure the correct subtraction of the NLO contributions that have already been taken into account through ∆ b effects in the bottom-quark Yukawa coupling. This is done by expanding the ∆ b correction to next-to-leading order in the subtraction term. For a certain value of the phase φ z of a complex parameter z, the virtual NLO amplitude A φ NLO (φ z ) can be approximated using for each of the lowest-order mass eigenstates φ ∈ {h, H, A}. Here A φ NLO (0) is the analytical result for the MSSM with real parameters, and A φ NLO (π) is the analytical result with z → −z. Using the factors cos φ 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 µ is already apparent in the lowestorder diagrams of gg → φ, the phase of M 3 only enters through the NLO virtual corrections. Besides the ∆ 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 µ. 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 Ã qii 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 π, 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 Section 6 and assign a conservative estimate of the uncertainty in our numerical analysis). We focus here on the gluon fusion amplitude withoutẐ factors, since in theẐ 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 µ in the LO cross section, which is fully taken into account. The dependence of the NLO amplitude on the phases of A q and µ 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Ẑ 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 O(α t α s ) were approximated with an interpolation before the full phase dependence of the corresponding two-loop corrections at O(α 2 t ) 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 Ã qii arise away from phases 0 and π, 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 σ e NLO and σ o NLO individually, following Ref. [21], by defining the NLO correction factors C e and C o : The amplitudes are given by A ha,e µ F denotes the factorisation scale, and β 0 = 11/2 − n f /3 with n f = 5. Note that the LO amplitudes A ha,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 The terms ∆σ 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 Ã qii , which are added to the CP-even components ∆σ e . The real corrections can be split in ∆σ e and ∆σ o since no interference terms arise.

The program SusHi and the extension SusHiMi
SusHi input file in SLHA-style

SusHiMi initialization
Higgs masses andẐ factors by FeynHiggs Interpolated results for squark contributions "CP-even" "CP-odd" 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 CP violation in the Higgs sector. Following our discussion in Section 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 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 "CP-even" calculation for σ e NLO and a "CP-odd" calculation for σ o NLO according to Eq. (37). Thus, the total gluon-fusion cross section is the sum of the two parts We obtain the result beyond LO QCD through 8 whereas σ  7), named ∆ b1 . σ t,e N k LO, EFT and σ t,o N k LO, EFT 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 β = 1) with mass m ha , respectively. This K-factor is subsequently reweighted with the exact LO cross section. For this purpose the employed LO cross sections σ t,e LO and σ t,o LO are again evaluated as discussed in Section 3 with fullẐ factors, but include only the top-quark contribution. They are multiplied with the K-factors in σ t,e N k LO, EFT and σ t,o N k LO, EFT , respectively. Due to their small numerical impact in σ t,e/o N k LO, EFT 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 and theẐ 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 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 GeV is usually assumed to have a dominant CP-even component, which is also the case in the scenarios which are considered in our numerical discussion. For the 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 σ t,e LO and σ t,o LO . 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 2 ha /s → 1, which we can take into account up to 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 CP-even component of the light Higgs boson and thus match the precision of the SM prediction. The electroweak correction factor δ lf EW multiplied in the "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 β enhanced sbottom effects in the LO QCD contributions. We add the results as follows 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 ∆ b1 . We add and subtract the LO QCD cross section using the full resummation according to Eq. (5), named ∆ 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 Section 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 β = 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 |Ẑ ah g h b +Ẑ aH g H b | 2 + |Ẑ aA g A b | 2 , which includes tan β-enhanced squark effects through ∆ b according to Eq. (7). This procedure equals the application of a K-factor on the full LO cross section includingẐ 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 Though, due to the similarity of both approaches we only discuss Higgs production through bottom-quark annihilation with simplified ∆ b resummation.

Numerical results
For our numerical analysis we slightly modify two standard MSSM scenarios introduced in Ref. [116], namely the m mod+ h 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 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). 9 For the m mod+ h -inspired scenario we choose for vanishing phases of the complex parameters: We use for the SM  The depicted value of α s is only used for the evaluations of FeynHiggs, for the cross sections the value of α 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 µ 0 R and µ 0 F , respectively, is (µ 0 R , µ 0 F ) = (m ha /2, m ha /2) for gluon fusion and (µ 0 R , µ 0 F ) = (m ha , m ha /4) for bottom-quark annihilation. More details are described in Section 6. Whereas for the m mod+ h -inspired scenario we pick heavy Higgs bosons through m H ± = 900 GeV with tan β = 10 and 40 for the study of ∆ b effects, we choose m H ± = 500 GeV with tan β = 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 CP-even Higgs boson, where they reduce the gluon-fusion cross section by about ∼ 90%. The Higgs boson masses and theẐ 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 β = 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φ A t and M 3 = mge iφ 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 µ 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 µ and tan β ≥ 10, and so X t ≈ 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 ∈ {u, d, s, c} and A l for l ∈ {e, µ}. 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 µ [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Ẑ factors) and its effect on production cross sections. Lastly we discuss ∆ b corrections in the context of the m mod+ h -inspired scenario with large tan β, 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 Section 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 Section 6.

Squark contributions in the light-stop inspired scenario
We start with a discussion of squark effects to the Higgs boson cross section σ(gg → h i ) for all three Higgs bosons h i in the context of the light-stop inspired scenario with m H ± = 500 GeV and tan β = 16. The variation of the light Higgs boson mass m h 1 as well as the stop masses is depicted in Fig. 3 (a) as a function of φ At . 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 308 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 φ At and φ 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.  The black, dot-dashed curve depicts the cross section with top-quark and bottom-quark con-tributions and electroweak corrections in the production amplitudes only, i.e. in the formulas of Section 4 we omit all squark contributions which enter either directly or through ∆ b . Note however that squark contributions are always part of theẐ factors. Due to the decoupling with large values of m H ± 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. 3 (b)) 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 ∆ 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 Section 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). Fig. 3 (b) also includes the cross section for h 1 normalised to the cross section of a SM Higgs boson with the same mass.
Here the ∼ 20% reduction due to squark effects is apparent once again, whereas the quarkinduced 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 φ M 3 leads to a very similar picture, even though the light Higgs mass variation is not as pronounced and the stop masses are unaffected. Moreover, in the comparison of the simplified and the full resummation of ∆ 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. 4 (a) and (b) we show the gluon-fusion cross sections of the heavy Higgs bosons h 2 and h 3 , respectively, as a function of φ At . 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, σ/σ 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 φ At , and therefore the associated phase space effect is small. For vanishing phase φ At = 0 it is known that squark effects are huge and reduce the cross section by ∼ 89% (h 2 ) and ∼ 22% (h 3 ) [119]. These squark effects are strongly dependent on the phase φ At and induce a large positive correction at phase φ At = π in case of h 2 . For h 3 the effects are not as pronounced, but still sizeable. The K-factor for both processes gg → h 2 and gg → 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 φ At 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 φ At is induced by the admixture of the two Higgs bosons throughẐ factors. We will discuss this feature in detail for the m mod+ h -inspired scenario in Section 5.2.
The phase dependence on φ 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Ẑ factors and through ∆ b induces a dependence on the gluino phase already for the LO TeV) light-stop, t β = 16, m H ± = 500 GeV 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 Section 6, since we obtain the largest relative interpolation uncertainty in the cross section variation with phases for the interpolation of the gluino phase φ M 3 .

h -inspired scenario
In this subsection we discuss the m mod+ h -inspired scenario with tan β = 10 and m H ± = 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 ±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 CP-violating interactions is taken into account, there can be a large mixing between the CP-even and CP-odd neutral Higgs states. This feature is clearly visible in Fig. 6. The dependence on the phase φ At is seen to be closely correlated to 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. 6 (b) and (c) that in the overall contribution there will be sizeable cancellations between the phase dependencies of the separate contributions.

h -inspired scenario
We finally discuss the impact of ∆ b effects, which we investigate for the two heavy Higgs bosons in the m mod+ h -inspired scenario with tan β = 40. In this scenario the admixture between the two heavy Higgs bosons is again sizeable both as a function of φ At and as a function of φ 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 CP-even mass eigenstate h e and the predominantly CP-odd mass eigenstate h o , as defined in Section 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 CP-character (defining h o and h e ) are shown in Fig. 7 as a function of φ At and φ M 3 . One can see that the states h 2 and h 3 drastically change their CP character upon variation of the phases φ At and φ M 3 , while on the other hand the state h e is almost purely CP-even and h o is almost purely CP-odd for the whole range of phase values. It should be kept in mind in this context that |Ẑ 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 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 Fig. 6 (a), the solid and dashed curves refer to h e and h o , respectively. and φ M 3 . In both cases the behaviour for the full prediction, including the squark contributions, is dominated by ∆ 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 π, however, the ∆ b corrections can also give rise to a significant enhancement of the cross section. In particular, for φ M 3 the quantity ∆ b changes sign between φ M 3 = 0 and φ M 3 = π, such that the bottom-Yukawa coupling is suppressed for small values of φ M 3 and enhanced for φ M 3 values close to π as a consequence of the resummation of the ∆ 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 φ At and φ M 3 follows a similar trend as the LO cross section. In the plot on the right, the asymmetric K-factor dependence on φ M 3 is related to the direct dependence of ∆ b on the phase φ 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 ∆ b contributions to the bottom-quark Yukawa coupling is taken into account. For the latter both the results for the full (∆ b2 , blue) and the simplified (∆ b1 , red) resummation are shown. While the the pure LO squark contributions are seen to have a moderate effect, it can be seen that the incorporation of the resummation of the ∆ 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 ∆ b contribution very well. The curves corresponding to ∆ b2 and ∆ b1 hardly differ from each   other both for the variation of φ At and φ M 3 . As before all curves include the sameẐ 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 ∆ 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 CP-violating effects relevant for Higgs boson production through gluon fusion in the MSSM -squark loop contributions, admixtures througĥ Z factors and resummation of ∆ b contributions -for completeness we also briefly discuss the bottom-quark annihilation cross section for the m mod+ h -inspired scenario with tan β = 40. The corresponding cross section is shown in Fig. 10 as a function of φ At and φ M 3 . For such a large value of tan β this cross section exceeds the gluon-fusion cross section by far. It shows a very significant dependence on the phases φ At and φ M 3 , which is mainly induced by the ∆ b contribution.

Remaining theoretical uncertainties
In the previous section we analysed our cross section predictions regarding CP-violating effects entering via squark loop contributions,Ẑ factors and ∆ 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+α 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+α s uncertainties: The fitted parton distribution functions (PDF) and the associated value of α 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+α s uncertainties are mostly a function of the Higgs boson mass m ha . 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 β = 1) as a function of its mass. For Higgs masses in the range between 50 GeV and 1 TeV the typical size of PDF+α s uncertainties for gluon fusion is ±(3 − 5)% following the prescription of Ref. [132]. They increase up to ±10% for Higgs masses up to 2 TeV. For bottom-quark annihilation they are in the range ±(3 − 8)% for Higgs masses between 50 GeV and 1 TeV and up to ±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 Section 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 β 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 lightstop inspired scenario, which is at the lower edge of the scale uncertainty depicted in Fig. 3 (b). Other uncertainties discussed in Ref. [119] are renormalisation and factorisation scale uncertainties and an uncertainty related to higher-order contributions to ∆ 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 (µ 0 R , µ 0 F ) = (m ha /2, m ha /2) for gluon fusion and (µ 0 R , µ 0 F ) = (m ha , m ha /4) for bottom-quark annihilation. We obtain the scale uncertainty by taking the maximal deviation from the central scale choice ∆σ obtained from the additional scale choices ( We perform this procedure individually for all three cross sections in Eq. (41) and then obtain the overall absolute uncertainty through 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 σ ± ∆σ scale .
• In order to display the propagation of an uncertainty arising from higher-order contributions to ∆ b to our cross section calculation, we vary the value of ∆ b obtained from FeynHiggs by ±10%. This variation by ±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 ∆σ resum and assign an uncertainty band of σ ± ∆σ 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 σ(φ z ) following Eq. (38) not only for the correct phase φ z in Eq. (35), but also leave the phase within Eq. (35) constant, i.e. fixed to 0 and π. We call the obtained cross sections σ(0) and σ(π). For each value of φ z we take the difference ∆σ int = sin 2 (φ z )|σ(0) − σ(π)|/2. It is reweighted with sin 2 (φ z ), since we know that our result is correct at phases 0 and π. The obtained uncertainty band is given by σ ± ∆σ 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 β = 16 and the resummation uncertainty for the m mod+ h -inspired scenario with tan β = 40. Fig. 11 shows the renormalisation and factorisation scale uncertainties ∆σ scale as before and in addition the above described interpolation uncertainty ∆σ int , which in case of the variation of φ 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 gluonfusion 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 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 π, where the interpolated result in the MSSM with complex parameters merges the known result of the MSSM with real parameters. For the variation of φ At 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 φ M 3 .
The described ∆ b uncertainties are depicted in Fig. 12. Since ∆ b crosses 0 as a function of φ 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. 12(b). Even for the large value of tan β chosen here the assigned ∆ b uncertainty of ±10% is much smaller than the scale uncertainty of the displayed cross sections. Despite the different behaviour with the phases φ At and φ 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 β, where the cross section of bottom-quark annihilation exceeds the gluon-fusion cross section.

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 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 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 × 3)-mixing of the lowest-order mass eigenstates of the Higgs bosons {h, H, A} into the loop-corrected mass eigenstates {h 1 , h 2 , h 3 } 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 β-enhanced sbottom contributions in ∆ 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 ∆ 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 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 CP-even component of the Higgs bosons mediated by light quarks, have been added to the 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 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 mod+ h . 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 ∆ 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 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 ∆ 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 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 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.
In SusHiMi the Higgs-(s)quark couplings are expressed in terms of the CP-even and CP-odd neutral gauge eigenstates φ g ∈ {φ 0 1 , φ 0 2 } and χ g ∈ {χ 0 1 , χ 0 2 }, respectively. In order to obtain the couplings of the squarks with the lowest-order mass eigenstates φ ∈ {h, H, A, G} the gauge eigenstates are rotated using the tree-level mixing matrix R as depicted in the following Feynman diagrams with v = 2m W /g = 1/ √ 2G F and the tree-level mixing matrix R given in Eq. (11). At the amplitude level the results will then also be multiplied with the correspondingẐ factor.
The couplings between the gauge eigenstates and the third generation quarks are g q L = 1/ cos β for φ 0 1 and χ 0 1 and g q L = 1/ sin β for φ 0 2 and χ 0 2 . For ∆ b corrections we refer to Eq. (5) and Eq. (7).
The couplings between the gauge eigenstates and the third generation squarks contain terms from the squark mass diagonalisation matrix Uq (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 They are obtained with MaCoR [133,134]. For the CP-even state φ 0 1 we have the stop couplings (using s β ≡ sin β, c β ≡ cos β, t β ≡ tan β, s W ≡ sin θ W and c W ≡ cos θ W ):