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

We discuss neutral Higgs boson production through gluon fusion and bottom-quark annihilation in the CP-conserving $\mathbb{Z}_3$-invariant Next-to-Minimal Supersymmetric Standard Model (NMSSM) at proton colliders. For gluon fusion we adapt known asymptotic expansions in supersymmetric particles for the inclusion of next-to-leading order contributions of squarks and gluinos from the Minimal Supersymmetric Standard Model (MSSM) and include electro-weak corrections involving light quarks. Together with the resummation of higher order sbottom contributions in the bottom-quark Yukawa coupling for both production processes we thus present accurate cross section predictions implemented in a new release of the code SusHi. We elaborate on the new features of an additional SU$(2)_L$ singlet in the production of CP-even and -odd Higgs bosons with respect to the MSSM and include a short discussion of theory uncertainties.


Introduction
After the discovery of a scalar boson at the Large Hadron Collider (LHC) [1,2] in 2012 an essential task of particle physicists is to reveal the nature of the Higgs-like state and thus the nature of electro-weak symmetry breaking. Apart from deviations from the Standard Model (SM) prediction of the properties of the found Higgs-like state, further work includes the search for additional less and/or more massive scalar bosons, which can nicely be accommodated in supersymmetric models. The Next-to-Minimal Supersymmetric Standard Model (NMSSM) extends the Minimal Supersymmetric Standard Model (MSSM) by an SU(2) L singlet and allows the dynamical generation of the µ-term through electro-weak symmetry breaking [3]. The latter singlet-doublet mixing term in the superpotential lifts the MSSM tree-level upper bound of the Higgs mass given by the Z-boson mass. Thus, the NMSSM can easily accommodate the SM-like Higgs boson with a mass close to 125 GeV. Whereas for the calculation of the NMSSM Higgs spectrum and branching ratios various spectrum generators are available and include higher orders in perturbation theory (see Section 3), the calculation of neutral Higgs production cross sections did not exceed leading order (LO) in quantum chromodynamics involving squarks and gluinos (SQCD) [4] and did not include electro-weak corrections.
It is therefore timely to present the missing ingredients and a code for the calculation of accurate neutral Higgs production cross sections in the NMSSM, where the five neutral Higgs bosons are predominantly generated through gluon fusion and bottom-quark annihilation at a proton collider. For this purpose we extend the code SusHi [5]. For the time being we restrict our implementation to the real NMSSM without additional CP violation, such that CP-even H 1 , H 2 and H 3 and CP-odd Higgs bosons A 1 and A 2 can be distinguished in the Higgs sector. Most recent efforts related to Higgs physics at the LHC are summarized in the reports of the LHC Higgs cross section working group [6][7][8]. The SM Higgs is mainly produced through gluon fusion, where the Higgs-gluon coupling is mediated through virtual top-and bottom-quarks [9]. Higher order QCD corrections at next-to-leading order (NLO) are of large importance [10][11][12]. In the effective theory of a heavy top-quark the inclusive cross section is known to next-to-next-to-leading order (NNLO) in QCD [13][14][15], in addition finite top-quark mass effects at NNLO were calculated [16][17][18][19][20]. Beyond NNLO QCD effects are accessible through resummation [21][22][23][24][25][26] and electro-weak corrections are known [27][28][29]. Meanwhile next-to-NNLO (NNNLO) QCD contributions were estimated in the so-called threshold expansion [30][31][32][33], but they are not further considered in this publication.
The SM results for Higgs production through gluon fusion can be adjusted to the MSSM and the NMSSM through a proper reweighting of the Higgs couplings to quarks. However, the gluon fusion process can also be mediated through their superpartners, the squarks. With respect to the MSSM the only generically new ingredient, which goes beyond the projection of the physical Higgs bosons onto the neutral components of the two Higgs doublets, are couplings of the NMSSM singlet to squarks, since no couplings of the singlet to quarks or gauge bosons are present in the tree-level Lagrangian. It is therefore of importance to include squark contributions to gluon fusion at the highest order possible, even though they decrease in size with increasing squark masses. For the pseudoscalars A i squark contributions to gluon fusion are only induced at NLO, which motivates to go beyond just LO squark contributions for all Higgs bosons. For this purpose we adapt the works of Refs. [34][35][36] for the MSSM  , where m φ denotes the Higgs mass and M a generic SUSY mass. In contrast to the MSSM we are at present only working in this expansion of inverse SUSY masses and do not include an expansion in the so-called VHML, the vanishing Higgs mass limit (m φ → 0) for the SQCD contributions, as implemented in evalcsusy [37][38][39] or discussed in Ref. [40]. In the latter limit higher order stop-induced contributions up to NNLO level are known [41][42][43] and were partially included in previous discussions of precise MSSM neutral Higgs production cross sections [44]. Although for a pure CP-odd singlet component NNLO stop induced contributions are the first non-vanishing contributions to the gluon fusion cross section, we leave an inclusion of these to future work. For completeness, we add that in the MSSM a numerical evaluation of NLO squark/quark/gluino contributions was also reported in Ref. [45,46], whereas Refs. [47][48][49] presented analytic results for the pure squark induced NLO contributions. Electro-weak contributions to the gluon fusion production process mediated through light quarks [28,29] can be adjusted from the SM to the MSSM [5] and similarly to the NMSSM and are known to capture the dominant fraction of electro-weak contributions for a light SM-like Higgs with a mass below the top-quark mass, whereas they are generically small for larger Higgs masses.
For large values of tan β, the ratio of the vacuum expectation values of the neutral components of the two Higgs doublets, the bottom-quark Yukawa coupling is enhanced, such that the bottom sector gets more important for gluon fusion, and the associated production with a pair of bottom-quarks pp → bbφ is significantly enhanced. SusHi includes bottom-quark annihilation bb → φ, which in case of non-tagged final state b-quarks is a good theoretical approach, since it resums logarithms through the b-parton distribution functions. The latter process is known as five-flavour scheme (5FS) up to NNLO QCD [50,51] and can be easily reweighted from the SM to the MSSM/NMSSM by effective couplings [52,53]. In the NMSSM the singlet does not couple to quarks at LO, however taking into account the singlet induced component into the resummation of higher order sbottom effects is mandatory, since also the singlet to sbottom couplings are enhanced by tan β.
The new release of SusHi thus provides gluon fusion cross sections at NLO QCD taking into account the third generation quarks and their superpartners, the squarks, for all the five neutral Higgs bosons of the NMSSM. The squark and squark/quark/gluino contributions are implemented in asymptotic expansions of heavy SUSY masses. Electro-weak corrections induced by light quarks through the couplings of the Higgs bosons to Z and W ± bosons can be added consistently like in the MSSM. Similarly, the NNLO top-quark induced contributions are included. In addition, sbottom contributions can be resummed into an effective bottom-quark Yukawa coupling, also taking into account the additional singlet to sbottom couplings. The latter also applies to the calculated bottom-quark annihilation cross section at NNLO QCD. All features SusHi provides for the MSSM are available for the NMSSM as well, in particular distributions with respect to the (pseudo)rapidity and transverse momentum of the Higgs boson under consideration can be obtained. Left for future work is a link to MoRe-SusHi [54] to allow for the calculation of momentum resummed transverse momentum distributions.
We proceed as follows: We start with a discussion of the theory background in Section 2, where we elaborate on the NMSSM Higgs sector and the calculation of the gluon fusion cross section. Then we present the NLO virtual amplitude for gluon fusion as well as the calculation of bottom-quark annihilation including the resummation of sbottom induced contributions to the bottom-quark Yukawa coupling in the NMSSM. Subsequently we comment on the implementation in SusHi in Section 3, before we investigate the phenomenological features of the singlet-like Higgs boson in the CP-even and CP-odd sector with regard to Higgs production in Section 4. We also include a short discussion of theory uncertainties. Finally, we conclude, and present the Higgs-squark-squark couplings in Appendix A.

Theory background
In this section we discuss the Higgs sector of the CP-conserving Z 3 -invariant NMSSM, before we proceed to the resummation of tan β enhanced sbottom contributions in the bottomquark Yukawa coupling. Subsequently we move to the discussion of the Higgs production cross section in gluon fusion, where we present the adapted formulas for the NLO SQCD virtual amplitude, and finally comment on the consequences of the additional singlet to bottom-quark annihilation.

2.1
The Higgs sector of the CP-conserving NMSSM Our notation of the Higgs sector of the CP-conserving Z 3 -invariant NMSSM closely follows Ref. [55]. For an extensive NMSSM review we refer to Ref. [3]. The superpotential can be written in the form where W MSSM equals the superpotential of the MSSM without µ-term.Ŝ denotes the additional SU(2) L singlet superfield compared to the MSSM with the two SU(2) L doublet su-perfieldsĤ d andĤ u . ǫ ab contracts the SU(2) L doublet components. Since the singletŜ is a neutral field, it induces one additional CP-even and one additional CP-odd neutral Higgs boson as well as one additional neutralino compared to the MSSM. The soft-breaking terms include the scalar components H d , H u and S of the superfields and are given by The soft-breaking mass m s can be derived from the minimization conditions of the tadpole equations (in addition to m 2 H d and m 2 Hu like in the MSSM), whereas A λ and A κ are usually considered input parameters. A λ can be alternatively replaced by the charged Higgs mass m H ± as input parameter. The neutral components of the Higgs fields are decomposed according to where v d , v u and v s denote the vacuum expectation values (vev) and the fields with indices R and I are the CP-even and CP-odd fluctuations around them. An effective µ term is generated through the vev of the singlet which will be further used within this article. We do not present the explicit form of the mass matrices here, but refer to Ref. [55]. We define the CP-even gauge eigenstate basis we perform a prerotation in the CP-odd sector to obtain the MSSM pseudoscalar A and the Goldstone G in the form H ′I = (G, A, S I ). The prerotation R G is given by the ratio of vacuum expectation values tan β = v u /v d , such that the mass eigenstates A i with i ∈ {1, 2} are obtained by where R P is a (3 × 3)-matrix, which however only consists of a (2 × 2)-mixing block, whereas R P i1 = R P 1i = 0 for i = 1 and R P 11 = 1. For Higgs production the Goldstone boson does not need to be considered. In the following we make use of the notation "singlet-like Higgs boson", which refers to the CP-even/odd Higgs boson with the dominant fraction of the singlet component S in gauge eigenstates. For this purpose we define the singlet character |R S i,3 | 2 for H i and |R P i+1,3 | 2 for A i . In our discussion of cross sections we denote the Higgs boson by the letter φ, which can be replaced by any of the physical Higgs bosons H i or A i . For picking viable scenarios for phenomenological studies we refer to Ref. [56] for a recipe to obtain positive eigenvalues for the singlet-like CP-even and -odd Higgs by varying A κ between a minimal and a maximal value.
Whereas the singlet component S does not couple to quarks, F -terms induce a coupling of the singlet-like Higgs to squarks, which is of relevance for Higgs production. We present the Higgs-squark-squark couplings to the third generation of squarks in Appendix A. We point out that the singlet component mixes with the Higgs doublets proportional to λ and also the couplings of the singlet component to squarks are proportional to λ. It is thus possible to mostly decouple the singlet component by lowering the value of the parameter λ. The couplings to quarks can be easily translated from the MSSM by the correct projection on the neutral doublet components H R d , H R u and the pseudoscalar A and yield relative to the SM with i ∈ {1, 2, 3} in the CP-even and i ∈ {1, 2} in the CP-odd Higgs sector. The relative strength g φ f enters the Yukawa couplings in the form

Resummation of higher order sbottom contributions
It is well-known in the MSSM that tan β enhanced sbottom corrections to the bottom-quark Yukawa coupling can be treated in an effective Lagrangian approach [57][58][59] to be resummed.
For the case of the NMSSM just taking into account SQCD corrections the effective Lagrangian can be written in the form [60] Ref. [60] additionally presents the inclusion of SUSY electro-weak corrections proportional to the soft-breaking parameter A t . The inclusion of electro-weak corrections into ∆ b does not harm our subsequent discussion of SQCD corrections and can thus always be included in the bottom-quark Yukawa coupling entering gluon fusion and bottom-quark annihilation. Apart from a coupling of the bottom-quarks to the gauge eigenstate H 0 u also an effective coupling to the singlet S can be induced at loop level, the latter being proportional to λv u instead of µ. The sbottom corrections can be absorbed into effective Yukawa couplings, which read [60] for the three CP-even Higgs field H i and for the two CP-odd Higgs fields A i .

Gluon fusion cross section
After our discussion of the CP-conserving NMSSM Higgs sector and the resummation of sbottom contributions in the bottom-quark Yukawa coupling, we present the gluon fusion production cross section for a Higgs boson φ, which can be written in the form [5] σ with τ φ = m 2 φ /s and the hadronic centre-of-mass energy s. The factor σ φ 0 includes the LO partonic cross section. C φ encodes NLO terms of singular nature in the limitŝ → m 2 φ with the partonic centre-of-mass energyŝ. The gluon-gluon luminosity is given by the integral The contributions ∆σ φ gg , ∆σ φ gq , and ∆σ φ qq are the regular terms in the limitŝ → m 2 φ in the partonic cross section and arise from gg, gq and qq scattering, respectively. The LO contribution σ φ 0 is obtained by the formulas presented in Ref. [5] using the NMSSM couplings of the Higgs bosons to quarks and squarks. Similarly the contributions ∆σ φ xy are obtained from the MSSM by a proper replacement of the involved Higgs boson to quark and squark couplings. The factor C φ can be decomposed in the form with β 0 = 11/2 − n f /3 and n f = 5 as well as the factorization and renormalization scales, ∞ is the LO (one-loop) virtual amplitude in the limit of large stop and sbottom masses. Φ (2l) is the NLO (two-loop) virtual amplitude and equals the form factors H (2l) i for the CP-even Higgs bosons and A (2l) i for the CP-odd Higgs bosons, which are presented in the following Section 2.4 to account for NLO virtual contributions from quarks and squarks in an appropriate way. In the MSSM limit they correspond to the form factors of Refs. [34][35][36] except from a constant factor of −3/4. Contrary to the case of the MSSM we do not employ evalcsusy [38,39] to obtain the NLO amplitude in the limit of heavy top-quark and stop masses for the light Higgs, but use the expanded form factors presented in the following sections instead. Accordingly our implementation does not (yet) include approximate NNLO stop contributions as presented in Ref. [44] for the MSSM. The NNLO top-quark contributions in the heavy top-quark effective theory making use of Refs. [13,61] are included according to Eq. (29) of Ref. [5].
Lastly we comment on the inclusion of the electro-weak corrections to the gluon fusion production cross section. Similarly to the MSSM the full SM NLO electro-weak (EW) corrections [27] can be added to the top-quark induced result only, assuming complete factorization of EW and QCD effects [62]. We recommend the latter procedure only for a SM-like Higgs boson. Contrary the inclusion of electro-weak corrections due to light quarks [28,29], where the Higgs boson couples to either the Z or W ± boson, can be adjusted to the MSSM and accordingly the NMSSM in an appropriate way [63]. For this purpose the generalized couplings of the i-th CP-even NMSSM Higgs boson to the heavy gauge boson V ∈ {W ± , Z} need to be inserted in the formulas of Ref. [5]. The missing projection R S i3 on the singlet component reflects the fact that the singlet does not couple to gauge bosons. The CP-odd Higgs bosons do not couple to gauge bosons either, such that electro-weak corrections due to light quarks are absent.

NLO virtual amplitude
Regarding the implementation of two loop contributions, we closely follow Refs. [34][35][36] for the MSSM, which can be translated to the NMSSM. Their calculation at NLO is based on an asymptotic expansion in the masses of the supersymmetric particles. We can project the form factors onto the ones in gauge eigenstates according to The individual contributions in gauge eigenstates are presented in Section 2.4.1 for the CPeven and in Section 2.4.2 for the CP-odd Higgs bosons. We included the constant factor between Refs. [34][35][36] and our work in the above equations.

CP-even Higgs bosons
In this subsection we present the form factors for the CP-even Higgs bosons in gauge eigen- which includes the effective µ parameter defined in Eq. (4). All functions in H and S R,(2l) can be directly taken over from Refs. [34,36] having in mind the different convention in the sign of the µ parameter. For on-shell (OS) parameters (see Refs. [34,36]) and thus for our implementation the contribution F 2l t is shifted according to Section 3.3 of Ref. [36] and F 2l b according to Ref. [34]. The shift also applies to F 2l t and F 2l b entering the singlet contribution S R , since the differences in the prefactors being µ, λv d or λv u are not renormalized when taking into account SQCD contributions and therefore do not contribute to the described OS shifts.
We are left with a discussion about the inclusion of resummed sbottom contributions into the bottom-quark Yukawa coupling within the virtual corrections to gluon fusion, where care has to be taken to avoid a double-counting of NLO SQCD contributions. The naive is incorporated in the same way as in case of the MSSM [34,44]. The resummation as presented in Section 2.2 instead needs the subtraction of the tan β enhanced contributions to G 2l b multiplied with the corresponding coupling correction, in detail for the three CP-even Higgs bosons H i with the factor K ′ i being All occurrences of the bottom-quark Yukawa coupling in the two-loop amplitude are multiplied with the factor K i =g H i b /g H i b using Eq. (10), such that the shift reported in Eq. (19) avoids double-counting of the purely SQCD induced contributions at the two-loop level. Employing the expansion in heavy SUSY masses the NLO virtual contributions to neutral CP-even Higgs production in the NMSSM are now fully presented.

CP-odd Higgs bosons
We now turn to the case of the two CP-odd Higgs bosons, where we present the form factor in the basis H ′I after a prerotation from gauge eigenstates 1 . At LO only diagrams involving quarks coupling to the pseudoscalar A exist, such that the form factor at LO only consists of the part H I,(1l) A , whereas S I,(1l) equals zero. At NLO however couplings of A and S I to squarksq iqj for i = j induce contributions to Higgs production. The two-loop form factor presented in Ref. [35] for the MSSM can therefore be translated to Whereas the individual contributions to H I,(2l) A can be taken from Ref. [35], we present the contributions to S I,(2l) separately: The functions K 1l and R 2 can be found in Ref. [35]. The functions R ′ t 1 and R ′ t 4 are given by: In the bottom sector the relevant function yields: The shifts of individual contributions in case of OS parameters can be taken over from the MSSM case. The inclusion of resummed sbottom contributions to the bottom-quark Yukawa coupling needs the following shift in the two-loop form factor for the CP-odd Higgs bosons using Again we point out that our sign convention with respect to µ is opposite to Ref. [35] and all occurrences of the bottom-quark Yukawa coupling in the two loop amplitude are multiplied with the factor K i =g A i b /g A i b using Eq. (11).

Bottom-quark annihilation cross section in the 5FS
The generalization of the calculation of bottom-quark annihilation cross sections in the fiveflavour scheme (5FS) from the MSSM to the NMSSM case is straightforward by using the appropriate couplings of Higgs bosons to bottom-quarks. For this purpose the resummation of sbottom contributions as described in Section 2.2 is taken into account. For the specific case of the singlet-like Higgs boson we point out that in case the coupling to the bottom-quark vanishes (due to cancellations in the mixing with the Higgs doublets) a priori the coupling to sbottom squarks can still be present. This is not taken into account by the resummation procedure.

Implementation in SusHi
In the current implementation of neutral Higgs production in the real NMSSM within the code SusHi the Higgs mixing matrices as well as the Higgs masses have to be provided as input in SUSY Les Houches Accord (SLHA) form [64,65] and can be obtained by spectrum generators for the NMSSM.  [74,77]. Special attention needs to be paid to the renormalization of the stop and sbottom sector, which in the ideal form should be identical in the calculation of Higgs masses and mixing and the calculation of Higgs production cross sections. For the time being, SusHi either relies on the internal calculation of on-shell stop and sbottom sectors as described in the manual [5] or on the specification of the on-shell masses mq 1 and mq 2 and mixing angles θq in the input file. For both cases input files can be found in the folder /example within the SusHi tarball. The user is asked to check the meaning of output parameters of spectrum generators, i.e. the chosen renormalization scheme. If the user specifies the on-shell squark masses and mixing angles together with the on-shell soft-breaking parameters A t and A b by hand, she/he should make sure that in the stop sector A t as well as the on-shell top-quark mass m t , the on-shell stop masses mt 1 and mt 2 and the mixing angle θt fit the formula In the sbottom sector on-shell and tree-level masses on the other hand differ by a shift in the (1, 1)-element, see ∆M 2 L in Ref. [5]. Moreover we employ the scheme which works with a dependent bottom-quark mass m b , whereas A b is defined to be on-shell, see e.g. Refs. [78][79][80]. To allow for maximal flexibility the specification of on-shell squark masses and mixing angles is now also possible in case of the MSSM. The Block RENORMSBOT is not of relevance in such input files, since m b is chosen as dependent parameter, whereas the squark mixing angle θ b and the soft-breaking parameter A b are understood as renormalized on-shell.
Two options for the pseudoscalar Higgs mixing matrix are accepted as input by SusHi, namely the full Higgs mixing matrix, which corresponds to the multiplication R P R G in the above notation, but instead also the rotation matrix R P can be used as input. Following SLHA2 [65] the full matrix (R P R G ) ij is provided in Block NMAMIX and asks for entries ij with i ∈ {2, 3} and j ∈ {1, 2, 3}. The matrix R P ij can be specified in Block NMAMIXR, which only asks for entries ij with {i, j} ∈ {2, 3}. We point out that in contrast to other codes the Goldstone boson remains the first mass eigenstate, such that Block NMAMIXR does not ask for entries with i = 1 or j = 1. The elements of the CP-even Higgs boson mixing matrix are specified in Block NMHMIX [65]. The Higgs masses need to be given in Block MASS using entries 25, 35 and 45 for the CP-even Higgs bosons and 36, 46 for the CP-odd Higgs bosons.
The block Block EXTPAR still contains the gluino mass as well as the soft-breaking parameters for the third generation squark sector. Entry 23 for the µ parameter is however replaced by entry 65, where the effective value of µ needs to be specified. Moreover entry 61 asks for the choice of λ. SusHi extracts the vev v s from µ and λ. Since the Higgs sectors including their mixing are provided, there is no need to provide the parameters κ, A κ , A λ (or m H ± ) in the SusHi input, since they do not enter the couplings relevant for Higgs production. planned, planed, planned The Block SUSHI entry 2 specifies the Higgs boson, for which cross sections are requested. The CP-even Higgs bosons are numbered 11, 12 and 13, the CP-odd Higgs bosons 21 and 22. Similarly the options 11, 12 and 21 also work in the 2HDM and the MSSM and 11 and 21 in the SM. We note that SusHi is still compatible with input files with 0 (light Higgs), 1 (pseudoscalar) and 2 (heavy Higgs) as options for entry 2. Output files however stick to the new convention.
For the time being we emphasize that SusHi is not strictly suitable for very low values of Higgs masses m φ < 20 GeV, where quark threshold effects start to become relevant and also electro-weak corrections are not implemented. This statement mostly applies to studies of a very light CP-odd Higgs boson, which is poorly constrained by LEP experiments in contrast to a light CP-even Higgs boson [81].

Phenomenological study
In this section we elaborate on the phenomenological consequences of the additional SU(2) L singlet in the NMSSM with respect to the MSSM for neutral Higgs production. Neglecting the squark induced contributions to gluon fusion, the only consequence of the additional singlet component is another admixture of the three CP-even/two CP-odd Higgs bosons. However, no generically new contributions to Higgs boson production arise. This differs when taking into account squark induced contributions to gluon fusion due to the additional singlet to squark couplings. In particular for the CP-odd Higgs bosons squark contributions are only induced at the two-loop level due to the non-diagonal structure of the CP-odd Higgs bosons to squark couplings. Subsequently we work with two scenarios, start with their definition, present the Higgs boson masses and admixtures and then discuss the behaviour of cross sections, including the squark and electro-weak corrections to the gluon fusion cross section. Our studies are performed for a proton-proton collider with centre-of-mass (cms) energy of √ s = 13 TeV, as planned for the second run of the LHC. Lastly we add a short discussion of renormalization and factorization scale uncertainties as well as PDF+α s uncertainties for one of the two scenarios. GeV. We vary κ between 0.15 and 0.80 and thus vary the mass of the singlet-like Higgs component in particular in the CP-even Higgs sector. We work out the characteristics for the singlet-like component in the following discussion. The relevant input for SusHi is obtained with NMSSMCalc 1.03, which incorporates the leading two-loop corrections O(α s α t ) to the Higgs boson masses calculated in the gaugeless limit with vanishing external momentum. We request NMSSMCalc to work with an on-shell renormalized stop sector and add local modifications to the NMSSMCalc input routines to read in on-shell parameters rather than DR renormalized parameters 2 . These modifications guarantee identical on-shell stop masses in NMSSMCalc and SusHi. The renormalization of the sbottom sector on the other hand is performed SusHi-internally.
We also choose a second scenario S 2 , in which we vary λ to decouple the singlet-like Higgs from the Higgs doublets. The detailed choice of parameters is M 1 = 150 GeV, M 2 = 300 GeV, The lower bound at λ = 0.04 is to avoid tiny cross sections for a heavy singlet-like Higgs boson and to keep its mass below the SUSY masses thresholds to justify the NLO SQCD expansion employed for the gluon fusion cross section calculation.
Both our scenarios come along with rather light third generation squark masses at the low TeV scale. Contrary to the Higgs mass calculations the squark contributions completely decouple from Higgs production for heavy SUSY spectra. Our scenarios are chosen to flash the phenomenology of an additional singlet-like Higgs boson and thus do not always include a SMlike Higgs boson with mass ∼ 125 GeV and are partially under tension from LEP searches [81] (for low CP-even Higgs masses below 110 GeV) or LHC searches [82][83][84][85][86][87][88].
We add for both scenarios the relevant SM input, which includes the MS renormalized bottom-quark mass m b (m b ) = 4.20 GeV, which is translated into a bottom-quark pole mass of m b = 4.92 GeV. In SusHi we choose the renormalization scheme, where the bottom-quark pole mass enters all occurrences of heavy bottom-quark masses in the loops and the bottomquark Yukawa coupling for the gluon fusion cross section. Bottom-quark annihilation is based on the running MS renormalized bottom-quark Yukawa coupling. As pointed out in Ref. [44] the gluon densities are hardly dependent on the bottom-quark pole mass fit value of the PDF fitting groups, emphasizing that there is no need to adjust the bottom-quark pole mass to the PDF fit value for the calculation of the gluon fusion cross section. The top-quark pole mass equals m t = 173.3 GeV. The strong coupling constant α s (m Z ) is set to 0.1172 for the calculation of running masses, and is obtained from the corresponding PDF set for the cross section calculation. We choose MSTW2008 [89] at the appropriate order in perturbation theory. Our central scale choices for gluon fusion are m φ /2 for both renormalization and factorization scale, µ 0 R and µ 0 F respectively, and µ 0 R = m φ and µ 0 F = m φ /4 for bottom-quark annihilation.  Subsequently we start with a discussion of the singlet admixture and the masses of the three CP-even and the two CP-odd Higgs bosons, which we obtain through a link to NMSSMCalc 1.03 as explained beforehand. For scenario S 1 the singlet component as a function of κ for the three CP-even Higgs bosons is shown in Fig. 1 (a). Clearly, for low values of κ the lightest Higgs H 1 is mainly singlet-like, whereas with increasing κ the dominant singlet fraction moves from H 1 to H 2 and for large values of κ to H 3 . The sum of all singlet components yields i |R S i3 | 2 = 1. The masses of the CP-even Higgs bosons can be found in Fig. 1 (b), where with increasing κ the in mass increasing singlet component can be clearly identified. Close to κ ∼ 0.35 H 2 shows the most dominant singlet fraction, which will be later visible in the gluon fusion cross section. Scenario S 1 includes for κ > 0.3 a SM-like Higgs boson H 1 with a mass of m H 1 ∼ 125 GeV. For very small values of κ the decay H 2 → H 1 H 1 opens and leaves a characteristic signature for the SM-like Higgs boson H 2 . Note that a light singlet-like Higgs boson lifts the mass of the SM-like CP-even Higgs through singlet-doublet mixing, which for our example equals m H 2 ∼ 153 GeV for κ = 0.1. The region of small κ and a light CP-even singlet-like Higgs boson H 1 is largely constrained by the LEP experiments [81]. Fig. 1 (c) and (d) show the behaviour of the singlet admixture and the masses for the two CP-odd Higgs bosons in scenario S 1 as a function of κ. We point again to the region in the vicinity of κ ∼ 0.35, where the light CP-odd Higgs boson A 1 is a pure singlet-like CP-odd Higgs boson contrary to the CP-even Higgs boson H 2 , for which H R d and H R u components remain. The coupling of A 1 to quarks vanishes, but the coupling to squarks is still present due to the relatively large value of λ = 0.62, which will be apparent when calculating the gluon fusion cross section.  For scenario S 2 Fig. 2 shows correspondingly the singlet character and masses for the CPeven and CP-odd Higgs bosons. Due to the fixed value of µ = 1 √ 2 λv s the singlet-like Higgs boson increases in mass (proportional to κv s ) with decreasing λ and thus for small λ H 3 as well as A 2 clearly decouple from the other Higgs bosons. We will later use this setup to show the decoupling behaviour of the cross sections. Below λ < 0.05 both H 3 and A 2 have a singlet character, which exceeds |R S/P 33 | 2 > 0.999.

Scenario S 1 : Inclusive cross sections for √ s = 13 TeV
In this subsection we investigate the gluon fusion σ gg and bottom-quark annihilation σ bb cross sections for scenario S 1 for √ s = 13 TeV for a proton-proton collider. The subsequent statements are however hardly dependent on the cms energy and thus hold for the 7/8 TeV LHC runs as well as for more energetic runs. Fig. 3 shows the cross sections for the three CP-even Higgs bosons. Naturally the cross sections are strongly dependent on the Higgs mass, which are in turn a function of κ. Thus, the cross section for the second CP-even Higgs bosons H 2 tends to decrease with increasing κ. Crucial is the singlet admixture of the Higgs boson under consideration. The larger the singlet component |R S i3 | 2 , the smaller the coupling to quarks becomes and thus the more sensitive is the cross section to squark and electro-weak contributions. For H 2 we observe a cancellation of quark contributions through the admixtures with the SU(2) L doublets around κ ∼ 0.35, where in turn due to the generally small cross section squark but also electro-weak corrections to the gluon fusion cross section are of large relevance, see Fig. 3 (c) and (d). For small values of κ the decay H 2 → H 1 H 1 opens in addition to the large gluon fusion cross section for the singlet-like CP-even Higgs boson H 1 . The region is therefore constrained by LEP experiments [81]. Much more smooth is the behaviour for the bottom-quark annihilation cross section, where the direct coupling to bottom-quarks is related to the non-singlet character of the Higgs under consideration. In an interval around κ ∼ 0.35 bottom-quark annihilation even exceeds the gluon fusion cross section for H 2 despite the small value of tan β = 2.
We show the effect of squark and electro-weak contributions to gluon fusion for the three CP-even Higgs bosons in Fig. 3 (c) and (d). σ q+q gg in Fig. 3 (c) includes stop-and sbottomquark induced contributions at NLO SQCD on top of the quark induced contributions without electro-weak contributions and compares to the pure quark induced cross section σ q gg without electro-weak contributions. All cross sections include NLO QCD quark contributions and the NNLO QCD top-quark induced contributions in the heavy top-quark effective theory. Fig. 3 (d) accordingly shows the effect of electro-weak contributions induced by light quarks following Eq. (15) in combination with Ref. [44] in comparison to the quark and squark induced cross section σ QCD gg = σ q+q gg . Note that in all our figures σ gg corresponds to σ QCD+EW gg . As expected for H 2 the region with small quark contributions induced by the admixture with the H R d and H R u components is in particular sensitive to squark corrections. For the other Higgs bosons the squarks corrections in this scenario are incidentally all of the order of O(−10%) and mostly independent of κ. We note that the squark corrections are mainly induced by stop contributions, whereas sbottom induced contributions only account for a small fraction. Interestingly the squark contributions show an interference-like structure with a maximum and minimum around κ ∼ 0.35, whereas the relative electro-weak corrections are always positive. This can be understood from a sign change in the real part of the quark induced LO and NLO amplitude for H 2 at κ ∼ 0.35, which is of relevance for the squark contributions, whereas the imaginary part, more relevant for the electro-weak contributions, does not change its sign. The size of the electro-weak corrections for H 2 follows from a suppression of the couplings of the second lightest Higgs H 2 to quarks in contrast to the couplings to gauge bosons. To obtain a pure singlet-like Higgs boson in the CP-even Higgs sector, which neither couples to quarks and gauge bosons, rarely happens due to the mixing between both S R and H R d as well as S R and H R u for large values of λ. For the SM-like Higgs boson with a mass below the top-quark mass the electro-weak corrections by light quarks are typically of the order of O(+5%) and cover most of the SM-electroweak correction factor. On the other hand, for Higgs masses above the thresholds m φ ≫ 2m W or 2m Z the electro-weak corrections by light quarks are small. The structure visible for H 2 in Fig. 3 (c) for κ < 0.3 is induced by the thresholds 2m W and 2m Z , which the Higgs mass m H 2 crosses between κ = 0.1 and 0.3. We leave the distortion of distributions, in particular transverse momentum distributions, for such a scenario to future studies.
Similarly we depict the gluon fusion and bottom-quark annihilation cross sections for the two CP-odd Higgs bosons in Fig. 4 (a) and (c). The pure singlet-like Higgs boson A 1 at κ ∼ 0.352 is clearly apparent, since both cross sections vanish. The corrections through squark contributions as shown in Fig. 4 (b) are very large around κ ∼ 0.352, since squark contributions are not suppressed through Higgs mixing, although they only appear at NLO SQCD. Electro-weak corrections induced through light quarks are absent for the CP-odd Higgs bosons. We note that in the range κ = 0.351 − 0.353, where the LO QCD gluon fusion cross section for the light CP-odd Higgs boson A 1 are tiny < 10 −5 pb, the prediction for σ q+q gg with squark induced NLO SQCD contributions for A 1 is unreliable, since SusHi calculates NLO QCD contributions through the multiplication of one-loop and two-loop contributions, where the latter tend to be significantly larger than the former and can thus even induce negative cross sections. However, in these regions the tiny cross sections are not of relevance for current searches. The cross section σ q+q gg and the relative correction to the vanishing only quark induced cross sections of more than 100% need to be taken with care for the CP-odd Higgs bosons. The fact that for a pure singlet-like CP-odd Higgs boson the gluon fusion cross section at NLO SQCD completely vanishes due to the absence of a LO contribution motivates to take into account NNLO SQCD stop contributions as it was done for the light CP-even Higgs boson in Ref. [44]. A first estimate yields tiny, positive cross sections, but we leave an inclusion in SusHi to future work. The CP-even Higgs bosons in contrast have a LO squark induced contribution, which leaves σ q+q gg mostly well-behaved. Only in rare cases, where LO squark and quark contributions cancel, similar difficulties can arise.

Scenario S 2 : Inclusive cross sections for √ s = 13 TeV
In this subsection we examine the inclusive cross sections for scenario S 2 , which includes low values of λ and thus reflects the decoupling of the singlet-like Higgs. The gluon fusion cross section together with squark contributions and electro-weak corrections through light quarks are shown in Fig. 5 for the three CP-even Higgs bosons. With decreasing λ the gluon fusion cross section for H 2/3 rapidly decreases, which has two reasons: First the increasing mass naturally decreases the cross section, but second the singlet-like Higgs boson also decouples from the other two Higgs bosons -being apparent in Fig. 2. Thus, the indirect coupling to quarks is suppressed. Moreover also the direct couplings to squarks are proportional to λ and thus decrease in size with decreasing λ. The relative correction induced by squark contributions remains rather constant, see Fig. 5 (b). The small interference structure visible around λ = 0.22 − 0.23 stems from the interchange of the dominant singlet character between H 1 and H 2 . Remark that for a SM Higgs the decrease of the gluon fusion cross section due to the increase in mass is between m H = 500 GeV and m H = 1200 GeV only a factor of ∼ 47, whereas we observe a decrease of more than five orders of magnitude, thus mainly driven by the decoupling. Electro-weak corrections by light quarks as depicted in Fig. 5 (c) are completely absent for a heavy singlet-like Higgs boson H 3 , but show a similar pattern at large λ for H 2 as for small κ in scenario S 1 . The reason is that the Higgs mass m H 2 in both scenarios crosses the 2m Z and 2m W thresholds for the electro-weak corrections by light quarks. The bottom-quark annihilation cross section shows a similar decoupling behaviour and is thus not explicitly shown.
For the CP-odd Higgs bosons we show the corresponding decoupling limit in Fig. 6, where the gluon fusion cross section and the bottom-quark annihilation cross section presented in Fig. 6 (a) and (c) respectively decrease dramatically for the singlet-like Higgs bosons A 2 . Again the squark corrections shown in Fig. 6 (b) remain rather constant, emphasizing also the decoupling of the squark contributions for the singlet-like Higgs boson A 2 .
We point out that the singlet-like Higgs boson H 3 and A 2 both approach SUSY masses thresholds with decreasing λ. Therefore we employ the lower bound of λ = 0.04, since for lower values of λ and thus larger masses of H 3 and A 2 we cannot guarantee the validity of the NLO SQCD contributions implemented in SusHi. Ref. [44] therefore assigned an additional theory uncertainty to the heavy SUSY masses expansion. In the decoupling regime we checked that the cross sections for H 1 , H 2 and A 1 coincide with the MSSM cross sections obtained for a mixing angle of α = −0.12347 with an accuracy of ∼ 10 −4 , which resembles the remaining singlet fraction of H 1 , H 2 and A 1 .

Theory uncertainties
In this section we shortly focus on theory uncertainties in the calculation of neutral Higgs boson production cross sections. Ref. [44] identified the most important theory uncertainties for the MSSM, which mostly apply to our discussion of the NMSSM as well. Apart from the well-known renormalization and factorization scale and PDF+α s uncertainties for cross sections at a proton-proton collider an additional uncertainty for gluon fusion cross section is the choice of a renormalization scheme for the bottom-quark Yukawa coupling, which is of particular relevance if the bottom-quark loop dominantly contributes. Secondly, the fact that NLO SQCD contributions are taken into account in an expansion of heavy SUSY masses induces an uncertainty, which grows for larger Higgs masses approaching SUSY particle masses thresholds. Thirdly and also relevant for bottom-quark annihilation are missing contributions in the resummation ∆ b , which induce an uncertainty, in particular in the limit ∆ b → −1. All of the above theory uncertainties as discussed in Ref. [44] apply to the NMSSM in a similar way. In contrast to the MSSM however phenomenological studies of the NMSSM focus on lower values of tan β, where both the uncertainty from the choice of the bottom-quark Yukawa coupling as well as the uncertainty induced from unknown contributions to ∆ b are of less importance. A detailed discussion in particular for the singlet-like CP-even and CP-odd Higgs boson is left for future work.
In the following we stick to the commonly studied renormalization and factorization scale uncertainties as well as the PDF+α s uncertainties. We present our results just for scenario S 1 , since no generically new features appear in other SUSY scenarios. We start with the scale uncertainty, where we follow the prescription employed in Refs. [6,44]. We thus consider seven combinations of renormalization and factorization scales defined as set C µ of pairs (µ R , µ F ) with µ R = {m φ /4, m φ /2, m φ } and µ F = {m φ /4, m φ /2, m φ } under the constraint 1/2 ≤ µ R /µ F ≤ 2 for gluon fusion. For bottom-quark annihilation the set is determined from The minimal and maximal cross sections are obtained according to which we present relative to the cross sections σ(µ 0 R , µ 0 F ) at the central scales µ 0 R and µ 0 F . They are µ 0 R = µ 0 F = m φ /2 for gluon fusion and µ 0 R = m φ and µ 0 F = m φ /4 for bottom-quark annihilation. We therefore define the relative uncertainties The scale uncertainties are shown for both the CP-even and CP-odd Higgs bosons in Fig Fig. 7 (a) and for A 1 in Fig. 7 (c).
For bottom-quark annihilation as shown in Fig. 7 (b) and (d) the scale uncertainty is mainly dependent on the Higgs mass, rather than the specific SUSY scenario. The large uncertainty for low Higgs masses reflects the need to move towards the four-flavour scheme (4FS) [90,91] in the description of the process.
As a last step we discuss PDF+α s uncertainties by applying the practical PDF4LHC recommendation [92,93] for the MSTW2008 [89] PDF sets in order to emphasize the findings of Ref. [44] for the case of the NMSSM. For this purpose we combine the results obtained with the 41 PDF sets of MSTW2008(n)nlo68cl with the α s uncertainties obtained by the PDF sets, which vary α s within the 68% confidence level interval. Fig. 8 shows the PDF+α s relative uncertainties ∆ ± PDF+αs with respect to the standard PDF+α s choice for the three CPeven Higgs bosons for gluon fusion and bottom-quark annihilation. The standard PDF+α s choice equals the zeroth PDF set of MSTW2008(n)nlo68cl together with the standard values α s = 0.120 at NLO and α s = 0.117 at NNLO QCD. Similar to the MSSM the uncertainties are mainly dependent on the Higgs mass and only slightly dependent on the specific SUSY scenario, even for the singlet-like Higgs boson. A very similar result applies to the CP-odd Higgs sector and is thus not explicitly presented. It therefore seems sufficient to take over the full relative PDF+α s uncertainties from a CP-even or CP-odd SM Higgs boson with the same mass, which is easily adjustable to future updates of the PDF4LHC recommendation. Taking into account the combination of the newest MMHT2014 [94], NNPDF 3.0 [95] and CT10 [96] PDF sets naturally results in larger PDF+α s uncertainties. However the simple recipe to obtain PDF+α s uncertainties just as a function of the Higgs mass is applicable for the combination of the PDF sets provided by the PDF fitting groups as well.

Conclusions
We presented accurate predictions for neutral Higgs production at proton colliders through gluon fusion and bottom-quark annihilation in the CP-conserving NMSSM. For gluon fusion we adapt the full NLO QCD and SQCD results from the MSSM to the NMSSM, based on an asymptotic expansion in heavy SUSY masses for squark and squark/quark/gluino twoloop contributions. Top-quark induced NNLO QCD contributions are added in the heavy top-quark effective theory. Electro-weak corrections to gluon fusion mediated through light quarks are taken into account and the resummation of sbottom contributions for large values of tan β can be translated from the MSSM to the NMSSM. The latter procedure also applies to bottom-quark annihilation. Our discussion comes along with an implementation of the neutral Higgs production cross section calculation in the code SusHi. The Higgs sector (obtained by an NMSSM spectrum generator) needs to be supplied through the SusHi input file. We shortly focused on the new features of the additional singlet-like CP-even or CP-odd Higgs boson for what concerns neutral Higgs production. Due to possible cancellations of quark induced contributions, squark and electro-weak corrections to gluon fusion can be of larger relevance than known in the MSSM, in particular for not too heavy third generation squark mass spectra. For a small singlet-doublet mixing term, which can be achieved by lowering the parameter λ, the singlet-like CP-even and -odd Higgs boson can both be decoupled from the remaining MSSMlike Higgs sector with regard to Higgs production. The renormalization and factorization scale uncertainties reflect the individual contributions to neutral Higgs production in case of gluon fusion, whereas scale uncertainties for bottom-quark annihilation as well as PDF+α s uncertainties for both production processes mainly remain a function of the Higgs mass.
We leave a more detailed investigation of theory uncertainties to future work. Moreover interesting for future studies is an expansion in a light Higgs mass rather than heavy SUSY masses for what concerns the inclusion of NLO and NNLO SQCD contributions, in particular since for pure singlet-like CP-odd Higgs bosons NNLO stop induced contributions are the first non-vanishing contributions to gluon fusion. Similarly a discussion of distributions and of the necessity of resummation for transverse momentum distributions is timely for the real NMSSM, but left for future work.

Acknowledgments
The author thanks Jonathan Gaunt, Robert Harlander, Hendrik Mantler and Pietro Slavich for very helpful comments on the manuscript. The author is moreover indebted to Pietro Slavich for help in the translation of the MSSM NLO SQCD corrections to the NMSSM and to Robert Harlander and Hendrik Mantler for their comments on the implementation of the NMSSM in SusHi. The author also thanks Kathrin Walz for help related to the renormalization of the stop and sbottom sectors within NMSSMCalc and Peter Drechsel for help with regard to the calculation of NMSSM Higgs boson masses. The author acknowledges support by the DFG through the Collaborative Research Center SFB 676 "Particles, Strings and the Early Universe".
A Formulas: Higgs-squark-squark couplings in the NMSSM In this section we present the squark couplings to the five neutral Higgs bosons φ of the NMSSM as implemented in the code SusHi. As denoted before the singlet component S does not couple to quarks, such that the couplings of the Higgs bosons to quarks can be taken over from the MSSM by replacing the mixing angle α and thus the projection on H 0 d and H 0 u by the proper mixing matrix elements R S /R P . Contrary the singlet component couples to squarks, for which we present the Feynman rules in the form with v = 1/ √ 2G F = v 2 d + v 2 u . The couplings g φ q,ij of squarks with indices {i, j} to CPeven Higgs bosons with k = {1, 2, 3} or CP-odd Higgs bosons with k = {1, 2} are subsequently presented in gauge eigenstates where in the CP-odd sector the prerotation with R G involving c β = cos β and s β = sin β is performed. They were obtained with the code MaCoR [97] and cross-checked against the formulas of Ref. [4] for what concerns the LO stop contributions to gluon fusion. The individual contributionsg Hk q,ij in the CP-even sector with k = {1, 2, 3} yield: In the CP-odd sector contributions with identical squark indicesg Ak q,ii do not exist. The remaining ones are given by: All occurrences of the soft-breaking parameters A t and A b were replaced by their relation to the squark mixing angles θ b and θ t . Trigonometric functions are abbreviated through s x = sin x, c x = cos x and t x = tan x.