Higgs pair production with SUSY QCD correction: revisited under current experimental constraints

We consider the current experimental constraints on the parameter space of the MSSM and NMSSM. Then in the allowed parameter space we examine the Higgs pair production at the 14 TeV LHC via $b\bar{b}\to hh$ ($h$ is the 125 GeV SM-like Higg boson) with one-loop SUSY QCD correction and compare it with the production via $gg\to hh$. We obtain the following observations: (i) For the MSSM the production rate of $b\bar{b} \to hh$ can reach 50 fb and thus can be competitive with $gg \to hh$, while for the NMSSM $b\bar{b} \to hh$ has a much smaller rate than $gg \to hh$ due to the suppression of the $hb\bar{b}$ coupling; (ii) The SUSY-QCD correction to $b\bar{b} \to hh$ is sizable, which can reach $45\%$ for the MSSM and $15\%$ for the NMSSM within the $1\sigma$ region of the Higgs data; (iii) In the heavy SUSY limit (all soft mass parameters become heavy), the SUSY effects decouple rather slowly from the Higgs pair production (especially the $gg\to hh$ process), which, for $M_{\rm SUSY}=5$ TeV and $m_A<1$ TeV, can enhance the production rate by a factor of 1.5 and 1.3 for the MSSM and NMSSM, respectively. So, the Higgs pair production may be helpful for unraveling the effects of heavy SUSY.


I. INTRODUCTION
The discovery of a Higgs boson at around 125 GeV has been announced by the ATLAS and CMS collaborations [1]. Up to now, the measurements of the Higgs boson properties are in good agreement with the Standard Model (SM) predictions except for the enhanced diphoton rate σ/σ SM = 1.65 +0. 34 −0.30 reported by the ATLAS collaboration. The future precise measurements will further test the SM and allow for a probe for new physics like supersymmetry (SUSY) which is a promising framework to accommodate such a 125 GeV Higgs boson [2][3][4][5][6]. Therefore, the intensive studies of the Higgs productions and decays are very important and urgent.
Among the productions of the Higgs boson at the LHC, the pair production is a rare process but quite important since it can be used to measure the Higgs self-couplings [7].
On the experimental side, the discovery potential of Higgs pair signal at the LHC has been studied by analyzing the decay channels hh → bbγγ/bbµ + µ − [8]. Recently, the jet substructure technique was applied to the Higgs pair production in the boosted final states [9], such as hh → bbτ + τ − /bbW + W − [10][11][12], which was found to be powerful in observing the events at the 14 TeV LHC with 600 fb −1 integrated luminosity [12]. On the theoretical side, in the SM the main pair production mechanism is found to be the gluon fusion gg → hh via heavy quark loops [13,14]. Numerous studies have also been performed for Higgs pair production in new physics models [15][16][17][18][19][20][21][22]. Note that although the bottom quark annilation bb → hh has a much smaller rate than the gluon fusion process in the SM [23,24], it can be significantly enhanced via the enlarged hbb coupling in new physics models like the Minimal Supersymmetric Standard Model (MSSM) [25].
In this work, we revisit the Higgs pair production in SUSY for two reasons. One is that the sizable SUSY-QCD correction must be considered for bb → hh, which has been presented in the MSSM but not in the NMSSM [26,27]. The other is that the studies should be updated by using the latest experimental constraints including the recent LHC Higgs data, the LHCb B s → µ + µ − data and the Planck dark matter relic density. It is also notable that the masses of the third generation sparticles involved in the SUSY-QCD correction to bb → hh have been pushed up to a few hundred GeV by the LHC direct searches [28]. So the size of such a correction will be quite different from the previous results in the literature [25,29]. This paper is organized as follows. In Section II we briefly review the Higgs sectors in the MSSM and NMSSM and give a description of the analytic calculation of the SUSY-QCD correction. Then in Section III we present the numerical results of Higgs pair production at the LHC and discuss the SUSY-QCD residual effects in the heavy sparticle limit. Finally, we draw the conclusion in Section IV.

II. A DESCRIPTION OF MODELS AND ANALYTIC CALCULATIONS
In the MSSM there are two complex Higgs doublets, H u and H d , which give rise to five physical Higgs bosons: two CP-even (h, H), one CP-odd (A) and a charged pair (H ± ).
Due to the µ term appearing in the superpotential, the MSSM suffers from the µ-problem.
Besides, in order to give a 125 GeV SM-like Higgs boson, large corrections to the Higgs mass from heavy stops is needed, which will lead to the little fine tuning problem . To overcome these difficulties, we can go beyond the MSSM. One alternative is the NMSSM, which introduces a singlet Higgs field. In the NMSSM the µ term does not appear in the superpotential. Instead, it is generated when the singlet Higgs field develops a vev. Also, the SM-like Higgs boson gets an extra tree-level mass from the mixing with the singlet field and thus the stops are not necessarily heavy to push up the Higgs mass, which alleviates the little fine-tuning problem [30][31][32]. In the NMSSM the singlet Higgs field mixes with the other two doublet scalars. Then the Higgs sector contains seven Higgs bosons, i.e., compared with the five Higgs bosons in the MSSM, the NMSSM contains one more CP-even and one more CP-odd Higgs bosons. In the following H 1,2 denote the real scalar components of H d,u in the MSSM and H 1,2,3 denote the real scalar components of H d,u,s in the NMSSM. tan β ≡ v u /v d is also used in our paper (here H d , H u and H s are the down-type, up-type and singlet Higgs fields, respectively). One can get the mass eigenstates from the CP-even states: where U 2 i1 +U 2 i2 = 1, V 2 i1 +V 2 i2 +V 2 i3 = 1 and the h i is aligned by mass. The singlet contribution is reflected by the rotation matrix elements V i3 via the formula h SM = V h SM 1 H 1 + V h SM 2 H 2 + V h SM 3 H 3 (a large V h SM 3 means that h SM has a considerable singlet component).
In our calculations, we follow the simplified ACOT prescription to deal with the b-quark mass [34][35][36]. By including the QCD and SUSY-QCD effects to the bottom Yukawa cou-plings, we can respectively obtain the effective h i bb couplings in the MSSM [37][38][39][40][41][42][43][44][45] and NMSSM [46]: NMSSM : where Here it should be noted that due to the contribution of the singlet field to the effective potential, an additional correction term The v d and v s are the VEVs of the Higgs fields H u and H d respectively. The auxiliary function I is defined as . The value of m DR b is related to the QCD-MS mass m MS b (which is usually taken as an input parameter [47]) by where n is the number of active quark flavors and m M S b (µ R ) is taken as When Q 2 > Q 1 , the evolution factor U n reads where Since the ∆ b -related corrections have already been included into the tree-level contribution, we need the following counter terms to subtract them to avoid double counting in the one-loop calculations [37] MSSM : NMSSM : For SUSY-QCD corrections to bb → hh, the sbottoms and gluino are involved in the loops. The sbottom mass matrix is given by [48] where After diagonalizing Eq.(13), we can obtain the sbottom masses mb 1,2 and the mixing angle θb: The Feynman diagrams for one-loop SUSY-QCD corrections to bb → hh has been represented in [29]. To preserve supersymmetry, we adopt the dimension reduction method to regulate the UV divergences in the gluino and squark loops. Then we use the on-shell renormalization scheme to remove these UV divergences.

A. A scan of parameter space
We use NMSSMTools [49] and LoopTools [50] to perform a random scan over the parameter space and loop calculations. For simplicity, we assume an universal parameter M L3 for the slepton sector and fix all irrelevant soft parameters for first two generation of the squark sector to be 1 TeV. We also set M D3 = M U 3 and A b = A t for the third generation of the squarks. Besides, we impose the grand unification relation of the gaugino masses, In our scan we consider the following experimental constraints: (ii) The constraints from the precision electroweak data [53] and flavor physics at 2σ level; (iii) The dark matter relic density from Plank at 3σ level and the limit of direct detection from XENON100 [54]; (iv) The explanation of muon g − 2 at 2σ level [55].  In our scan, for each experimental data which has a central value, we require the samples to agree with the experimental data at 2σ level, except for the dark matter relic density which is required to agree with the measured value at 3σ level (we made such a choice just in order to be consistent with the analysis in the literature). For the LEP and Tevatron direct search bounds on sparticle masses, we just require the samples to satisfy such bounds.
For the LHC Higgs search of H/A → τ τ and H ± → τ ν τ , we require the samples to satisfy the upper limits. The scan ranges of the parameters are large, we keep the samples survived various experimental constraints as stated above. Besides, we further require gluino mass larger than 1 TeV to avoid multi-jets search on SUSY [56]. However, we did not impose other LHC direct limits on sparticles for the following reasons. First, we required the first and second generations of squarks to be 1 TeV and the gluino beyond 1 TeV. But the latest LHC search results gave more stringent constraints on such squark and gluino mass (the most stringent bound is for the CMSSM, which is mg > 1.7 TeV in case of mg ≃ mq and mg > 1.1 TeV in case of mq ≫ mg). Actually, our results are not sensitive to these masses. Second, the current LHC limit is about 500-600 GeV for stop and 400-600 GeV for sbottom [57].
However, such limits were obtained in some simplified model or by assuming a certain decay branching ratio to be 100%. In our case the stop and sbottom decays are quite complicated, which will weaken the LHC limits. Further, for electroweak gauginos and sleptons, the current LHC limits will also be weakened in our case for the same reason. After that we also require surviving samples to avoid Landau singularity at GUT scale and we checked that all of our surviving samples satisfy √ λ 2 + κ 2 < 0.75 in NMSSM. We note that a large tan β exist in the surviving samples of the MSSM, this is because that a 125 GeV neutral Higgs mass is guaranteed by a large A t (which provides X t /M s close to √ 6) even for tan β as large as 40. As for the flavor constraints, we projected our samples onto the tan β versus the charged Higgs mass plane and found that when tan β increases the charged Higgs mass grows dramatically (especially, for tan β close to 40, the charged Higgs mass is heavier than 700 GeV) and thus can satisfy the flavor constraints. For the samples surviving the above constraints (i)-(iv), we further perform a fit by using the available Higgs data at the LHC.
We define the Higgs signal strength µ i as where p is the Higgs boson production mode and i stands for the measured channels by Tevatron, ATLAS and CMS collaborations. For each production mode p, its contribution to the channel i can be determined by the selection efficiency ǫ i p [58]. We summarize all experimental signal strength µ exp i with their 1σ error-bars and selection efficiencies in Fig.1.
We can see that most measurement results are consistent with the SM predictions. The We use the combined Higgs mass M exp h = 125.66 ± 0.34 GeV [60]. The χ 2 definition in our fit is where σ i and σ M h only denote the experimental errors.

B. The cross section of bb → hh with SUSY-QCD correction
We use CTEQ6L1 and CTEQ6m [61] for the leading order and SUSY-QCD calculation, respectively. The renormalization scale µ R and factorization scale µ F basically can vary between M h /2 and 2M h . In order to compare our results with [29] where µ R = µ F = M h /2 is assumed, we also made this assumption in our calculation. The input parameters of the SM are taken as [62] m In Fig. 2, we display the parameter space satisfying the experimental constraints (i-iv), showing the cross sections of the SM-like Higgs pair productions via bb annihilation (with SUSY QCD correction) and gg fusion versus M A at the 14 TeV LHC in MSSM and NMSSM.
In this paper we aim to investigate the property of the bb → hh production by including the SUSY QCD corrections. For the gg → hh production, we only calculate its cross section at one-loop level, not including the SUSY QCD corrections due to its small relative correction [63] comparing the SUSY QCD correction on bb → hh process. We used our own codes and combined them with Looptools to do our calculation. We checked our results with [27] and found good agreement.
We checked that our results agree with [29] for bb → hh and with [26] for the gluon fusion process.We can see that due to the constraints from the LHC and B-physics, such as , the values of m A must be larger than about 300 GeV. In the MSSM the maximal cross section can still reach 50 fb at 14 TeV LHC, which can be competitive with gg → hh. However, we also notice that the hadronic cross section proceeding through bb → hh deceases when m A or tan β becomes large. The reason can be understood as follows. On the one hand, for a moderate m A , the dominant contribution to bb → hh comes from the resonant production bb → H → hh. With the increase of M A , the mass of H gets heavy and then the production rate of bb → hh is suppressed. Besides, the coupling of hhH will approach to zero for a large m A and also leads to the reduction of the cross section. On the other hand, for a small tan β, H has a large branching ratio into a pair of Higgses hh [65], for a large tan β, the production rate of bb → H can be enhanced but the branch ratio of H → hh is highly suppressed. So the total production rate of bb → hh will become small.The decoupling behavior of the cross section proceeding through gg → hh can be understood with the following considerations: To predict a 125 Gev Higgs boson, a large A t is required, which induces a sizable SUSY effect for the process gg → hh. M A affects the process gg → hh mainly through the Higgs mass m h . So when we require m h in the range of 123-127 GeV, the process gg → hh is not sensitive to M A . Further, since gg → hh is dominated by the stop loops, the value of tan β affects this process through the coupling ht itj . Because this coupling is not sensitive to tan β for our surviving points, our results depend weakly on tan β.
In NMSSM the SM-like Higgs boson h with mass around 125 GeV can be either h 1 or h 2 .
However, we focus on the h = h 2 scenario that is more welcomed by the naturalness. From Fig.2 we can see that the maximal cross section of bb → hh can only reach about 4 fb, which is much smaller than gg → hh. We find that the suppression of bb → hh in NMSSM mainly has two reasons. One is that in NMSSM the tan β value is around 3-5 which is much smaller than in MSSM which is always larger than 10. To further investigate the influence of the Higgs data in Fig.2 on the SUSY-QCD effect in bb → hh, we define the relative SUSY-QCD correction δ SQCD as In our calculation we use the α LO s for the LO cross-section and α N LO s for the NLO crosssections, respectively. In Fig. 3 we show the dependence of δ SQCD for the bb → hh on the SUSY parameters M A , tan β, the lightest sbottom mass (mb 1 ) and gluino mass (mg) in the MSSM. In this figure the samples satisfying the experimental constraints (i-iv) are further classified according to the Higgs data: We use the χ 2 and the degree of freedom to calculate the p-value for each point and plot the points whose p-values are larger than 0.045 (2σ) and From the lower panel we note that for heavy mb 1 and mg, the SUSY-QCD effects decouple slowly. This behavior is because that the SUSY-QCD corrections depend on the ratio of the SUSY parameters. For example, in the triangle diagrams, the SUSY-QCD correction to the vertex hbb is proportional to M 2 EW /M 2 A and M 2 EW /M 2 b [42,66]. So only when all the sparticles and m A are heavy, the SUSY-QCD effect can completely decouple from the process of bb → hh.
The relative SUSY-QCD corrections for the bb → hh in the NMSSM are presented in Fig.4. It can be seen that the maximal SUSY-QCD correction can reach 15% for the samples in 1σ range. From the upper panel we can see that δ SQCD becomes small with the increase of λ or m h 3 . The reason is that with the increase of the λ, the m h 3 gets heavy and its contribution to the cross section becomes small. From the lower panel we see that, due to the residual effects of the sparticles, the SUSY-QCD corrections can still reach about 9% for heavy sbottom and gluino.
In Fig.5 we show the total cross section of the Higgs pair production at the 14 TeV LHC (via both bb annihilation and gg fusion) for the samples in the 1σ and 2σ ranges of the Higgs data. We can see that in the 1σ range the total cross section can be maximally enhanced by a factor of 2.7 and 2.2 in the MSSM and NMSSM, respectively.
Finally, considering the null results of the direct search for sparticles at the LHC, we investigate the SUSY-QCD effect in Higgs pair production in the limit of heavy sparticles.
For simplicity, we assume a common mass M SU SY for all relevant SUSY mass parameters: Fig.6  We can see that for M SU SY = 1 TeV, the ratios will maximally reach 3 and 2 in the MSSM and NMSSM, respectively. When M SU SY goes up to 5 TeV, the enhancements become weak but can still reach 1.8 and 1.4 in the MSSM and NMSSM, respectively. So the effects of heavy sparticles decouple quite slowly from the Higgs pair production. We checked that the SUSY effects decouple quickly in bb → hh but slowly in gg → hh.

IV. CONCLUSION
We considered the current experimental constraints on the parameter space of the MSSM and NMSSM. Then in the allowed parameter space we examined bb → hh (h is the 125 GeV SM-like Higg boson) with one-loop SUSY QCD correction and compared it with gg → hh.
We obtained the following observations: (i) For the MSSM the production rate of bb → hh (with one-loop SUSY QCD correction) can reach 50 fb and thus can be competitive with gg → hh, while for the NMSSM bb → hh has a much smaller rate than gg → hh due to the suppression of the hbb coupling ; (ii) The SUSY-QCD correction to bb → hh is sizable, which can reach 45% for the MSSM and 15% for the NMSSM within the 1σ region of the Higgs data; (iii) In the heavy SUSY limit (all soft mass parameters become heavy), the SUSY effects decouple rather slowly from the Higgs pair production, which, for M SUSY = 5 TeV and m A < 1 TeV, can enhance the production rate by a factor of 1.5 and 1.3 for the MSSM and NMSSM, respectively. Therefore, the Higgs pair production may be helpful for unraveling the effects of heavy SUSY.