1 Introduction

In the next-to-minimal supersymmetric (NMS) Standard Model (SM), it is possible for either one of the additional singlet-like scalar and pseudoscalar Higgs bosons to be almost degenerate in mass with the  GeV SM-like Higgs state. In the real NMSSM (rNMSSM), when the mass difference between two scalar states is comparable to their individual total decay widths, the quantum mechanical interference, due to the relevant diagonal as well as off-diagonal terms in the propagator matrix, between them can become sizeable. This possibility invalidates usage of the narrow width approximation (NWA) to compute the cross section for the production of a di-photon pair with a given invariant mass via resonant Higgs boson(s) in the gluon fusion process at the Large Hadron Collider (LHC). When, motivated by the baryon asymmetry of the universe, CP-violating (CPV) phases are explicitly invoked in the Higgs sector of the NMSSM, all the interaction eigenstates mix to give five CP-indefinite physical Higgs bosons. In this scenario, the interference effects due to the off-diagonal terms in the Higgs mass matrix that mix the pseudoscalar-like state with the SM-like one can also become significant, when these two are sufficiently mass-degenerate. We perform a detailed analysis, in both the real and complex NMSSM, of these interference effects, when the full propagator matrix is taken into account, in the production of a photon pair with an invariant mass near 125 GeV through gluon fusion. We find that these effects can account for up to % of the total cross section for certain model parameter configurations. We also investigate how such mutually interfering states contributing to the  GeV signal observed at the LHC can be distinguished from a single resonance.


Two Higgs bosons near 125 GeV in the NMSSM: beyond the narrow width approximation

Biswaranjan Das, Stefano Moretti, Shoaib Munir and Poulose Poulose

[0.25cm] Department of Physics, IIT Guwahati, Guwahati, Assam 781039, India

[0.25cm] School of Physics and Astronomy, University of Southampton,

Highfield, Southampton SO17 1BJ, UK

[0.25cm] School of Physics, Korea Institute for Advanced Study, Seoul 130-722, Republic of Korea

1 Introduction

The discovery of a Higgs boson [1, 2] at the LHC provides convincing evidence of spontaneous electro-weak (EW) symmetry breaking (SB) through the Higgs mechanism. It is also intriguing that the subsequent measurements of its properties have shown their remarkable agreement with the expectations from the SM of particle physics. These measurements include the signal rates and coupling strengths in the various Higgs boson production and decay channels that have so far been analysed, as well as its spin and parity. However, the shortcomings of the standard Higgs mechanism, including primarily the stability of the mass of the Higgs boson against large quantum corrections, need to be addressed properly in order to completely understand the dynamics of EWSB. Putting this together with other unresolved issues in the SM, such as its inability to explain the mass of neutrinos, the nature of dark matter (DM) and the large baryon asymmetry of the universe, compel us to believe that the elementary particle spectrum could be richer than the minimal one embedded in the SM.

Supersymmetry (SUSY), proposed originally as a remedy for some of the above problems faced by the SM, presents an appealing explanation for the stability of the Higgs mass, while providing also a natural candidate for DM. However, its minimal manifestation, known as the minimal supersymmetric Standard Model (MSSM), is becoming increasingly constrained by the LHC measurements of the Higgs boson properties, besides becoming more and more fine-tuned in explaining null searches for its own direct signatures (i.e., of sparticle states). The MSSM is also troubled by issues arising purely from naturalness considerations, like the presence of a quadratic Lagrangian term with a new mass parameter, , which is phenomenologically required to lie at the EW scale but has no theoretical ground to do so [3, 4]. The NMSSM was proposed to take care of this so-called -problem through the introduction of a new singlet scalar superfield [5]. The presence of this superfield leads to two additional neutral Higgs mass eigenstates in the NMSSM (see, e.g., [9, 10] for reviews) compared to the MSSM. When all the parameters in the Higgs and sfermion sectors are assumed to be real and hence CP-conserving (CPC), one of these two additional states is a scalar and the other a pseudoscalar. With five neutral Higgs bosons in total, the real NMSSM (rNMSSM) provides some unique possibilities for collider phenomenology.

In multi-Higgs models like the MSSM and NMSSM, the observed signature near 125 GeV can in principle be explained in terms of a single Higgs resonance or, alternatively, two or more Higgs resonances that cannot be individually resolved by the experiment as yet. In the MSSM, the possibility of the two scalars both lying near 125 GeV is ruled out by the fact that such a mass for one of the scalars generally necessitates the other scalar to be essentially decoupled, and hence much heavier or lighter (see, e.g., [11]). In the rNMSSM, in contrast, it is still a plausible scenario, as discussed in [12, 15]. The alternative possibility of a pair of  GeV scalar and pseudoscalar has also been studied in [17].

In the NMSSM, to address the baryonic asymmetry of the universe, CP-violation can be invoked directly and explicitly in the Higgs sector at the tree level, unlike in the MSSM, where it is only radiatively induced into the Higgs sector beyond the Born approximation. This is done by taking the Higgs sector trilinear couplings to be complex parameters, hence we refer to this version of the model as the complex NMSSM (cNMSSM) here. For non-zero CPV phases of these parameters, instead of the distinct CP-even and CP-odd Higgs bosons, the model contains five CP-indefinite neutral states. This CP mixing in the Higgs sector can get additional contributions from the complex Higgs-sfermion-sfermion couplings as well as the phases in the neutralino-chargino sector, as in the CPV MSSM. Consequently, depending on the sizes of these phases, the mass spectrum and production/decay rates of the Higgs states can get considerably modified compared to the CPC case [18], similarly to the MSSM [19, 26, 27]. The phenomenology of a single CPV Higgs boson near 125 GeV in the cNMSSM has been studied for a variety of possible underlying scenarios in [18, 30, 31]. The specific case of two mass-degenerate Higgs resonances near 125 GeV within the cNMSSM has also been considered in [32], where it was shown that these can give a better fit to the LHC Run-I data, performed using the program HiggsSignals [33], compared to a single resonance.

In all the above-mentioned analyses of two  GeV Higgs bosons in the NMSSM, it was assumed that each of them is produced on-shell and decays subsequently into any of the observed final states. However, for very strong mass degeneracy, it is possible that the two Higgs states produced in gluon fusion oscillate into each other before they decay. This could be perceived as being facilitated through quantum corrections of the propagator with two different mass eigenstates at the two ends. Such effects, coming into play beyond a Breit-Wigner (BW) resonance, in general contexts as well as in specific scenarios like the CPC MSSM, have been considered in some studies [34]. The specific case of the CPV MSSM with one-loop effects in the full propagator was treated in [26].

The purpose of the present work is to explore this possibility of quantum mechanical interference between two Higgs states near 125 GeV in the NMSSM, both real and complex. We shall demonstrate here that the inclusion of such effects enhances the span of the model solutions mimicking the LHC observation. We shall then investigate possible ways to identify signatures of such a coupled system of Higgs bosons through a shape study of the profile of the resonance in the invariant mass distribution of the di-photon decay products. The analysis is carried out by first performing numerical scans of the model parameter space to identify specific Benchmark Points (BPs) where two of the Higgs bosons are degenerate with mass around 125 GeV, within the uncertainty allowed by present LHC measurements. For these BPs, we then calculate the cross section for the production of a di-photon pair with invariant mass near 125 GeV via Higgs resonance(s) in the gluon fusion process at the LHC, using a Monte Carlo (MC) integration code developed in-house. This cross section is calculated assuming three different approaches: the full Higgs propagator matrix in the amplitude; only diagonal terms in the propagator matrix; for the two Higgs bosons individually without any mutual interference. A comparison of these three cross sections shows significant effects of interference, with the full propagator case deviating by up to about 38% compared to the sum of the two individual Higgs boson contributions, along with a hint of a distortion in the line-shape of the differential distribution.

The article is organised as follows. In the next section we will briefly revisit the Higgs sector of the NMSSM. In section 3 we will derive the analytical expression for the cross section that includes the full Higgs propagator matrix. In section 4 we will provide details of our methodology for the numerical analysis. In section 5 we will discuss the results of our analysis and in section 6 we will present our conclusions.

2 The NMSSM Higgs sector

The NMSSM is defined by the superpotential


where , and are the quark and lepton Yukawa coupling constants, and are the left-handed quark and lepton doublet superfields, , and are the right-handed up-type, down-type and electron-type singlet superfields, respectively, and the charge conjugation is denoted by the superscript . Furthermore, and in the above superpotential are Higgs doublet superfields with opposite hypercharge, , as in the MSSM, and is a Higgs singlet superfield. Here, and are dimensionless trilinear coupling constants.

The fourth term on the right hand side of Eq. (1) replaces the Higgs-higgsino mass term, , present in the MSSM superpotential. The Higgs singlet superfield acquires a non-zero vacuum expectation value (VEV), , after EWSB. This is naturally of the order of the SUSY-breaking scale (herein operationally defined as , where and are the physical masses of the two stops), thus solving the -problem of the MSSM by generating an effective -term,


The absence of a term, however, results in a symmetry, which is explicitly broken here by the last term in Eq. (1), thus introducing instead a discrete symmetry and making the NMSSM superpotential scale-invariant as well.

The Higgs potential, derived from the above superpotential, is given as


where and are the and gauge coupling constants, respectively, and . Here, and are soft SUSY-breaking Higgs trilinear couplings, while , and denote the soft Higgs masses. The fields , and are expanded about their respective VEVs, , and , as


For correct EWSB, the , rewritten in terms of these expanded fields, should have a minimum at non-vanishing , and , implying


which leads to six ‘tadpole conditions’ (see, e.g., [31]).

Taking the second derivative of at the vacuum yields the tree-level neutral Higgs mass matrix-squared, , in the basis . It can be expressed in the general form


where the block corresponds to the CP-even interaction states (, , ), the block to the CP-odd states (, , ), while is responsible for mixing between the CP-even and -odd states.

In the rNMSSM, where all the Higgs sector trilinear coupling parameters are real, is a null matrix. One can therefore simply rotate only the submatrix to isolate the massless Nambu-Goldstone boson field, , which can then be dropped to yield a mass matrix . This mass matrix receives higher order corrections, , from various sectors of the model, and thus gets modified as


The most dominant of these corrections can be found in [9, 38]. Some further two-loop corrections have been calculated in [39] and [40]. After the inclusion of these corrections, the submatrices and are separately diagonalised to obtain the three CP-even mass eigenstates, (with ), and the two CP-odd physical Higgs bosons, (with ).

On the other hand, one can also invoke CP-violation explicitly by assuming , , and . This leads to non-zero entries in the submatrix, implying that the CP-even and CP-odd interaction eigenstates also mix mutually. In this cNMSSM, the state is first separated out through a rotation of the entire by ,


and dropped before calculating the higher order corrections to the resulting . The complete expressions for the tree-level and the dominant one-loop contributions to from the (s)quark and gauge sectors in the cNMSSM were studied in [41, 43, 44] in the renormalisation group equations-improved effective potential approach. The corrections from the gaugino sector were added in [31] and, more inclusively, in [45]. In the Feynman diagrammatic approach, the complete one-loop Higgs mass matrix was derived in [30] and the contributions to it were calculated recently in [46].

The physical Higgs mass eigenstates of the cNMSSM are then obtained from the interaction states through another rotation by ,


resulting in the diagonalised squared mass matrix,


Here , , , and are the five neutral CP-indefinite Higgs bosons, ordered in terms of increasing mass.

Note here that only the physical phase combination appears in the cNMSSM Higgs sector at the tree-level, as the other possible phase combinations involving and are determined by it, up to a twofold ambiguity, through the tadpole conditions. Beyond this level, the CPV phases of the gaugino mass parameters, , and of the soft trilinear couplings, , of the Higgs boson to the sfermions also get radiatively induced into this sector.

3 Di-photon production via gluon fusion

We now turn our attention to the process under scrutiny, i.e., di-photon production from gluon fusion via Higgs states at the LHC. The squared amplitude for the process, with collectively denoting the five neutral CP-indefinite Higgs bosons, can be written as [47]


where are the gluon and photon helicities, respectively, and is the Higgs boson propagator, with being the squared center-of-mass (CM) energy of the incoming gluons. The amplitude for the production part is given as


where the scalar and pseudoscalar form-factors are [48]


with and the loop function being


Note that and in Eq. (13) are the couplings of to quarks and squarks , respectively, which depend on the elements of the Higgs mixing matrix, , noted in Eq. (10) above. The exact forms of these couplings can be found in [44].

The amplitude for the decay part is similarly given as


with the form-factors being


where is the colour factor for (s)quarks and charged (s)leptons, respectively, with being the corresponding electric charge. Finally, the full propagator in Eq. (11) is a matrix111Assuming negligible off-resonance self-energy transitions of any of the five Higgs bosons to the would-be Goldstone boson, ., given as


where , with being the absorptive parts of the Higgs self-energies, for . The diagonal absorptive parts are equivalent to the widths of the corresponding Higgs states, . The explicit expressions for are given in the Appendix. Note here that in our numerical analysis we will only focus on and having masses very close to 125 GeV. This implies that essentially the propagator matrix elements with are the only ones that contribute to the production of the di-photon pair with invariant mass near 125 GeV, assuming that the remaining three Higgs bosons are relatively heavy. Moreover, in the case of the rNMSSM, since the CP-even-odd mixing terms in the Higgs mass matrix vanish, it is sufficient for our purpose to consider only the propagator matrix corresponding to the CP-even states instead of the complete one above.

When the splitting between the Higgs boson masses is much larger than the sizes of the absorptive parts in Eq. (19), the NWA can be applied to the th Higgs boson propagator as


so that the partonic cross section becomes


The total cross-section for the process is then written as


where and are the parton distribution functions (PDFs) of the two gluons. By substituting in the above equation as


where is the total CM energy of the system, and performing the integration over , one gets


In contrast, when two (or more) Higgs bosons of the model are almost degenerate in mass near a given , the sizes of the corresponding absorptive parts can become comparable to their mass difference. As a result, the th Higgs state can undergo resonant transition to the th state through quantum corrections, as shown in Fig. 1. In such a scenario, the NWA is no longer valid, and the total cross section is given as


where is the propagator matrix given in Eq. (19). From the above equation, one obtains the differential cross section (recall that ) as

Figure 1: Illustration of the effect of mixing in the propagator induced by quantum corrections.

4 Numerical analysis

We first performed numerical scanning of the parameter space of the NMSSM, requiring and to lie within the range222The extended range of Higgs boson masses around the actual measured experimental value of  GeV is to allow for up to  GeV uncertainty from unknown higher order corrections in their model prediction.. Our first scan corresponded to the rNMSSM, wherein sufficient mass degeneracy near 125 GeV between the two lightest scalars can generally be obtained for large values of the couplings and and a relatively small , which results in maximal mixing between the doublet- and singlet-like states, as noted in some earlier studies [12]. In the rNMSSM, while it is also possible for to lie near 125 GeV [17], it does not mix with the SM-like when the coupling parameters are all real. Therefore, the corresponding off-diagonal absorptive parts in the propagator matrix given in Eq. (19) are zero. When the complex phases are turned on though, all the Higgs states become CP-indefinite, and any of the off-diagonal terms in the full propagator matrix can be non-zero and contribute to the interference effects. Therefore, either one of the scalar-singlet-like and pseudoscalar-singlet-like states can have strong mass-degeneracy with the  GeV SM-like state and interfere with it.

As stated earlier, at the tree level, only the phase combination appears in the Higgs sector of the cNMSSM. Furthermore, several studies [30, 43, 44] have shown that, out of all the individual phases, including those appearing beyond the Born approximation, the phase is virtually unconstrained by the measurements of the fermion Electric Dipole Moments (EDMs). Therefore, after setting all the other phases to , we performed two separate parameter space scans of the cNMSSM also, with the value of fixed to in one and to in the other. In Tab. 1 we list the scanned ranges of the free parameters (input at the EW scale), which assume the following universality conditions:

where and are the soft masses of the sfermions, those of the gauginos and the soft trilinear couplings. These ranges are consistent across the three scans and correspond to the parameter space region that was noted to yield maximally mass-degenerate and in a previous study [32], where more details about the scanning methodology can also be found. It was additionally pointed out in that study that for larger values of it gets increasingly difficult to obtain both and near 125 GeV in the cNMSSM.

Parameter Scanned range
 (GeV) 800 – 2000
 (GeV) 100 – 500
  (GeV) – 0
2 – 8
0.58 – 0.7
0.3 – 0.6
 (GeV) 100 – 200
 (GeV) 200 – 1000
 (GeV) – 0
Table 1: NMSSM parameters and their scanned ranges.

For each parameter space input point generated by the scanning algorithm, the masses as well as branching ratios (BRs) of the Higgs bosons were calculated with the public code NMSSMCALC v2.00 [48]. The Supersymmetric Les Houches Accord [49] output file produced by NMSSMCALC for a scanned point was then passed to HiggsBounds v4.3.1 [51] to check for the consistency of each Higgs boson with the direct Higgs search results from LEP, Tevatron and LHC. We further made sure that a point only got through the scan if it satisfied the limits from measurements of the EDMs, computed intrinsically by NMSSMCALC. Finally, the CMS and ATLAS collaborations have performed measurements of the total width of the  by analysing its off-shell production and subsequent decays in the and channels [55, 56]. The most recent observed 95% confidence level upper limit for the two channels combined is 13 MeV. Therefore, we also require each of the and in a given scan to observe this constraint, unless stated otherwise for exceptional scenarios, which may well be plausible, as such a limit presumes an underlying BW resonance for the signal [58, 59], which is not the case here.

Next, from the points collected in each scan, we selected BPs satisfying certain specific criteria, which will be explained later. In order to perform the numerical calculation of the cross sections for these BPs, we implemented the expressions given in Eqs. (25) and (26) in a locally developed fortran program. This code is linked to the LAPACK package v3.6.0 [60] for propagator matrix inversion, as well as to a locally modified version of the VEGAS routine [61] to perform the 2-dimensional numerical integration. As a test of the reliability of our results, for a given model point, we calculated the cross section in the NWA for each of the two Higgs bosons with our code and compared it with the gluon fusion cross section computed using the publicly available code SusHi v1.6.0 [62] multiplied by its di-photon BRs obtained from NMSSMCALC. We found that the two results agreed within 5% or better in all cases. The various Higgs boson couplings for a given parameter space point, needed for the calculation of the absorptive parts of the propagator matrix as well as of the production and decay form-factors in our code, were also obtained from NMSSMCALC.

Note that our program calculates the total cross section only at the Leading Order (LO), since implementing Higher Order (HO) corrections, e.g., as included in SusHi, is a highly involved task beyond the scope of this work, which is aimed at comparing the effects of including the full propagator against the simplest approach of two separate BWs on the total cross section. In fact, since these HO corrections apply only to the production process, they should have exactly the same impact in both approaches, hence including their effect is simply tantamount to rescaling the cross section by a ‘ factor’ (defined as , with HO implying the perturbative order at which the cross section is to be evaluated), which can be obtained from a dedicated public tool. For a few test points corresponding to the rNMSSM, using SusHi, we thus also estimated the Next-to-Next-to-LO (NNLO) factor, , by calculating the gluon-fusion cross section at both LO and NNLO in QCD. We found this to always lie very close to 3 for both and . However, since SusHi is not yet compatible with the cNMSSM, the factor cannot be estimated exactly for the points corresponding to the cNMSSM. Therefore, for an approximate evaluation of the NNLO cross section in our discussion below, we will multiply the LO cross section obtained with our tool for a given point, both in the real and the complex NMSSM, with a constant . (In fact, for an almost decoupled SUSY sector, as is generally the case here, the factor is indeed a constant for a SM-like  GeV Higgs boson, which can be obtained from, e.g., [65].)

Finally, we used the CT10 [66] PDF set for gluons in our cross section calculation, with renormalisation/factorisation scale set to the  mass, i.e., 125 GeV. We have, however, verified that the gross features of our analysis are independent of the choice of the PDF set as well as of the fixed numerical inputs for the SM and NMSSM parameters, for which the default NMSSMCALC values were retained.

5 Results and discussion

In Fig. 2 we show those points obtained from the three scans for which is smaller than one (or both) of the widths, and , of the two lightest Higgs bosons. The top panel corresponds to the rNMSSM, while the bottom panels to the cNMSSM with (left) and with (right). We note in the figure that, given the parameter space in Tab. 1, for the vast majority of points, and tend to lie within 3–4 MeV of each other in the rNMSSM. For , the size of splitting between and can range from very small to fairly large across the points collected. However, for , no points appear along the diagonal for  MeV in the bottom right frame, as the splitting between these two widths starts growing beyond this value.

From each of these scans we selected a few BPs to study the cross section for the production of a di-photon pair with an invariant mass, , near 125 GeV via resonant Higgs boson(s) in gluon fusion at the LHC with  TeV. More specifically, we studied the distributions, with respect to the partonic CM energy (which is the same as at LO), of the differential cross section, calculated such that:
Case 1: the and terms in the propagator matrix, Eq. (19), each contribute alternatively to two amplitudes, which are squared and then summed, implying two independent Higgs boson resonances;
Case 2: both and contribute to the amplitude which is then squared, thus allowing for interference between and but without any mixing effects;
Case 3: besides and , the off-diagonal terms and also contribute to the amplitude before squaring, leading to additional interference effects arising from the mixing of the two Higgs states.

Figure 2: Points obtained from the parameter scans of the rNMSSM (top) and of the cNMSSM with (bottom left) and with (bottom right). For all the points shown, (colour map) is always smaller than (x-axis) and/or (y-axis).

The distributions obtained for the Cases 1, 2 and 3 (colour-coded in red, green and blue, respectively) are plotted in Fig. 3 for each of the three selected BPs corresponding to the rNMSSM, with the integrated cross section for each curve also given in the legend. For these distributions, a bin width of 2 MeV has been used. As noted above, for most of the points in the rNMSSM. Therefore, in order to illustrate the dependence of the interference effects on the mass difference and relative widths of the two Higgs bosons, we selected BP1 such that , BP2 such that and BP3 with . One sees, going from the top panel to the bottom right one in the figure, that these effects are always positive and grow larger as decreases compared to , as expected. Also, the interference effects due to the mixing terms in the propagators matrix (Case 3) are notably larger than those due only to the diagonal terms (Case 2) for each of the three BPs. The deviation in the total cross section with the full propagator matrix compared to the Case 1 for BP3 at an inclusive level is about 38%, clearly indicating that the interference effects can be quite sizeable. We point out here that although the BP3 represents maximal enhancement of these effects among all the points collected in our scans, it is possible that they can be even slightly larger for certain other parameter combinations in the vicinity of this BP. Note also that the total integrated cross section (obtained at NNLO in QCD, as mentioned earlier) is generally consistent with the fiducial one, as estimated for the SM-like Higgs boson near 125 GeV [67], or measured by the ATLAS and CMS collaborations for at  TeV.333We also point out here that a considerable discrepancy exists between the ATLAS measurement, which reads  fb [68], and the CMS one,  fb [69], besides a fairly large error in each of these itself. This renders an accurate estimation of the total cross section of little significance here and justifies our use of an approximate constant NNLO factor.

The values of the input parameters for all the selected BPs can be found in Tab. 2, and the masses and widths of and as well as the total cross sections corresponding to the three Cases for each of the BPs are given in Tab. 3.

Figure 3: Distribution of the differential cross section as a function of the di-photon invariant mass (assumed equal to ) for the three benchmark points in the rNMSSM. The red, green and blue curves correspond to the Cases 1, 2 and 3, respectively, discussed in the text.
BP tan
1 1380.9 458.51 4.39 0.6970 0.4594 423.23 113.60
2 1598.3 471.51 4.34 0.6907 0.4823 402.53 110.86
3 1498.2 379.87 3.91 0.6969 0.4538 385.05 117.92
4 1366.6 426.35 3.92 0.6878 0.4657 361.11 112.79
5 1476.6 363.81 4.67 0.6725 0.4304 485.87 120.41
6 1427.1 249.93 4.53 0.6852 0.3360 610.69 147.10
7 1350.2 23.24 4.50 0.6630 0.3053 618.04 148.83
8 1270.6 176.67 3.96 0.6781 0.4501 538.70 168.65
9 1491.9 167.11 5.22 0.6920 0.4599 797.56 175.84
10 1378.0 173.35 3.99 0.6877 0.4483 564.66 172.87
11 1416.6 170.40 4.45 0.6684 0.3853 687.11 177.72
12 1429.0 168.46 4.71 0.6562 0.4303 689.40 173.02
Table 2: Values of the input parameters for all the selected BPs. All dimensionful parameters are in units of GeV.

The enhancement in the interference effects for a larger difference between and is further confirmed by the distributions shown in the top panels of Fig. 4 for BP4 and BP5, in the cNMSSM with . For both these BPs, the interference is again constructive, as in the rNMSSM. It is, however, also possible for the interference to be destructive in this scenario. This is the case for BP6, shown in the bottom left panel of the figure, for which the cumulative cross section is smaller for Case 2 compared to Case 1. Turning on the mixing terms in the propagator matrix further contributes negatively to lower the cross section, although the overall effect is hardly per cent level for this particular parameter space point.

Figure 4: As in Fig. 3, for the BPs corresponding to the cNMSSM with .
BP (fb)
(GeV) (GeV) (MeV) (MeV) (MeV) Case 1 Case 2 Case 3
1 125.3688 125.3782 9.4 10.7 9.7 50.36 52.41 59.78
2 124.9498 124.9562 6.4 10.1 9.1 53.58 57.54 69.23
3 126.1641 126.1667 2.6 10.1 9.3 53.10 58.36 73.33
4 125.3960 125.4052 9.2 9.6 9.5 48.11 50.06 56.16
5 124.6742 124.6757 1.5 9.1 8.4 56.90 59.53 71.28
6 125.6285 125.6393 10.8 11.1 5.9 44.68 43.96 43.54
7 124.2607 124.2625 1.8 2.5 2.3 97.10 100.00 89.10
8 124.9873 124.9968 9.5 10.3 3.0 46.94 48.38 48.89
9 124.9669 124.9742 7.3 10.6 3.0 45.22 46.54 47.31
10 125.1874 125.1924 5.0 10.3 2.9 46.56 49.55 50.62
11 125.1826 125.1828 2.0 10.1 2.6 50.14 51.42 52.22
12 124.7542 124.7604 6.2 10.3 2.7 47.96 47.12 49.03
Table 3: The masses and total widths of and in the selected BPs. Also listed for each BP is the cross section for the process calculated in the three different ways considered.
Figure 5: As in Fig. 3, for the BPs corresponding to the cNMSSM with .

For BPs 4–6 above, the and are scalar-like, which is the case for almost all the points obtained in the scan for this scenario. We note here that the singlet-like Higgs boson near 125 GeV is classified as scalar (pseudoscalar)-like if its coupling to the gauge bosons are significantly larger (smaller) than those of the singlet-like , which itself also lies fairly close in mass.444Evidently, both and cannot be singlet-like, since in that case both of them will have highly reduced couplings to the SM particles and resultantly the di-photon production cross section will be extremely suppressed. Moreover, in such a scenario ought to have SM-like couplings and would therefore be ruled out by the collider limits tested against using HiggsBounds. While one of the main reasons for considering the cNMSSM was to explore the possibility of interference of a  GeV pseudoscalar-like Higgs boson with the SM-like one, only a handful of such points were found by our scan, wherein the widths of the and are always relatively small.

BP7 in Fig. 4 illustrates the scenario with a pseudoscalar-like , with its total width, as well as that of the SM-like , being smaller than 3 MeV. An intriguing feature of this point is that, while the overall interference effect is positive for Case 2, enhancing the cross section compared to Case 1, it becomes negative when the complete propagator matrix is included. The destructive interference here is much stronger than in BP6 above and has the desirable consequence of bringing the total cross section down to a level consistent with the LHC measurements noted earlier. The reason why the cross section for Case 1 for this point is considerably higher than that seen for the other BPs so far can be ascribed to the fact that the is much more SM-like here compared to the points where the scalar-like gets closer in mass to it owing to large singlet-doublet mixing.

In Fig. 5 we display the distributions for the five BPs selected in the cNMSSM with . As noted from the corresponding scatter plot above, this value of allows for a much larger splitting between and , compared particularly to the rNMSSM. This makes it possible to test the impact of interference when , instead of being smaller than both the widths, lies somewhere in between them. Thus, for each of these points  MeV and  MeV, with varying from about 9.5 MeV for BP8 to 2 MeV for BP11. While the behaviour of the interference is similar to what has been observed earlier, i.e., it grows larger as the gap between and increases, its overall size remains comparatively small, reaching about 9% for BP10, for which is only slightly larger than . However, when is lowered below also, as in BP11, the interference effects get reduced instead of continuing to grow. BP12 is another example of mutually opposite contributions to the interference effects from the diagonal and off-diagonal elements of the Higgs propagator matrix. As opposed to BP7 though, here the negative interference comes from the diagonal elements in the propagator matrix, while the mixing effects contribute positively to again raise the cross section slightly for Case 3.

A closer look at the curves for BP6 and BP10 above reveals very small kinks near in addition to tall peaks near . These kinks result from the large splitting between the and , coupled with the fact that , while still being sufficiently smaller than the to cause notable interference effects, is larger than the bin size of 2 MeV adopted for plotting the distribution. Thus, for these points not only does the values of the inclusive (i.e., integrated) cross section change between the respective Cases 1 and 3, but also the shape of the distribution for the exclusive (i.e., differential) one. However, a bin width of 2 MeV is about three orders of magnitude smaller than the current experimental resolution of around 1 GeV [70]. Evidently, any differences between their shapes corresponding to each of the three Cases for a given BP, which could prove crucial for mutually distinguishing them and hence unraveling the interference effects, would not appear had a realistic bin width of 1 GeV been used. It is nevertheless interesting to study whether these differences persist to some extent once the differential distributions are convolved with a Gaussian distribution emulating detector effects. We performed the convolution using the ListConvolve function [71] in Mathematica by supplying the differential distribution data for a point as a list, as well as specifying the mean and width of the Gaussian.

In the top frames of Fig. 6 we display the result of the convolution of the distributions corresponding to Cases 1 and 3 for BP10 with a Gaussian of width (resolution) 1 GeV, which is also used as the size of the bins for first reproducing these distributions with our MC integrator. We use two prospective integrated luminosities at the LHC: 300 fb (left), which is expected by the end of the machine run in standard configuration; and 1000 fb (right), foreseen for the High-Luminosity (HL) LHC option [72].555The higher luminosity only serves to reduce the sizes of the error bars in the figure, which refer only to the statistical error in fact, as we are not able to account for the systematic one. It is worth appreciating in the figures that, while the kinks have expectedly disappeared and the distributions are smoother, there exists some marginal scope at the LHC to distinguish at the differential level the simplistic scenario generally assumed (Case 1) and the one rigorously predicted (Case 3), as the difference in the heights of the red and blue curves is not constant for all the bins in . However, evidently this requires knowledge of additional model inputs (e.g., the Higgs boson couplings) as such difference could also be generated by a different point in the parameter space. This difference is even more pronounced in the bottom panels of the figure, which correspond to the convolution with a Gaussian of resolution 300 MeV (clearly an unrealistic value at present, yet potentially within the reach of a detector upgrade). In fact, the histograms referring to these two predictions could eventually be statistically separable, again though, with the aforementioned caveat that one ought to have established additional model parameters elsewhere.

Figure 6: Convolution of the distributions 1 and 3 for BP10 with Gaussians of width 1 GeV (top) and 300 MeV (bottom). An integrated luminosity of 300 fb is assumed in the left panels and of 1000 fb in the right panels.

The difficulty to separate the Cases 1 and 3 for BP10 (as it would be for all other BPs) with present machine and detector conditions at the LHC is ultimately related to the enforcement of the  MeV constraint from off-shell Higgs measurements in our selection of the BPs. We therefore consider a Test Point (TP) 1, where this constraint is dismissed and only the milder  MeV constraint, as obtained from a global fit to the on-shell Higgs boson signal strength measurements [73], is imposed. This is all the more important in light of the fact that some critiques have been drawn about the model-independence of such a measurement, see, e.g., [74] and [75], or its stability against theoretical uncertainties [76]666Recall also that the mentioned experimental measurement of the (off-shell) SM-like Higgs width suffers from a small signal yield and large backgrounds.. The top right frame of Fig. 7 shows the convoluted distributions 1 and 3 for TP1 with a Gaussian of 1 GeV width, illustrating again the fact that also in this case there exists some scope in separating the Cases 1 and 3, even for current di-photon mass resolutions, so long that sufficient luminosity is accrued. The bottom frames of the figure illustrate that this scope gets further enriched if the mass resolution is improved to 300 MeV.

Figure 7: Top: The differential distributions for TP1 without convolution (left) and after convolution with a Gaussian of width 1 GeV for an integrated luminosity of (right). Bottom: TP1 distributions after convolution with a Gaussian of width 300 MeV for an integrated luminosity of (left) and (right).

In fact, one could ignore constraints on the total width altogether in order to estimate the minimal mass splitting that could be potentially detectable. While this exercise may appear academic (i.e., to dismiss a crucial experimental constraint), it is worth noting that the current procedures adopted to extract the Higgs boson properties inevitably work with the underlying assumption that only one resonance is produced around 125 GeV. This implies that the allowed intrinsic widths of a pairs of degenerate Higgs states need not relate directly to the currently fitted value.

In Figs. 8 and 9 we thus display the distributions of Cases 1 and 3 for two more TPs, 2 and 3, respectively, convolved, again, with Gaussians of 1 GeV and 300 MeV widths for two prospective integrated luminosities. In both these TPs, the is very wide, (100) MeV, while is a few 10s of MeVs, as seen in Tab. 4.777The input parameters for the three TPs are provided in Tab. 5. But since is only about 70 MeV for TP2, while it is larger than half of for TP3, the interference effects are highly enhanced for the latter (about 46%) compared to the former (%). These figures more effectively bring home the point that a very large (as noticeable in the top-left frames) does not impact significantly the quality of the fit to what, in the end, looks like a single object shape (as visible in the other three frames). Though, clearly, the difference between the Cases 1 and 3 is much more pronounced here than for TP1 (and all the BPs). This difference may potentially be established experimentally within the next few years, more likely so the wider (one of) the Higgs states.

Figure 8: As in Fig. 7, for the TP2.
Figure 9: As in Fig. 7, for the TP3.
TP (fb)
(GeV) (GeV) (MeV) (MeV) (MeV) Case 1 Case 2 Case 3
1 124.7928 124.8158 23 10.8 38.3 4.61 4.76 4.94
2 123.8696 124.1991 329.5 400.2 73.5 0.353 0.385 0.458
3 123.4590 123.7876 328.6 704.9 39.2 1.09 1.45 1.58
Table 4: Higgs boson masses and widths as well as the cross sections corresponding to the three Cases for the three selected TPs.
TP tan
1 1438.0 255.53 4.80 0.6935 0.3287 653.21 145.77
2 1405.7 154.63 5.63 0.6861 0.4602 546.59 108.58
3 1895.2 115.14 1.76 0.6524 0.5752 74.865 105.95
Table 5: Values of the input parameters for the three TPs considered. All dimensionful parameters are in units of GeV.

6 Conclusions

In summary, we have scrutinised in detail the proton-proton to di-photon process, through which a 125 GeV resonance consistent with the SM Higgs boson has been discovered at the LHC. Indeed, this is the signature for which the Higgs mass resolution is highest amongst all those accessible at the CERN collider. Measurements of its cross section, at both the inclusive and exclusive level, however, do not exclude the possibility of non-SM explanations. Amongst these, particularly intriguing are those invoking two Higgs bosons produced via gluon fusion, with such a small mass difference that they cannot be resolved by the experimental apparatus. This scenario can emerge only in non-minimal realisations of SUSY, such as the NMSSM, wherein (unlike the MSSM) two coexisting Higgs bosons can contribute to the 125 GeV signal (in as well as other final states). In this case, an accurate treatment is required of the propagation of the two states, which not only goes beyond the NWA but also allows for full interference between these. Hence, we have studied the quantitative impact of interference between two Higgs states near 125 GeV, with and without mixing effects, relative to the simplistic approach where the two resonant objects are treated independently of each other. For a full treatment, including the possibility of complex couplings as well, we have considered both real and complex NMSSM.

Our analysis involved scanning of the parameter space of the model for finding possible solutions consistent not only with the LHC exclusion limits on the additional Higgs bosons but also with the constraints from EDM measurements. These scans further collected only model solutions yielding two Higgs bosons with masses lying within the uncertainty of the measurements of the 125 GeV resonance. This was followed by a dedicated computation, performed with the help of a locally developed MC program, producing both integrated and differential cross sections for the full process . We have found that the aforementioned interference effects can be sizeable, with some of the selected BPs providing a difference of around 40% in inclusive rates between the standard approach consisting in treating the two resonances as separate BW functions and the full propagator one including all non-trivial quantum effects.

We then considered the possibility of a shape analysis of the emerging profile, which could reveal the presence of multiple resonances, assuming realistic, current and prospective, di-photon mass resolutions of the LHC detectors. This revealed some long-term potential to see the difference between the generally exploited simplistic case of assuming two separate resonances and the one where the two nearly mass-degenerate states interfere due to the inclusion of the complete propagator matrix in the amplitude calculation. These differences are in fact more visible with a smaller di-photon mass resolution and a larger data sample, both of which can only be achieved with upgraded detectors and/or machine. Finally, in attempting to distinguish the two approaches, we have also noted a tension in the underlying dynamics. Any distortion effect of a single BW shape can only be exploited when the mass difference is sufficiently larg