---
abstract: 'The next-to-minimal supersymmetric standard model (NMSSM) more naturally accommodates a Higgs boson with a mass of approximately 125 GeV than the minimal supersymmetric standard model (MSSM). In this work, we assume that the next-to-lightest CP-even Higgs boson $h_2$ is the SM-like Higgs boson $h$, whereas the lightest CP-even Higgs boson $h_1$ is dominantly singlet-like. We discuss the $h_1h_1$, $h_2h_2$, and $h_1h_2$ pair production processes via gluon-gluon fusion at the LHC for an collision energy of 14 TeV, and we consider the cases in which one Higgs boson decays to $b\bar b$ and the other one decays to $\gamma\gamma$ or $\tau^+\tau^-$. We find that, for $m_{h_1} \lesssim$ 62 GeV, the cross section of the $gg \to h_1 h_1$ process is relatively large and maximally reaches 5400 fb, and the production rate of the $h_1h_1\to b\bar b \tau^+\tau^-$ final state can reach 1500 fb, which make the detection of this final state possible for future searches of an integrated luminosity of 300 and 3000 $fb^{-1}$. This is mainly due to the contributions from the resonant production process $pp\to h_2\to h_1h_1$ and the relatively large branching ratio of $h_1\to b\bar b$ and $h_1\to\tau^+\tau^-$. The cross sections of the $pp \to h_2h_2$ and $pp \to h_1 h_2$ production processes maximally reach 28 fb and 133 fb, respectively.'
author:
- 'Zhaoxia Heng, Xue Gong, Haijing Zhou'
title: 'Pair production of Higgs boson in NMSSM at the LHC with next-to-lightest CP-even Higgs boson being SM-like'
---

Introduction
============

In July 2012, the ATLAS and CMS collaborations at the CERN Large Hadron Collider (LHC) [@1207ATLAS-CMS] reported the observation of a new boson with a mass of approximately 125 GeV. After such discovery, which marked a milestone in particle physics, the extensive amount of data accumulated by both experiments [@LHC-experiment] have established that the properties of the new boson are consistent with those predicted by the standard model (SM). Since the discovery of the Higgs boson, massive work has been done to study its properties by measuring all the possible production channel cross sections and decay rates. In particular, it has been demonstrated that, among all the possible production modes, the double Higgs production is a sensitive probe for new physics [@hh-newphysics; @hh-1617; @hhyang; @Ellwanger:2013ova; @gghh-cao2013; @hh-susy].

At leading order in the SM, the Higgs pair production at the LHC happens through the gluon–gluon fusion $gg\to hh$ process [@Djouadi:1999rca; @gghh]. The SM cross section of such a process receives contributions from loops mediated by heavy quarks, which induce box and triangle diagrams ($gg\to h^*\to hh$). The amplitudes of these two types of diagrams interfere with each other destructively, thus leading to a rather small cross section. In new physics models, such as low energy supersymmetric models, both box and triangle diagrams receive additional contributions from the loops mediated by the third generation of squarks. Moreover, the triangle diagrams are correlated with the trilinear Higgs coupling, which plays an important role in the reconstruction of the Higgs potential [@SM-NLO-35fb]. Therefore, the $gg\to hh$ process is particularly sensitive to new physics contributions.

Among new physics models, the next-to-minimal supersymmetric standard model (NMSSM) has attracted considerable attentions [@Ellwanger:2009dp], as it more naturally accommodates a Higgs boson with a mass near 125 GeV than the minimal supersymmetric standard model (MSSM). Compared with the MSSM, the NMSSM introduces an additional singlet Higgs field $\hat S$. Because of the coupling between singlet and doublet Higgs fields in the superpotential, the mass of the SM-like Higgs boson can be enhanced at the tree-level. Furthermore, the mixing between singlet and doublet Higgs fields can also increase the mass of the Higgs boson when the next-to-lightest CP-even Higgs boson corresponds to the SM-like Higgs boson [@Cao:2012fz; @Kang:2012sy; @Cao-natural]. These two contributions can easily increase the mass of the Higgs boson to 125 GeV without any large radiative corrections. In this work, we only consider the scenario in which the next-to-lightest CP-even Higgs boson acts as the SM-like Higgs boson $h$.

Due to the presence of the additional superfield $\hat S$, the NMSSM contains three CP-even Higgs particles $h_i\ (i=1,2,3)$, which respect the convention $m_{h_1}<m_{h_2}<m_{h_3}$. In this work, we assume that $h_2$ acts as the SM-like Higgs boson $h$, and $h_1$ is dominantly singlet-like considering the constraints for Higgs searches from LEP, Tevatron and LHC. We study the main di-Higgs production process at the LHC, the gluon–gluon fusion production of $h_1h_1$, $h_2h_2$, and $h_1h_2$ Higgs pairs. At present, the possible most promising detection channels for Higgs pair production at the LHC are $hh\to b\bar b\gamma\gamma$, $hh\to b\bar b \tau^+\tau^-$ and $hh \to b\bar b W^+W^-$ [@hhsearch]. However, due to the weak interaction between the lightest CP-even Higgs $h_1$ and $W^+W^-$, we will focus on the $b\bar b\gamma\gamma$ and $b\bar b \tau^+\tau^-$ channels.

This work is organized as follows. In section II we briefly introduce the Higgs sector of the NMSSM. In section III we scan over the NMSSM parameter space by considering various experimental constraints; then we calculate the cross sections of $h_1h_1$, $h_2h_2$, and $h_1h_2$ di-Higgs gluon-fusion production processes and we also discuss the production rates for the $b\bar b\gamma\gamma$ and $b\bar b \tau^+\tau^-$ final states. Finally, we present our conclusions in section IV.

Higgs sector in NMSSM
=====================

As in the MSSM, the NMSSM contains two doublet Higgs superfields $\hat H_u$ and $\hat H_d$; however, in addition to the MSSM, the NMSSM introduces one gauge singlet superfield $\hat S$. As a consequence, the Higgs sector has a richer structure in the NMSSM than in the MSSM. To avoid the occurrence of dimensionful parameters in the superpotential, a discrete $Z_3$-symmetry is imposed. The corresponding superpotential and soft-breaking terms in the Higgs sector of the $Z_3$-invariant NMSSM are given by Equations (1) and (2) [@Ellwanger:2009dp; @NMSSM]: $$\begin{aligned}
 \label{NMSSM}
 W_{\rm NMSSM}&=&Y_u\hat{Q}\cdot\hat{H_u}\hat{U} -Y_d\hat{Q}\cdot\hat{H_d}\hat{D} -Y_e\hat{L}\cdot\hat{H_d}\hat{E}\nonumber\\
 &&+ \lambda\hat{H_u} \cdot \hat{H_d} \hat{S}
 + \frac{1}{3}\kappa \hat{S^3},\\
V_{\rm soft}^{\rm NMSSM}&=&\tilde m_u^2|H_u|^2 + \tilde m_d^2|H_d|^2
+ \tilde m_S^2|S|^2 \nonumber\\
&& +(A_\lambda \lambda SH_u\cdot H_d
+\frac{A_\kappa}{3}\kappa S^3 + h.c.).\end{aligned}$$ where $\hat{Q}$, $\hat{U}$, $\hat{D}$, $\hat{L}$, and $\hat{E}$ are chiral superfields with $Y_u, Y_d, Y_e$ being their Yukawa coupling coefficients; $\lambda$ and $\kappa$ denote dimensionless coupling coefficients; and $\tilde{m}_{u}$, $\tilde{m}_{d}$, $\tilde{m}_{S}$, $A_\lambda$, and $A_\kappa$ are the soft-breaking parameters. When the electroweak symmetry is broken, the parameters $\tilde{m}_{u}$, $\tilde{m}_{d}$, and $\tilde{m}_{S}$ are expressed in terms of $m_Z$, $\tan\beta\equiv v_u/v_d$, and $\mu\equiv \lambda v_s$ with $v_u$, $v_d$, and $v_s$ being the vacuum expectation values (VEV) of the Higgs fields $H_u$, $H_d$, and $S$. Thus, at the tree level, the Higgs sector of the NMSSM can be represented by six independent parameters given by Equation (3): $$\begin{aligned}
\label{parameter}
\lambda, \quad \kappa, \quad \tan \beta, \quad \mu, \quad  M_A, \quad A_\kappa,\end{aligned}$$ with $M_A^2= \frac{2 \mu (A_\lambda + \kappa v_s)}{\sin 2 \beta}$.

In the NMSSM, $H_u$, $H_d$, and $S$ are given by Equation (4): $$\begin{aligned}
H_u = \left ( \begin{array}{c} H_u^+ \\
       v_u +\frac{ \phi_u + i \varphi_u}{\sqrt{2}}
        \end{array} \right),~~
H_d & =& \left ( \begin{array}{c}
             v_d + \frac{\phi_d + i \varphi_d}{\sqrt{2}}\\
             H_d^- \end{array} \right),~~
S  =  v_s + \frac{1}{\sqrt{2}} \left(\sigma + i \xi \right).\end{aligned}$$ By diagonalizing the corresponding mass matrices, we get the Higgs mass eigenstates as in Equation (5): $$\begin{aligned}
 \left( \begin{array}{c} h_1 \\
h_2 \\ h_3 \end{array} \right) = U^H \left( \begin{array}{c} \phi_u
\\ \phi_d\\ \sigma\end{array} \right),~ \left(\begin{array}{c} a_1\\
a_2\\ G^0 \end{array} \right) = U^A \left(\begin{array}{c} \varphi_u
\\ \varphi_d \\ \xi \end{array} \right),~ \left(\begin{array}{c} H^+
\\G^+ \end{array}  \right) = U^C \left(\begin{array}{c}H_u^+\\ H_d^+
\end{array} \right).  \label{rotation}\end{aligned}$$ Here, $h_i$ (with $i=1,2,3$) denotes the physical CP-even neutral scalars with $m_{h_1}<m_{h_2}<m_{h_3}$; $a_i$ (with $i=1,2$) denotes the physical CP-odd neutral scalars with $m_{a_1}<m_{a_2}$; $H^+$ is the physical charged Higgs boson; and $G^0$ and $G^+$ are the Goldstone bosons absorbed by $Z$ and $W$ bosons, respectively.

Because of the coupling $\lambda\hat{H_u} \cdot \hat{H_d} \hat{S}$ in the superpotential, the mass of the SM-like Higgs boson $h$ at the tree level has the following dependencies: $$m^2_{h,tree}= m_Z^2\cos^22\beta+\lambda^2v^2\sin^22\beta.$$ As shown by Equation (6), the SM-like Higgs boson ($h$) mass increases as $\lambda$ increases and $\tan\beta$ decreases. Moreover, the mixing between singlet and doublet Higgs fields can also significantly alter the $h$ mass. In particular, the mixing will pull down the mass when $h_1$ acts as $h$, whereas it will push up the mass when $h_2$ acts as $h$ [@Cao:2012fz; @Kang:2012sy]. Therefore, when $h_2$ plays the role of the SM-like Higgs boson, both the additional tree-level contributions and the effect of Higgs doublet-singlet mixing can increase the $h$ mass; such phenomenon makes the large radiative corrections from top-squark loops unnecessary to predict a Higgs boson with a mass near 125 GeV.

Calculations and numerical results
==================================

Scanning strategy of the parameter space
----------------------------------------

We perform a comprehensive scan over the parameter space of the NMSSM by using the package NMSSMTools [@NMSSMTools]. We constrain the mass of the gluino and the soft breaking parameters in the first two generations of the squark sector to be equal to 2 TeV. Moreover, we assume a common value denoted by $m_{\tilde{l}}$ for all soft breaking parameters in the slepton sector. The parameter space considered in the scan is represented by Equation (7): $$\begin{aligned}
\label{NMSSM-scan}
&& 0 <\lambda,\kappa\leq 0.75, \quad  2 \leq \tan{\beta} \leq 60, 100{\rm ~GeV}\leq \mu \leq 1 {\rm ~TeV},
     \quad 50 {\rm ~GeV}\leq M_A \leq 2 {\rm ~TeV}, \nonumber \\
&&  |A_{\kappa}| \leq 2 {\rm TeV},
\quad  100{\rm ~GeV}\leq m_{\tilde{l}} \leq 1 {\rm ~TeV}, 100{\rm ~GeV}\leq M_{Q_3}\leq 2 {\rm ~TeV}, \nonumber\\
&& 100{\rm ~GeV}\leq M_{U_3}(=M_{D_3}) \leq 2 {\rm ~TeV}, |A_{t}(=A_{b})|\leq {\rm min}(3 \sqrt{M_{Q_3}^2 + M_{U_3}^2}, 5 {\rm TeV}), \nonumber\\
&& 20 {\rm GeV} \leq M_1 \leq 500 {\rm GeV},100 {\rm GeV} \leq M_2 \leq 1 {\rm TeV}.\end{aligned}$$ In addition to the constraints implemented in the package NMSSMTools we also consider the following:

-   Constraints from the direct searches for the Higgs boson at LEP, Tevatron, and LHC. We adopt the package HiggsSignals for the fit of the 125 GeV Higgs data [@HiggsSignal] and the package HiggsBounds for the search of non-standard Higgs bosons at colliders [@HiggsBounds].

-   Constraints from direct searches of sparticles at the LHC RunI and Run-II. We use the packages FastLim [@FastLim] and SModelS [@SModelS] to implement such constraints. In particular, the results of the experimental searches for gluino and squarks are implemented in the package FastLim, whereas those of sleptons and electroweakinos searches are implemented in the SModelS package. Efficiencies or upper bounds on the sparticle productions can be obtained to limit the parameter space of the NMSSM.

-   Constraints from the searches for electroweakinos and top-squarks at the LHC Run-II. The surviving samples that satisfy the constraints from FastLim and SModelS are further tested by detailed Monte Carlo simulations. In details, we use MadGraph/MadEvent [@madgraph] to simulate events at the parton level and Pythia [@pythia] to do parton shower and hadronization, then we put the events into Delphes [@delphes] to do detecor simulation. Finally, we implement the ATLAS and CMS searches within CheckMATE [@checkmate] to decide whether each parameter point in our scenaio is allowed or not.

Numerical results
-----------------

In our numerical calculation we take $m_t$ = 173.2GeV, $m_b$ = 4.18GeV, $m_Z$ = 91.19GeV [@PDG] and take the renormalization and factorization scales to be the invariant mass of the Higgs pair. We select the samples for 122 GeV $\leq m_h \leq$ 128 GeV and require the next-to-lightest CP-even Higgs boson ($h_2$) to be the SM-like Higgs boson $h$.

![Feynman diagrams of Higgs pair production in the NMSSM with a,b=1,2 and c=1,2,3.[]{data-label="fig1"}](fig1.eps){width="15cm"}

According to the NMSSM scenario, the Higgs pair production processes possibly available at the LHC happen via gluon–gluon fusion, as shown in Fig. \[fig1\]. Differently from the SM, in the NMSSM the gluon–gluon fusion receives additional contributions from the loops mediated by the third generation of squarks. We checked that we can reproduce the results in the SM and MSSM presented in [@Djouadi:1999rca] and [@gghh-MSSM]. In the numerical calculations, for an LHC collision energy of 14 TeV, we find that the cross section of the $gg\to hh$ process in the SM scenario is 18.7 fb at $m_h =$ 125 GeV.

Figure \[fig2\] shows the cross sections of the $h_1h_1$ (left), $h_2h_2$ (middle), and $h_1h_2$ (right) pair production processes. The cross section of the process $pp \to h_1 h_1$ (left panel) assumes values mainly divided into two regions: for $m_{h_1} \lesssim$ 62 GeV, the cross section is relatively large and it can maximally reach 5400 fb. This is because the rare decay $h_2\to h_1h_1$ is open when $m_{h_1} \lesssim$ 62 GeV. As can be clearly seen from Fig. \[fig3\], which shows the cross section of the $pp\to h_1h_1$ process as a function of the branching ratio of the decay $h_2\to h_1h_1$, the $pp\to h_1h_1$ cross section increases as the $h_2\to h_1h_1$ branching ratio increases.

![Cross sections of Higgs pair production processes in the final states $h_1h_1$, $h_2h_2$, and $h_1h_2$.[]{data-label="fig2"}](fig21.eps "fig:"){width="6.3cm"} ![Cross sections of Higgs pair production processes in the final states $h_1h_1$, $h_2h_2$, and $h_1h_2$.[]{data-label="fig2"}](fig22.eps "fig:"){width="6.3cm"} ![Cross sections of Higgs pair production processes in the final states $h_1h_1$, $h_2h_2$, and $h_1h_2$.[]{data-label="fig2"}](fig23.eps "fig:"){width="6.3cm"}

![Cross section of the Higgs pair production process in the final state $h_1h_1$ as a function of the $h_2\to h_1h_1$ branching ratio.[]{data-label="fig3"}](fig3.eps){width="9cm"}

The plot in the middle panel of Fig. \[fig2\] indicates that the cross section of the $pp\to h_2h_2$ process can maximally reach 28 fb, which is 1.5 times larger than the SM prediction. The enhancement mainly originates from the contributions of the diagrams (3), (4), and (5) in Fig. \[fig1\]. For lighter $m_{\tilde t_1}$ and larger values of $A_t$ the cross section is larger, which is in agreement with the results reported in [@gghh-cao2013].

From the right panel of Fig. \[fig2\] it is visible that, for the majority of the surviving samples, the cross section of $pp\to h_1h_2$ is less than 10 fb; this mainly because $h_1$ is dominantly singlet-like. However, for a few surviving samples, the cross section can exceed 100 fb and reach 133 fb at maximum. In these cases, the heaviest CP-even Higgs $h_3$ is relatively light with a mass of approximately 500 GeV, and $h_3$ is produced on-shell with the branching ratio to $h_1h_2$ being about 10%. It is important to notice that the mass of the heavy Higgs boson near 500 GeV for $\tan\beta\lesssim 5$ in the MSSM is not excluded by experimental results based on the LHC Run-II data [@heavy; @Higgs-ATLAS; @heavy; @Higgs-CMS].

Since the discovery potential of the Higgs boson depends on the decay spectrum of the Higgs boson, we now consider the $h_1h_1$, $h_2h_2$, and $h_1h_2$ di-Higgs production processes with one Higgs boson decaying to $b\bar b$ and the other one decaying to either $\gamma\gamma$ or $\tau^+\tau^-$. In Fig. \[fig4\] and Fig. \[fig5\] we show the production rates of the final states $b\bar b \gamma\gamma$ and $b\bar b \tau^+\tau^-$, respectively. As it can be seen, the production rate of the process $h_1h_1\to b\bar b \tau^+\tau^-$ is relatively large for $m_{h_1}\lesssim$ 62 GeV, and it can maximally reach 1500 fb. This is because, for most of the surviving samples, the $h_1\to b\bar b$ and $h_1\to\tau^+\tau^-$ branching ratios are close to 90% and 8%, respectively.

![Production rates of Higgs pair production processes in the final states $h_1h_1\to b\bar b \gamma\gamma$, $h_2h_2\to b\bar b \gamma\gamma$, $h_1h_2\to b\bar b \gamma\gamma$, and $h_2h_1\to b\bar b \gamma\gamma$.[]{data-label="fig4"}](fig41.eps "fig:"){width="8.5cm"} ![Production rates of Higgs pair production processes in the final states $h_1h_1\to b\bar b \gamma\gamma$, $h_2h_2\to b\bar b \gamma\gamma$, $h_1h_2\to b\bar b \gamma\gamma$, and $h_2h_1\to b\bar b \gamma\gamma$.[]{data-label="fig4"}](fig42.eps "fig:"){width="8.5cm"} ![Production rates of Higgs pair production processes in the final states $h_1h_1\to b\bar b \gamma\gamma$, $h_2h_2\to b\bar b \gamma\gamma$, $h_1h_2\to b\bar b \gamma\gamma$, and $h_2h_1\to b\bar b \gamma\gamma$.[]{data-label="fig4"}](fig43.eps "fig:"){width="8.5cm"} ![Production rates of Higgs pair production processes in the final states $h_1h_1\to b\bar b \gamma\gamma$, $h_2h_2\to b\bar b \gamma\gamma$, $h_1h_2\to b\bar b \gamma\gamma$, and $h_2h_1\to b\bar b \gamma\gamma$.[]{data-label="fig4"}](fig44.eps "fig:"){width="8.5cm"}

![Production rates of Higgs pair production processes in the final states $h_1h_1\to b\bar b \tau^+\tau^-$, $h_2h_2\to b\bar b \tau^+\tau^-$, $h_1h_2\to b\bar b \tau^+\tau^-$, and $h_2h_1\to b\bar b \tau^+\tau^-$.[]{data-label="fig5"}](fig51.eps "fig:"){width="8.5cm"} ![Production rates of Higgs pair production processes in the final states $h_1h_1\to b\bar b \tau^+\tau^-$, $h_2h_2\to b\bar b \tau^+\tau^-$, $h_1h_2\to b\bar b \tau^+\tau^-$, and $h_2h_1\to b\bar b \tau^+\tau^-$.[]{data-label="fig5"}](fig52.eps "fig:"){width="8.5cm"} ![Production rates of Higgs pair production processes in the final states $h_1h_1\to b\bar b \tau^+\tau^-$, $h_2h_2\to b\bar b \tau^+\tau^-$, $h_1h_2\to b\bar b \tau^+\tau^-$, and $h_2h_1\to b\bar b \tau^+\tau^-$.[]{data-label="fig5"}](fig53.eps "fig:"){width="8.5cm"} ![Production rates of Higgs pair production processes in the final states $h_1h_1\to b\bar b \tau^+\tau^-$, $h_2h_2\to b\bar b \tau^+\tau^-$, $h_1h_2\to b\bar b \tau^+\tau^-$, and $h_2h_1\to b\bar b \tau^+\tau^-$.[]{data-label="fig5"}](fig54.eps "fig:"){width="8.5cm"}

Both the ATLAS and CMS collaborations performed searches for 125GeV Higgs pair production in the final states $b\bar b b\bar b$, $b\bar b\tau^+\tau^-$ and $b\bar b \gamma\gamma$. Based on the data collected at $\sqrt{s}$ = 13 TeV with an integrated luminosity of 35.9 $fb^{-1}$, the observed 95% CL upper limit for $b\bar b b\bar b$ final state is 147fb corresponding to 13 times the SM prediction [@bbbb-atlas], for $b\bar b\tau^+\tau^-$ final state is 75.4fb corresponding to 30 times the SM prediction [@bbtau-CMS] and for $b\bar b \gamma\gamma$ final state is 1.67fb corresponding to 19.2 times the SM prediction [@bbgamma-CMS]. For the present integrated luminosity of LHC, the Higgs pair production in the considered scenario of NMSSM cannot be directly detected. Due to the large production rate of the process $pp\to h_1h_1\to b\bar b\tau^+\tau^-$, for future searches of an integrated luminosity of 300 and 3000 $fb^{-1}$, Higgs pair production for $h_1h_1\to b\bar b\tau^+\tau^-$ final state may be detected.

Conclusion
==========

Because of the mixing effect of Higgs doublet-singlet fields, and the coupling effect between Higgs doublet and singlet fields in the superpotential, the NMSSM more naturally accommodates a Higgs boson with a mass of approximately 125 GeV. We assume that the next-to-lightest CP-even Higgs boson $h_2$ plays the role of the SM-like Higgs boson. In this scenario the lightest CP-even Higgs boson $h_1$ is dominantly singlet-like and it can be much lighter than 125 GeV. In this work, we discussed the $h_1h_1$, $h_2h_2$, and $h_1h_2$ pair production processes via gluon-gluon fusion at the LHC for an collision energy of 14 TeV. In addition, we studied the production rate of such di-Higgs events in which one Higgs boson decays to $b\bar b$ and the other one to either $\gamma\gamma$ or $\tau^+\tau^-$. We found that, for $m_{h_1}\lesssim$ 62 GeV, the cross section of the $gg \to h_1h_1$ production process is relatively large and maximally reaches 5400 fb, and that the production rate of the $h_1h_1\to b\bar b \tau^+\tau^-$ final state can reach 1500 fb. These results are mainly due to the contributions from the resonant production process $pp\to h_2\to h_1h_1$and the relatively large branching ratio of the $h_1\to b\bar b$ and $h_1\to\tau^+\tau^-$ decays. The cross sections of the $pp \to h_2h_2$ and $pp \to h_1 h_2$ production processes maximally reach 28 fb and 133 fb, respectively. Therefore, for future searches of an integrated luminosity of 300 and 3000 $fb^{-1}$, Higgs pair production for $h_1h_1\to b\bar b\tau^+\tau^-$ final state may be detected.

Acknowledgement {#acknowledgement .unnumbered}
===============

We thank Prof. Junjie Cao and Dr. Liangliang Shang for helpful discussions in the revision of the work. This work was supported in part by the National Natural Science Foundation of China (NNSFC) under grant No. 11305050 and 11705048.

[99]{}

G. Aad [*et al.*]{} \[ATLAS Collaboration\], Phys. Lett. B [**716**]{}, 1 (2012); S. Chatrchyan [*et al.*]{} \[CMS Collaboration\], Phys. Lett. B [**716**]{}, 30 (2012). G. Aad [*et al.*]{} \[ATLAS Collaboration\], Eur. Phys. J. C [**76**]{}, no. 1, 6 (2016) \[arXiv:1507.04548 \[hep-ex\]\]; V. Khachatryan [*et al.*]{} \[CMS Collaboration\], Eur. Phys. J. C [**75**]{}, no. 5, 212 (2015) \[arXiv:1412.8662 \[hep-ex\]\]; S. Chatrchyan [*et al.*]{} \[CMS Collaboration\], Phys. Rev. Lett.  [**110**]{}, no. 8, 081803 (2013) \[arXiv:1212.6639 \[hep-ex\]\]; G. Aad [*et al.*]{} \[ATLAS Collaboration\], Phys. Lett. B [**726**]{}, 120 (2013) \[arXiv:1307.1432 \[hep-ex\]\]. M. J. Dolan, C. Englert and M. Spannowsky, Phys. Rev. D [**87**]{}, no. 5, 055002 (2013) \[arXiv:1210.8166 \[hep-ph\]\]; S. Dawson, A. Ismail and I. Low, Phys. Rev. D [**91**]{}, no. 11, 115008 (2015) \[arXiv:1504.05596 \[hep-ph\]\]; H. J. He, J. Ren and W. Yao, Phys. Rev. D [**93**]{}, no. 1, 015003 (2016) \[arXiv:1506.03302 \[hep-ph\]\]; G. D. Kribs and A. Martin, Phys. Rev. D [**86**]{}, 095023 (2012) \[arXiv:1207.4496 \[hep-ph\]\]. M. Carena, Z. Liu and M. Riembau, arXiv:1801.00794 \[hep-ph\]; J. H. Kim, Y. Sakaki and M. Son, arXiv:1801.06093 \[hep-ph\]; L. Bian, N. Chen and Y. Jiang, Int. J. Mod. Phys. A [**32**]{}, no. 34, 1746002 (2017) \[arXiv:1712.01632 \[hep-ph\]\]; J. Ren, R. Q. Xiao, M. Zhou, Y. Fang, H. J. He and W. Yao, arXiv:1706.05980 \[hep-ph\]; I. M. Lewis and M. Sullivan, Phys. Rev. D [**96**]{}, no. 3, 035037 (2017) \[arXiv:1701.08774 \[hep-ph\]\]; S. Dawson and M. Sullivan, arXiv:1711.06683 \[hep-ph\]; K. Nakamura, K. Nishiwaki, K. y. Oda, S. C. Park and Y. Yamamoto, Eur. Phys. J. C [**77**]{}, no. 5, 273 (2017) \[arXiv:1701.06137 \[hep-ph\]\]. C. Han, X. Ji, L. Wu, P. Wu and J. M. Yang, JHEP [**1404**]{}, 003 (2014) \[arXiv:1307.3790 \[hep-ph\]\]; X. F. Han, L. Wang and J. M. Yang, Nucl. Phys. B [**825**]{}, 222 (2010) \[arXiv:0908.1827 \[hep-ph\]\]; X. F. Han, L. Wang and J. M. Yang, Mod. Phys. Lett. A [**31**]{}, no. 31, 1650178 (2016) \[arXiv:1509.02453 \[hep-ph\]\]; J. Cao, D. Li, L. Shang, P. Wu and Y. Zhang, JHEP [**1412**]{}, 026 (2014) \[arXiv:1409.8431 \[hep-ph\]\]; Z. Heng, L. Shang, Y. Zhang and J. Zhu, JHEP [**1402**]{}, 083 (2014) \[arXiv:1312.4260 \[hep-ph\]\]. U. Ellwanger, JHEP [**1308**]{}, 077 (2013) \[arXiv:1306.5541 \[hep-ph\]\]. J. Cao, Z. Heng, L. Shang, P. Wan and J. M. Yang, JHEP [**1304**]{}, 134 (2013) \[arXiv:1301.6437 \[hep-ph\]\]. P. Huang, A. Joglekar, M. Li and C. E. M. Wagner, arXiv:1711.05743 \[hep-ph\].

A. Djouadi, W. Kilian, M. Muhlleitner and P. M. Zerwas, Eur. Phys. J. C [**10**]{}, 45 (1999) \[hep-ph/9904287\]. T. Plehn, M. Spira and P. M. Zerwas, Nucl. Phys. B [**479**]{}, 46 (1996) Erratum: \[Nucl. Phys. B [**531**]{}, 655 (1998)\] \[hep-ph/9603205\]; E. W. N. Glover and J. J. van der Bij, Nucl. Phys. B [**309**]{}, 282 (1988). J. Baglio [*et al.*]{} JHEP [**1304**]{}, 151 (2013) \[arXiv:1212.5581 \[hep-ph\]\]; D. Y. Shao, C. S. Li, H. T. Li and J. Wang, JHEP [**1307**]{}, 169 (2013) \[arXiv:1301.1245 \[hep-ph\]\]. U. Ellwanger, C. Hugonie and A. M. Teixeira, Phys. Rept.  [**496**]{} (2010) 1 \[arXiv:0910.1785\].

J. J. Cao, Z. X. Heng, J. M. Yang, Y. M. Zhang and J. Y. Zhu, JHEP [**1203**]{}, 086 (2012) \[arXiv:1202.5821 \[hep-ph\]\]. Z. Kang, J. Li and T. Li, JHEP [**1211**]{}, 024 (2012) \[arXiv:1201.5305 \[hep-ph\]\]. J. Cao, Y. He, L. Shang, W. Su and Y. Zhang, JHEP [**1608**]{}, 037 (2016) \[arXiv:1606.04416 \[hep-ph\]\]; J. Cao, Y. He, L. Shang, W. Su, P. Wu and Y. Zhang, JHEP [**1610**]{}, 136 (2016) \[arXiv:1609.00204 \[hep-ph\]\]. U. Baur, T. Plehn and D. L. Rainwater, Phys. Rev. D [**69**]{}, 053004 (2004) \[hep-ph/0310056\]; A. Papaefstathiou, L. L. Yang and J. Zurita, Phys. Rev. D [**87**]{}, no. 1, 011301 (2013) \[arXiv:1209.1489 \[hep-ph\]\]; T. Huang [*et al.*]{}, Phys. Rev. D [**96**]{}, no. 3, 035007 (2017) \[arXiv:1701.04442 \[hep-ph\]\]; A. Adhikary, S. Banerjee, R. K. Barman, B. Bhattacherjee and S. Niyogi, arXiv:1712.05346 \[hep-ph\]. M. Maniatis, Int. J. Mod. Phys. A [**25**]{}, 3505 (2010) \[arXiv:0906.0777 \[hep-ph\]\]. U. Ellwanger, J. F. Gunion and C. Hugonie, JHEP [**0502**]{}, 066 (2005) \[hep-ph/0406215\]; U. Ellwanger and C. Hugonie, Comput. Phys. Commun.  [**175**]{}, 290 (2006) \[hep-ph/0508022\]. P. Bechtle, S. Heinemeyer, O. Stal, T. Stefaniak and G. Weiglein, Eur. Phys. J. C [**74**]{} (2014) no.2, 2711 \[arXiv:1305.1933 \[hep-ph\]\]; JHEP [**1411**]{}, 039 (2014) \[arXiv:1403.1582 \[hep-ph\]\]; O. Stal and T. Stefaniak, PoS EPS [**-HEP2013**]{}, 314 (2013) \[arXiv:1310.4039 \[hep-ph\]\]. P. Bechtle, O. Brein, S. Heinemeyer, G. Weiglein and K. E. Williams, Comput. Phys. Commun.  [**181**]{}, 138 (2010) \[arXiv:0811.4169 \[hep-ph\]\]; Comput. Phys. Commun.  [**182**]{}, 2605 (2011) \[arXiv:1102.1898 \[hep-ph\]\]. M. Papucci, K. Sakurai, A. Weiler and L. Zeune, Eur. Phys. J. C [**74**]{}, no. 11, 3163 (2014) \[arXiv:1402.0492 \[hep-ph\]\]. S. Kraml, S. Kulkarni, U. Laa, A. Lessa, W. Magerl, D. Proschofsky-Spindler and W. Waltenberger, Eur. Phys. J. C [**74**]{}, 2868 (2014) \[arXiv:1312.4175 \[hep-ph\]\]. J. Alwall [*et al.*]{}, JHEP [**1407**]{}, 079 (2014) \[arXiv:1405.0301 \[hep-ph\]\]; J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, JHEP [**1106**]{}, 128 (2011) \[arXiv:1106.0522 \[hep-ph\]\].

T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP [**0605**]{}, 026 (2006) \[hep-ph/0603175\].

J. de Favereau [*et al.*]{} \[DELPHES 3 Collaboration\], JHEP [**1402**]{}, 057 (2014) \[arXiv:1307.6346 \[hep-ex\]\].

D. Dercks, N. Desai, J. S. Kim, K. Rolbiecki, J. Tattersall and T. Weber, Comput. Phys. Commun.  [**221**]{}, 383 (2017)\[arXiv:1611.09856 \[hep-ph\]\]; J. S. Kim, D. Schmeier, J. Tattersall and K. Rolbiecki, Comput. Phys. Commun.  [**196**]{}, 535 (2015) \[arXiv:1503.01123 \[hep-ph\]\]; M. Drees, H. Dreiner, D. Schmeier, J. Tattersall and J. S. Kim, Comput. Phys. Commun.  [**187**]{}, 227 (2015) \[arXiv:1312.2591 \[hep-ph\]\].

C. Patrignani [*et al.*]{} \[Particle Data Group\], Chin. Phys. C [**40**]{}, no. 10, 100001 (2016). A. Belyaev, M. Drees, O. J. P. Eboli, J. K. Mizukoshi and S. F. Novaes, Phys. Rev. D [**60**]{}, 075008 (1999) \[hep-ph/9905266\]. M. Aaboud [*et al.*]{} \[ATLAS Collaboration\], JHEP [**1801**]{}, 055 (2018)\[arXiv:1709.07242 \[hep-ex\]\].

CMS Collaboration, CMS PAS HIG-17-020.

M. Aaboud [*et al.*]{} \[ATLAS Collaboration\], arXiv:1804.06174 \[hep-ex\].

A. M. Sirunyan [*et al.*]{} \[CMS Collaboration\], Phys. Lett. B [**778**]{}, 101 (2018) \[arXiv:1707.02909 \[hep-ex\]\].

CMS Collaboration \[CMS Collaboration\], CMS-PAS-HIG-17-008.
