Interpreting 750 GeV Diphoton Excess in Plain NMSSM

NMSSM has enough ingredients to explain the diphoton excess at 750 GeV: singlet-like (pseudo) scalar ($a$) $s$ and higgsinos as heavy vector-like fermions. We consider the production of the 750 GeV singlet-like pseudo scalar a from a decay of the doublet-like pseudo scalar $A$, and the subsequent decay of $a$ into two photons via higgsino loop. We demonstrate that this cascade decay of the NMSSM Higgs bosons can explain the diphoton excess in the 13 TeV data consistently with the absence of a significant excess in the 8 TeV data.


Introduction
Recently ATLAS and CMS have reported excesses in the diphoton mass distribution around m γγ 750 GeV in their 13 TeV data. The local significance assuming narrow width is ∼ 3.6 σ for ATLAS [1] and ∼ 2.6 σ for CMS [2]. ATLAS and CMS have presented their updated analyses at Moriond conference. With the improved analyses, the local significance has increased to ∼ 3.9 σ and ∼ 3.4 σ for ATLAS and CMS, respectively [3,4]. Fitting that excess with a narrow resonance around 750 GeV, CMS reports for the cross section times branching ratio, σ 13TeV · BR γγ , the value 2.6 ÷ 7.7 fb at 1 σ and 0.85 ÷ (11 − 12)fb at 2 σ (see Fig.10 of [4]). The CMS fit of the excess around m γγ 750 GeV in the 8 TeV data gives 0.31 ÷ 1.00 fb at 1-σ and 0.06 ÷ 1.45 fb at 2-σ [4] 1 . The ATLAS collaboration has not provided such a detailed analysis for a narrow resonance hypothesis. A fit reported in ref. [5] gives for σ · BR γγ the values 4 ÷ 7 fb and 0 ÷ 0.42 fb at 1 σ at 13 TeV and 8 TeV, respectively. No information about the 2 σ regions is available.
The possible interpretation and implications of the excess has been intensively studied. Most such studies introduce new particles to account for the excess without asking about their UV origin, and interpretation within the known models in particular Minimal Supersymmetric Starndard Model (MSSM) and Next-to-Minimal Supersymmetric Standard Model (NMSSM) is rare. 2 In this paper we study the possibility to explain the diphoton excess within the framework of NMSSM without introducing additional particles.
One of the most straightforward interpretations of the excess is to consider a direct production of a scalar or pseudoscalar 750 GeV particle, X, decaying to two photons: αβ → X → γγ, where α, β are the initial state partons. If the model is renormalizable, X → γγ suggests the existence of electrically-charged vector-like fermions (or scalars) coupled to X , which generate the effective operator XF µν F µν (F µν ). Such fermions should be heavier than m X /2 375 GeV, otherwise the diphoton rate is strongly suppressed because X predominantly decays into the vector-like fermions on shell. Similar argument disfavours the possibility to identify X as the heavy higgs bosons in the MSSM or 2HDM, 3 because in such models X predominantly decays into tt and/or bb [71]. In general, in such scenarios the decay branching ratios of X are strongly correlated with the production cross section.
Another possibility is to consider the production of X from a decay of a heavy resonance Y r associated with another particle [10,[72][73][74][75]. This topology has two advantages. First, BR(X → γγ) is independent of the production cross section of the resonance. This is not the case for the previous topology, because a large production cross section leads to a large rate of the inverse decay process X → αβ, which suppresses BR(X → γγ). Second, the mass of Y r has to be larger than m X 750 GeV, and the 13 1 In Fig.10 of [4], CMS rescaled the fitted cross section of the 8 TeV result to 13 TeV assuming the gg initial state. We rescale this back into 8 TeV.
2 For R-parity violating (RPV) MSSM see [6,7] and for NMSSM with pp → H → aa → (γγ)(γγ) see [8,9]. 3 See however [69,70]. TeV production cross section of Y r is more enhanced with respect to the 8 TeV cross section, compared to the previous topology. In this context we notice that, while there is no big tension between 8 and 13 TeV data in the CMS fits interpreted as a direct production of a 750 GeV resonance, the fit of ref. [5] to the ATLAS data shows such a tension well above 2 σ level. For instance, if the initial partons are gluons, αβ = gg, translating the results of that fit for 8 TeV, interpreted as a direct production of the 750 GeV resonance, to 13 TeV clearly shows the problem. Thus, the cascade topology may slightly help to reconcile the ATLAS data at 8 and 13 TeV and the results of both experiments. This topology can be relatively easily realised in the NMSSM by identifying Y r = A, Y d = s and X = a: αβ → A → sa, a → γγ, as shown in Fig. 1, where A is the doublet-like pseudo scalar and (a) s is the singlet-like (pseudo) scalar. In NMSSM a → γγ is induced by a higgsino loop diagram also shown in Fig. 1. The Y d = h is disfavoured because non-zero Aha coupling requires doublet-singlet mixing in the pseudo-scalar sector (Aa mixing), suppressing a → γγ branching ratio. In our scenario, s predominantly decays into bb through a mixing with H. Although the current data would not have enough sensitivity to discriminate these extra jets from other jets with QCD origin, this scenario can be tested by looking at these b-jets in the future analysis.
The paper is organised as follows. In section 2 we demonstrate our scenario in a simplified framework in which the mixing between singlet and doublet states are ignored. In section 3 we consider how our scenario can be realised in the NMSSM taking the effect of mixing into account. We conclude this paper in section 4.

Interpretation with pure states
We first discuss our scenario in a simplified framework where the resonance A is pure doublet state and the lightest CP even and odd Higgs bosons, s and a, are exclusively originated from   the singlet field S. The signal of the diphoton excess is given by where the cross section σ(pp → A) depends on the centre of mass energy of the proton-proton collision. Fig. 2 shows the NLO production cross section of A from the bb (red) and gg (blue) initial states as a function of m A for √ S = 13 (solid) and 8 (dashed) TeV. In the left (right) panel of Fig. 2, the thick and thin lines correspond to tan β = 50 and 30 (1.5 and 3), respectively. The cross sections are calculated using SusHi v.1.5.0 [76][77][78][79][80][81][82][83]. As can be seen, the √ S = 13 TeV production cross section for large (small) tan β values is dominated by bb (gg) initial state. It can be as large as 400 (200) fb for tan β = 50 (30) at m A ∼ 850 GeV. The cross section enhances from 8 TeV to 13 TeV by factor of 5 for gg and 6.7 for bb initial states. We impose the m A dependent upper limit on σ(bb → A) · BR(A → τ + τ − ) and σ(gg → A) · BR(A → τ + τ − ) obtained from the 8 TeV CMS search for the neutral Higgs boson decaying to di-tau [84]. We found the m A < ∼ 840 GeV is excluded by this constraint for bb initial state at tan β = 50 and this region is not shown in Fig. 2. For tan β = 30 or the gg initial state, the whole region with m A > 750 GeV is allowed.
We define the interaction between A-s-a as With the coupling g Asa , the partial decay rate of A → sa is given by whereλ In what follows we assume m s = 65 GeV and 815 ≤ m A ≤ 875 GeV. In this parameter region, h → ss and A → ha are kinematically forbidden while A → sa is allowed. The A → sa decay mode competes with A → bb and A → tt in the large and small tan β regimes, respectively. The partial decay rates are given by (4) The decay modes into gauge bosons are highly suppressed due to the CP property. Fig. 3 shows the branching ratio of A → sa for m A = 850 GeV, m s = 65 GeV as a function of tan β and |g Asa |/(246 GeV). At a fixed g Asa , BR(A → sa) is maximised around tan β ∼ 7. This is because the decay rate of A → ff is minimised in this region. For small ( < ∼ 2) and large ( > ∼ 30) tan β, We focus on the process in which a decays to two photons through higgsino loop. 4 If a is pure singlet and the gauginos are decoupled, BR(a → γγ) does not depend on the higgsino mass nor the ah +h− coupling, and is entirely determined by quantum numbers of higgsinos. The branching ratios are given as fit of ref. [5] to the ATLAS data are not on equal footing, we do not average their results and discuss them in turn. The results for the coupling g Asa based on the CMS analysis are summarized in the left panel of Fig. 4. The blue region is favoured by the 13 TeV excess at 1 σ level, (σ · BR) signal 13TeV ∈ [2.6, 7.7] fb, and the yellow one by the 2 σ range [0. 85,12] fb. The green region is favoured by the excess in the 8 TeV data at 1 σ level, (σ · BR) signal 8TeV ∈ [0.31, 1.00] fb. The grey region corresponds to (σ · BR) signal 8TeV > 1.45 fb which is disfavoured at 2-σ at 8 TeV. As can be seen, there exist two favoured regions, (a) small ( < ∼ 2) tan β region and (b) large ( > ∼ 20) tan β region. This is because the production cross section, pp → A, is maximised for these two regions. In the small tan β region gg → A via the top-quark loop dominates the production processes, whereas bb → A is dominant in the large tan β region. The enhancement in the cross section compensates the slight suppression in BR(A → sa) (see Fig. 3). For moderate tan β, the signal event rate cannot be large enough to be within the 1 σ regions due to the small cross section even for |g Asa |/(246 GeV) > ∼ 1.5 where the BR(A → sa) is already saturated BR(A → sa) ∼ 1 and increasing g Asa further does not help to enhance the signal rate. As can be seen, both favoured regions require relatively large g Asa coupling. For large and small tan β regions, the 1 σ region requires |g Asa |/(246 GeV) > ∼ 1 and 2, respectively. In the right panel of Fig. 4 we show the results for the coupling g Asa based on the fit of ref. [5] to the ATLAS data. We see that the tension between the 13 and 8 TeV data does not disappear even with the cascade decay topology, where the primary object has the mass of 850 GeV, and remains at the level of approximately 2 σ. 5 The excess at 13 TeV requires, at 1 σ, somewhat larger values of the coupling g Asa .
In the simplified framework discussed so far, the dominant decay mode of s becomes s → γγ because other gauge boson final states are not kinematically allowed. This will cause a strong tension with the fact that ATLAS and CMS did not observe extra photons other than the diphoton excess with m γγ 750 GeV. However, this problem can be easily circumvented by introducing a mixing between s and H. With this mixing s will dominantly decay into bb.

Realisation in NMSSM
The superpotential and soft SUSY breaking Lagrangian of the NMSSM are given by (c.f. [86]) where we assume all couplings are real. 6 Notice that the MSSM µ-term, W MSSM ⊃ µ MSSM H u H d , can be removed by redefining S by a constant shift. We fix S in this way, hence µ MSSM = 0. We first rotate the doublet Higgs bosons H u and H d by the angle β and define the new field basisĤ = sin βH 0 dR − cos βH 0 uR ,ĥ = cos βH 0 dR + sin βH 0 uR ,ŝ = S R , A = sin βH 0 dI + cos βH 0 uI ,Ĝ = cos βH 0 dI − sin βH 0 uI ,â = S I .
By this rotation,Ĥ does not have the vacuum expectation value, andĜ becomes the Goldstone boson eaten by Z. The scalar mass eigenstates, denoted by h i (with h i = h, H, s where h is the SM-like Higgs), are expressed in terms of the hatted fields with the help of the diagonalization matrixS: The pseudoscalar mass eigenstates, A and a, are related to the hatted fields,Â andâ, by a rotation by angle θ Aa . TheÂ-ŝ-â interaction is given by the F-term of S as | ∂W In the previous section we mentioned that one should allow theĤ-ŝ update for this number has been given after Moriond conference and one cannot infer it from the fit of ref. [5]. We note that for 850 GeV resonance produced from bb initial state the increase of the cross-section from 8 to 13 TeV is about 40% bigger than that for 750 GeV resonance produced from gg initial state. In the right panel of Fig. 4 we see that, indeed, the compatibility between the fits of ref. [5] to 8 and 13 TeV data sets is somewhat better at large tan β than at its small values. 6 We use general NMSSM Lagrangian without imposing Z 3 or scale invariance. This version of NMSSM has various phenomenological advantages. See e.g. [87,88]. mixing in order to suppress unwanted s → γγ decay. Neglecting A-a mixing, the coupling g Asa is given as In the previous section we have shown that |g Asa /v SM | > ∼ 1 (2) is required for tan β > ∼ 20 ( < ∼ 2) (See Fig. 4.). Clearly, one needs the product |λκ| > ∼ 1 (2) for large (small) tan β to explain the excess. Such large values of λ and/or κ indicate the Landau pole at the scale µ UV much below the GUT scale. In Fig. 5 we show the constraint from the Landau pole. If our topology is responsible for the observed diphoton excess, this indicates the existence of the UV cut-off typically of the order of 100 TeV.
Dropping the Goldstone mode, the entries of the mass matrix for the pseudo-scalar sector (A, a) are given by where B eff ≡ A λ + κv s andm 2 3 ≡ m 2 3 + λ(µ v s + ξ F ). The m 2 3 is the soft breaking mass term L MSSM soft ⊃ m 2 3 H u H d , v s ≡ s ≡ µ/λ and ∆ 2 AA/aa/Aa are the radiative corrections.
The mixing between A and a is determined by where we used m A = 850 GeV, m a = 750 GeV. This mixing strongly affects the BR(a → γγ) because it introduces a → bb and tt modes through the mixing. The reduction of the signal strength can be parameterized by r as For | sin θ Aa | 1, r can be written as where Γ A ff is the sum of the partial decay rates in Eq. (4) at m A = 750 GeV and Γ a V V is the sum of the partial decay rates of the pure state a into W + W − , ZZ, Zγ and γγ, which can be written as The factor |λ| 2 can be understood because ahh coupling is given by λ √ 2 . The f (mh) is obtained from the higgsino loop diagram and we find f (mh) 1.5 · 10 −2 GeV for mh = |µ| 375 GeV. The condition r > ∼ 0.5 can be translated as for large ( > ∼ 10) or small ( < ∼ 2) tan β. This puts a strong constraint on the parameters appearing in Eq. (15).
In the scalar sector (Ĥ,ĥ,ŝ), the elements of the mass matrix are given by where Λ ≡ B eff + κv s + µ = A λ + 2κv s + µ and (δm 2 h ) rad is the radiative correction induced by the stop loop. Typically, for large tan β this scenario requires heavy stops (mt ∼ O(10) TeV) depending on the size of the stop mixing parameter X t in order to achieve m h = 125 GeV. The ∆ 2 HH/hh/ss/Hh/Hs/hs are the radiative corrections contributing to the NMSSM Higgs boson mass matrices.
The elements of the diagonalization matrixS must respect various phenomenological constraints. The LEP limit on the e + e − → Z * → Zs (s → bb) process for the 65 GeV scalar gives the boundS sĥ · BR(s → bb) < ∼ 0.16 [89], where BR(s → bb) depends in principle onS sĤ mixing and tan β [90]. The measurements of the properties of the SM-like Higgs boson at the LHC also give constraints on the mixing angles. The deviation of its coupling to the gauge bosons is now constrained up to ∼ 20 % at 95 % CL [91,92]. This translates into the constraint on theS entries asS hĤ ,S hŝ < ∼ 0.2. In the parameter space relevant for our model, the elementsS sĤ andS Hŝ remain unconstrained and may be large. Neglecting the small mixing elements they may be approximated byS where for future convenience we introduced the mixing angle θ sH satisfying In the last equality we have used m H = 850 GeV, m s = 65 GeV. The two small off-diagonal entries ofS may be approximated as follows The elementsS hĤ andS hŝ are related to the above ones by orthogonality ofS: S hŝ − cos θ sHS sĥ + sin θ sHS Hĥ ,S hĤ − cos θ sHS Hĥ − sin θ sHS sĥ .
Clearly, the values of the Higgs boson masses and the constraints on the mixing angles would select some regions of the NMSSM parameter space. However, the complexity of the NMSSM Higgs potential make a full quantitative analysis of our scenario, with radiative corrections included, challenging and premature. Merely for the illustration purpose, we attempt to find the NMSSM parameters that satisfy the above conditions using approximate forms of the 1loop radiative corrections. Some attention has to be paid to the magnitude of the radiative corrections. Indeed, we note that some of the 1-loop radiative correction terms are proportional to the 3rd power of λ or κ and can be as large as the tree level terms for |λ|, |κ| > ∼ 1 [93]. The 2-loop corrections may also be large [94] in this region. 7 For large λ and κ, neglecting the corrections proportional to the gauge and Yukawa couplings, the leading terms of the radiative corrections to the off-diagonal mass matrix elements are given by 8 where It is easy to find solutions for the parameters of the model satisfying the constraints m H = m A = 850 GeV, m s = 65 GeV , µ = 375 GeV, vanishing Aa mixing (θ Aa = 0) and smallS sĥ . We used the following procedure: The scalar mass squared matrix has 6 independent parameters. We choose them as 3 eigenvalues (m 2 h , m 2 H , m 2 s ) and 3 off-diagonal entries of the diagonalization matrix (S sĤ ,S sĥ ,S hĤ ). Using this parameterization we calculate the off-diagonal elements of the scalar mass squared matrix and compare them with the same elements expressed by the parameters of the model in eqs. (23)- (25). One of the parameters, µ , is fixed by the requirement of vanishing A-a mixing: µ = A λ − 2µκ/λ. Then, for some fixed values of the elements (S sĤ , S sĥ ,S hĤ ), we are left with the set of three equations for three parameters: λ, κ and A λ . In general there is a discrete set of solutions.
In the actual numerical calculations we had to modify this simple prescription. In order to compare our results with the experimental constraints illustrated in Fig. 4 we were fixing the value of g Asa given by eq. (11). This fixes one combination of the parameters λ, κ and A λ . Thus, only two mixing elements (chosen to beS sĤ ,S sĥ ) remain as input for our calculations while the third one (S hĤ ) is obtained as output. Numerical iteration procedures are used to find solutions.
One of the input mixing elements,S sĥ , is quite strongly constrained by the LEP data. Thus, after fixing the values of the scalar masses and tan β,S sĤ remains the only input quantity which may be changed in a relatively wide range. The dependence of the results onS sĤ is shown in Fig. 6 for the example with tan β = 20,S sĥ = 0 and |g Asa | = 0.6. One can see that λ increases withS sĤ while |κ| has a minimum. The behavior of λ follows from the fact that for bigger mixingS sĤ one needs bigger M 2 Hŝ which grows with λ (at least the tree contribution, see eq. (24)). Then the behavior of κ follows from relation (11). The leading (in λ and κ) loop correction to the Higgs mass is a quite complicated function of the parameters. From the right panel in Fig. 6 one can see that it may even vanish for some combination of λ and κ but generally is an increasing function of the input mixing parameterS sĤ . Examples presented in Fig. 6 (and in Table 1) were obtained forS sĥ = 0. We checked that the results do not change substantially for the values ofS sĥ allowed by the LEP data.  Table 1. For large tan β the values ofS sĤ are chosen to give |κ| close to the smallest possible (for a given set of other parameters) value in order to get the Landau pole scale as big as possible. For small tan β we have to choose much smaller S sĤ in order to avoid huge tree level contribution to the Higgs mass (value of λ increases with S sĤ ). The first example in Table 1 shows thatS sĤ > ∼ 0.1 can easily lead to too large m tree h for tan β = 2. The last two columns of Table 1 show the SM-like Higgs boson at the tree level, m tree h , and with the leading (for large λ and κ) loop corrections, m h (but before including the radiative correction from the scalar top loop). An interesting observation is that in the parameter range selected by the constraints of very smallĥ-ŝ andÂ-â mixings the radiative corrections to the Higgs potential from the NMSSM Higgs bosons are actually small, in spite of the sizable values of λ and, particularly, κ. This is related to the fact that some of potentially large contributions are proportional to appropriate mixing elements and are small in the limit of small mixings. Thus, the values of m h given in Table 1 are almost entirely controlled by the tree-level effects. The mixing elements, other thanS sĤ ≈ −S Hŝ , are small onceS sĥ is taken to be small (to fulfill the LEP constrains).S Hĥ is suppressed by m 2 H (see eqs. (28)) and typically is below 0.01. The two remaining off-diagonal elements are also small due to relations (29).S hŝ ≈ −S sĥ up to small corrections while |S hĤ | < 0.1 (< 0.01 in most cases). All these mixing elements are well below present experimental bounds. The numbers given in the table illustrate the expected order of magnitude for the soft mass parameters necessary to explain the di-photon excess in our scenario and indicate that it will be fine-tuned at the level of 1 per mille.
Finally we comment on the constraint from electroweak precision tests. It has been pointed out [95][96][97]    higgsinos to the T -parameter [98] as a consequence of violation of SU(2) custodial symmetry. However, generically, in the selected region, λ < 1. Moreover, ref. [95] shows that even for λ = 2 there are strips around the singlino mass parameter |µ s | = |µ + (κ/λ)µ| 750 ÷ 800 GeV where the higgsino contribution to the T -parameter vanishes or is negative independently of tan β and weakly dependent on the value of µ. It is not difficult to find solutions with the singlino mass in the above range, as for instance the last example in Table 1. We, therefore, expect the higgsino contribution to the T parameter not to be a problem for our scenario. One can also expect some cancellation between the higgsino contribution and the contributions from NMSSM Higgs bosons. We leave a detailed numerical analysis for future work.

Conclusions
We demonstrate that the plain NMSSM can explain the observed diphoton excess at m γγ 750 GeV as a decay of a single particle into two photons at the price of a relatively low UV cut-off (around 100 TeV) and of a certain fine tuning of the parameters. The mechanism behind this scenario is production of a doublet-like pseudo scalar A, decaying into a singlet-like pseudo scalar a, which subsequently decays via the vector-like higgsino loop into two photons. The predicted width of a is very small, much below the experimental resolution. The two-photon signal should be associated with b-quark jets coming from the decay A → as, with s decaying dominantly into a pair of b quarks. The pseudo scalar a decays also into other channels with the branching ratios given by Eq. (5).
The topology proposed in this paper is the only one that can explain the 750 GeV excess in the plain NMSSM due to a single particle decay. Another possibility for the NMSSM, recently proposed, is to explain the observed signal by the decays of two light pseudo scalars, to two collimated photons each. The latter interpretation could explain a broad peak at 750 GeV, if confirmed experimentally. The width of the signal will give a crucial discrimination between different proposed interpretations, in particular between perturbative and non-perturbative scenarios.