The NMSSM lives: with the 750 GeV diphoton excess

We propose an NMSSM scenario that can explain the excess in the diphoton spectrum at 750 GeV recently observed by ATLAS and CMS. We show that in a certain limit with a very light pseudoscalar one can reproduce the experimental results without invoking exotic matter. The 750 GeV excess is produced by two resonant heavy Higgs bosons with masses ∼750\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim }750$$\end{document} GeV, which subsequently decay to two light pseudoscalars. Each of these decays to collimated photon pairs that appear as a single photon in the electromagnetic calorimeter. A mass gap between heavy Higgses mimics a large width of the 750 GeV peak. The production mechanism, containing a strong component via initial b quarks, ameliorates a possible tension with 8 TeV data compared to other production modes. We also discuss other constraints, in particular from low-energy experiments. Finally, we discuss possible methods that could distinguish our proposal from other physics models describing the diphoton excess in the Run-II of the LHC.


Introduction
The ATLAS [1] and CMS [2] experiments at the Large Hadron Collider (LHC) have both reported an excess in the diphoton channel at an invariant mass of about 750 GeV corresponding to a local (global) significance of 3.6 σ (2.0 σ ) and 2.6 σ (1.2 σ ), respectively. The result is of course not conclusive, but if the excess were confirmed, this would be the first sign of new physics at the terascale energies.
The simplest explanation requires the production of an s-channel spin-0 or spin-2 resonance according to the Yang-Landau theorem [3,4]. The observed cross section of roughly O(10) fb is relatively large and thus it is natural to assume that the new resonance is produced via the strong interaction and have a large decay rate into diphotons. A light quark initiated resonance would be in severe tension with the LHC Run-I, a e-mail: krolb@fuw.edu.pl since the parton luminosity ratio between √ s = 13 TeV and 8 TeV is relatively small for light quark initial states. As a consequence, the resonance would have been observed at LHC Run-I if it were produced via quark-antiquark initial states. For a gluon induced resonance the tension with 8 TeV is reduced but still significant [5]. On the contrary, associated production with b quarks does not suffer from 8 TeV constraints. Moreover, the reported event topology is consistent with the single production of a resonance, i.e. non-resonant production of the 750 GeV particle in a cascade decay [6] is disfavored since no additional activity was observed in the peak-region events. 1 Finally, the apparently large width of around 45 GeV, preferred by ATLAS, points to large couplings to its daughter particles. However, strict constraints exist on decays of heavy resonances into electroweak gauge bosons and light Standard Model (SM) fermions and thus the resonance should decay into final states which evade all current experimental searches, implying e.g. a large invisible decay rate or decays to quarks and gluons. Another way out (which we actually consider in this paper) would be the presence of two overlapping resonances with narrow widths which allow one to explain the large width within the current experimental accuracy [8].
The observed diphoton rate cannot be explained with a SM-like Higgs boson because its tree-level decays into third generation quarks and/or gauge bosons are too large compared to the loop induced decays into diphoton final states. However, simple extensions of the SM Higgs sector such as a singlet extension or Two-Higgs-Doublet Model (2HDM) are also plagued with too small diphoton rates and the way out is to introduce new vector-like fermions: see e.g. [9,10]; see [11] for an overview. There are only a few phenomenologically viable explanations within the framework of supersymmetry (SUSY) [12]. It seems to be impossible to find a solu-tion within the Minimal Supersymmetric Standard Model (MSSM) [13]. New vector-like fermions and singlets have to be added to the particle spectrum of the MSSM in order to explain the diphoton excess [14][15][16][17][18]. Other SUSY solutions either involves R-parity violation [19,20] or assume a very low SUSY breaking scale [21][22][23].
The Next-to-Minimal Supersymmetric Standard Model (NMSSM) -see [24,25] for a review -is a well-motivated supersymmetry-inspired extension of the SM. Beyond the elegant features of supersymmetric models in view of the hierarchy problem or one-step unification, and their potential in terms of Dark Matter (DM), the original purpose of the NMSSM rests with the 'μ-problem' [26] of the simpler MSSM: this issue is addressed via the addition of a singlet superfield to the matter content of the MSSM, the 'μ'-parameter then being generated dynamically when the singlet takes a vacuum expectation value (v.e.v.). Additionally, the NMSSM has received renewed attention ever since the Higgs discovery in the Run-I phase of the LHC [27,28], due to its interesting features in terms of a supersymmetric interpretation of the observed Higgs signals -see [29] for a recent analysis and list of references. While several versions of the NMSSM can be formulated, we will focus here on the simplest one, characterized by a Z 3 -symmetry and CPconservation. Let us stress here that we will not include new exotic matter but rely strictly on the simple matter content of this model.
Our purpose in this paper is to present a phenomenologically viable scenario accounting for the diphoton excess at ∼750 GeV in the context of the NMSSM. This explanation rests on the possibility that a pp → → 2( → γ γ ) process -see Fig. 1 -could not be distinguished from a diphoton signal in the experimental searches [30][31][32][33][34][35][36]. The NMSSM Higgs sector then offers a suitable framework to embed this topology: can be identified with a very light CP-odd singlet decaying dominantly into a diphoton pair, while heavy CP-even doublet and the CP-even singlet have to be combined to mimic a 750 GeV resonance with adequate properties. Note that contrarily to the proposals which we mentioned earlier, the mechanism that we consider does not rely on the ad hoc inclusion of additional matter (e.g. vectorlike fermions) and uses only the existing features and degrees of freedom of a rather simple and well-motivated model. 2 While our project was in its finalizing stages, we became aware of another recent proposal to explain the diphoton signals within the NMSSM [37], which shares some traits with our interpretation but also differs in several respects. The diphoton decay of the pseudoscalar in [37] relies on a substantial 'quasi-mixing' of this Higgs state with the η-meson: 2 We observe also that this mechanism can be transposed to an apparently simpler -in fact theoretically less motivated -singlet extension of the (Type II) 2HDM.  Fig. 1 The resonant production of followed by the decay to two scalars and photons. The final-state photons are pairwise collimated this requirement induces substantial limits from ϒ-decays and results in a quite constrained regime. In our case, we consider the mixing with the π 0 and estimate this effect more quantitatively using the chiral perturbation theory for pions. Moreover, we shall propose a wider selection of benchmark points, illustrating the flexibility of the mechanisms that we employ. Furthermore, we shall analyze in further detail how our scenario compares to the experimental data and study complementary signatures.
In the following section, we will detail how the NMSSM provides an acceptable framework for the ∼750 GeV signal and address several phenomenological issues which constrain the parameter space. We will thus propose specific examples of NMSSM points satisfying these requirements before discussing their relevance in fitting the diphoton excess in Sect. 3. Then we consider other experimental signatures that are specific to our model and have promising prospects in the Run-II, including the decays of the light pseudoscalar to e + e − pairs. We also discuss possible experimental signatures from the higgsino, slepton, and heavy Higgs bosons sector. We will conclude with a brief summary in the last section.

Embedding the 750 GeV diphoton excess in the parameter space of the NMSSM
The NMSSM Higgs sector -see e.g. [24,25] -consists of two doublets, 174 GeV, G F representing the Fermi constant. Moreover, in order to provide a more physical handle on the parameter space, we define μ ≡ λs -which is analogous to its MSSM 'μ' counterpart and also sets the tree-level mass of the higgsino states -and M 2 A ≡ 2λs sin 2β (A λ + κs) -which sets the scale of an approximate Higgs SU (2)-doublet, generally heavy, i.e. M A > 125 GeV, if one wants to avoid the complications of light non-SM-like doublet states. 3 In the following, we shall employ the parameter set (M A , tan β, μ, λ, κ, A κ ).

Identifying the light state ' '
In order to fit the diphoton excess, we want to employ the pp → → 2( → γ γ ) topology. Obviously, would be a CP-even (in order to allow a decay to two lighter identical states) Higgs state at ∼750 GeV, while could be in principle CP-odd or CP-even: given the limited pool of Higgs states in the NMSSM and the fact that we will have another use for the CP-even singlet (see below), this state will be identified, however, with the singlet-like pseudoscalar. For its diphoton decay to be indistinguishable from a single photon, this Higgs state would have to be light enough: m 0.5 GeV [32]. We observe that this scenario with a light CP-odd state is phenomenologically viable, in view of the current status of collider constraints. However, severe limits intervene at low mass from, e.g. flavor observables or the properties of the SM-like Higgs boson H SM -a sizable branching fraction H SM → would be at odds with the measurements of ATLAS and CMS in the LHC Run-I phase. Note that the state cannot be dominantly doublet in view of the consequences for the spectrum: this would indeed imply the presence of a light charged Higgs (likely in contradiction with limits from top decays [38,39]) and a light CP-even doublet state (in tension with LEP [40,41] and likely to open up sizable unconventional decays of H SM ).
We complete this discussion with a remark concerning the naturalness of a light CP-odd Higgs in the NMSSM: in two specific limits, this particle appears as the pseudo-Goldstone boson of an approximate symmetry of the Higgs sector: • for κ → 0, i.e. κ λ, the Higgs potential is invariant under a U (1) Peccei-Quinn symmetry, under which 3 A similar quantity would be the mass of the charged-Higgs state: : this is the R-symmetry limit since the Higgs charges for the specific choice Q R S = −2/3 coincide with those that these fields would receive under a genuine R-symmetry (also broken by the gaugino masses) of the NMSSM with unbroken supersymmetry. Note that with our choice of parameters, We will see that, in the scenario under consideration, the factor 1 − λM 2 A 2κμ 2 sin 2β has to be small, so that the Rsymmetry limit can be invoked.
Most of the characteristics of the NMSSM which intervene in the interpretation of the signal at ∼750 GeV and the subsequent constraints can be understood from the relations at tree level. This follows from the fact that the relevant physics is driven by comparatively heavy or singlet-like states, for which the radiative corrections are relatively mild in proportion. We shall therefore propose a discussion at tree level in the following. Note, however, that loop corrections play a crucial role for the mass of the SM-like Higgs state so that we will employ tools including leading radiative effects in the numerical analysis (see below). Of course, the exact correlations at tree level are slightly displaced by loop corrections so that small adjustments in the choice of parameters will prove necessary and the relations that we derive below should not be understood as rigid constraints, but rather as qualitative guidelines/trends.
In the base (A D , A S ), the tree-level squared-mass matrix for the NMSSM CP-odd sector reads (1) We observe that the singlet squared-mass m 2 A S is largely determined by the choice of A κ (since the other terms are small). In particular, low values of A κ ensure that the singlet is light. Note also that in this case M 2 A m 2 A S -as we will see later, M A ∼ 750 GeV. The subdominant doublet component of the light CP-odd mass state A 1 = 1 − P 2 d A S + P d A D A S + P d A D can be obtained approximately: From now on, we identify this light state A 1 to the state of the pp → → 2( → γ γ ) topology. Now we will address the decays of A 1 . For the mass range under consideration, m A 1 0.5 GeV, the interplay of the pseudoscalar with the strongly interacting sector is non-trivial. In a naive partonic approach, the diphoton decay is significantly suppressed -O(10 −5 ) -for a usual NMSSM pseudoscalar, as long as the competition of the hadronic and muonic final states remains open. Then, below three times the pion mass (the decay to two pions would violate parity), hadronic final states seem inaccessible so that the pseudoscalar would then mainly decay to muons or, below the μ + μ − threshold, to γ γ -O(70 %) at ∼210 MeV, and the e + e − final state would eventually dominate at lower masses ( 160 MeV). Moreover, in such a regime, the pseudoscalar would appear as relatively long lived below the dimuon threshold: the diphoton and e + e − final states are essentially mediated by the small doublet component of A 1 , resulting in the following (approximate) widths (for m A 1 210 MeV): The corresponding total width is thus of order A 1 ∼ (10 −13 GeV) P 2 d tan 2 β ∼ 10 −14 GeV for P d ∼ 3 · 10 −2 , tan β ∼ 10, and m A 1 ∼ 200 MeV, and narrower at lower masses. The decay length of A 1 at rest is then of order 2 cm. Including the boost factor in the considered topology would lead to a decay length of about 40 m at the LHC, i.e. to an invisible pseudoscalar.
This picture is over-simplistic, however. As [37] already noticed, hadronic effects could substantially affect the decays of A 1 . If one were to disregard the tripion threshold for strong-interacting decays, the ss and gg final states would completely dominate the partonic decays of A 1 , which highlights the sensitivity of the pseudoscalar to the strong sector. In particular, considering the couplings of the pseudoscalar to mesons, one can write the following operators: Such terms induce a mixing of the Higgs pseudoscalars with the pseudoscalar mesons or, in other words, open the possibility for A 1 → (π 0 ) * or A 1 → (η) * decays (the latter being the choice of [37]). The magnitude of the mixing elements δm 2 A 1 π 0 /η can be assessed by rewriting the partonic couplings of A 1 in terms of the axial-flavor currents and their expression in the pion model (chiral perturbation theory) [42][43][44]: where 2sγ μ γ 5 s) and the divergences of these currents in the pion model are determined by the equations of motion: Here, η 1 denotes the Goldstone field associated with the U (1) axial-flavor symmetry while π 3 and π 8 correspond to the gen- axial-flavor symmetry. Dismissing the refinements of the π 3 , π 8 , and η 1 mixings, we can identify π 3 ∼ π 0 , π 8 ∼ η, and η 1 ∼ η . This determines: We observe that these mass-mixing parameters are small: with f π 93 MeV and the typical values P d ∼ 0.03 and tan β ∼ 10, δm 2 A 1 π 0 ∼ (4 · 10 −5 m 2 π ) and δm 2 A 1 η ∼ (2·10 −5 m 2 η ). Therefore, they intervene only in the immediate mass-vicinity of the mesons, and a mixing between A 1 and the mesons of order 10 −2 thus requires a proximity in mass at the MeV level. In particular, following our analysis, the decay width mediated by a η-state for m A 1 ∼510 MeV, as considered in [37], should be completely superseded by e.g. muon decays of A 1 . On the other hand, multi-pion decays are still kinematically open for m A 1 ∼510 MeV, which may result in multi-photon final-states. We shall not elaborate further on this possibility and we now focus on the mixing of A 1 with π 0 , i.e. assume m A 1 ∼135 MeV. At first order in δm 2 A 1 π 0 , the mixing angle θ reads θ δm 2 A 1 π 0 /(m 2 π − m 2 A 1 ) and the mass shift driven by mixing is m π θ 2 2 1 − panel of Fig. 2. As an example, for P d = 0.035, tan β 10, and |m A 1 − m π | < 0.3 MeV we find θ = 10 −2 . The consequences for the decays of A 1 are sizable: via its small π 0 component, the pseudoscalar acquires pion-like decays with values about essentially to γ γ (BR[π 0 → γ γ ] ∼ 99 %), which overrule the 'pure-Higgs' decays of order 10 −15 GeV. The decay length in the lab frame for the process of Fig. 1 is then reduced to ∼ 0.5 m, placing most of the time the decay vertex to photons before A 1 reaches electromagnetic calorimeter. 4 Even shorter decay lengths are possible for larger A 1 -π 0 mixing. For the pion, on the other hand, the consequences are minimal as its natural decays are slightly decreased at the level of 0.01 % (far below the theoretical or experimental precisions) and slightly shifted by a Higgs-like width at the level of 10 −16 -10 −17 GeV (via interference terms). The 2 σ experimental error on the pion decay width to two photons of 4 % defines our upper limit on the mixing. 5 In the right panel of Fig. 2 we show contours of the constant mixing angle θ in the m A 1 -X d plane. With the above constraints, the regions shaded in gray are excluded. The central part around m A 1 135 MeV is excluded by the upper limit of the measurement of the width [π 0 → γ γ ], while the outer part by A 1 decay length at the LHC, which would be larger than ∼0.5 m. It should be noted that the Particle Data Group [47] reports BR[π 0 → e + e − ] = (6.46 ± 0.33) · 10 −8 leading to a width [π 0 → e + e − ] 5 · 10 −16 GeV, known at the 6 % level. Inclusion of radiative corrections shifts this value to BR[π 0 → e + e − ] = (7.48 ± 0.38) · 10 −8 [48], while the recent theoretical calculation gives (6.2 ± 0.1) · 10 −8 [49], so that the inclusion of new physics at this level may even be welcome: see e.g. [50]. In our case, the red contour in Fig. 2 denotes the mixing angles that could account for the discrepancy in π 0 → e + e − (considering the central values only).
In the following, we will assume that this mixing mechanism between A 1 and π 0 is responsible for the apparent A 1 width-resulting in a decay length of order 0.5 m for the topology of Fig. 1 -and a BR[A 1 → γ γ ] 99 %.

Identifying the heavy state ' '
We now turn to the CP-even sector. In the base ( the tree-level squared-mass matrix reads Footnote 5 continued a change in the measured mean decay length one could directly probe parameters of A 1 -π 0 mixing. It should be noted that, as long as M A v, the doublet sector is approximately diagonal and that the mass of the heavy doublet state H D falls close to M A . Moreover, keeping in mind that A κ is small in order to ensure a light A 1 , we observe that the mass of the CP-even singlet H S is dominated by the term 2κ λ μ, as long as this quantity is larger than O(v). In such a regime, the lightest CP-even Higgs state H 1 can thus be identified with H SM (approximately): where the singlet component is very small: S 13 1. Consequently, the H 1 behaves SM-like, in agreement with the experimental results [51].
We now want to identify the state of the pp → → 2( → γ γ ) topology with one of the CP-even states of the NMSSM. The associated mass should fall close to ∼750 GeV. Both H D and H S are potential candidates. However, the state should be sizably produced at the LHC, which implies a significant coupling to SM particles, such as gluons or quarks. In the case of a dominantly singlet state, such couplings are typically suppressed in proportion to the small doublet component of the state. On the other hand, H D has sizable couplings to tops (though tan −1 β suppressed with respect to the SM) and to bottom quarks (tan β-enhanced), so that a measurable production in gluon-gluon fusion (ggf) or in association with b's (bbh) is plausible. We observe that the production cross section of a 750 GeV H D in ggf at 13 TeV is well approximated by (0.6 pb)/ tan 2 β for tan β 15, and the bbh cross section, by (4·10 −4 pb)×tan 2 β. The production cross section would thus point toward an identification ∼ H D , which would imply M A ∼ 750 GeV.
Yet, another consideration is that should have a large decay into a pair of A 1 's. Naively, the SM fermionic final states would offer a sizable competition: with m 750 GeV and the relative fermionic couplings to their SM counterparts On the other hand, the decay width into 2 A 1 can be estimated as (still for m ∼ 750 GeV): For the heavy doublet state, the fermionic couplings are approximately determined by tan β: The fermionic channels will thus typically give a width of order 1 GeV at least, for moderate tan β = O (10). On the other hand, the leading terms in the couplings of a CP-even doublet with a pair of CP-odd singlets read O(GeV), so that the associated branching ratio is suppressed in view of the fermionic channels. On the other hand, for H S , the decay channels into SM particles are naturally suppressed while the decay into light CP-odd singlets is large: If this state H S were at ∼750 GeV, then 2 κ λ μ 750 GeV, while A κ is negligible (from the low mass of the CP-odd singlet). We thus obtain that [H S → A 1 A 1 ] is of order 1 GeV, as long as κ 0.25. Therefore, the branching ratio 6 → pleads for the identification ∼ H S , hence 2 κ λ μ 750 GeV. These two apparently conflicting requirements, ∼ H D for the production and ∼ H S for the decay, can actually be reconciled if one keeps in mind that the mass states H 2 and H 3 are in fact admixtures (essentially) of H D and H S . Provided the mixing is large, then the mass states will combine the properties of both their doublet and singlet CP-even components. Two interpretations are then possible for : • The first is that the diphoton excess corresponds to only one of the two states H 2 or H 3 , while the other would give a subdominant effect due to the details of the singletdoublet mixing. Considering the somewhat large width of order O(45 GeV) which seems associated with the excess, this interpretation tends to imply a very large , which could only be achieved with κ ≥ 1, that is, outside the perturbative regime. We will not discuss this scenario any further. 7 • The second possibility is that both states are sufficiently close in mass to appear as a single excess. The very large binning makes such a scenario plausible, even for mass differences of order O . Then, whatever the width of the physical states H 2 and H 3 , the associated signal would look like a broad resonance; this will be discussed in Sect. 3.1. From previous considerations, we see that the actual widths of H 2 and H 3 in the considered regime would be of order O(1 GeV). Their minimal mass difference (for a maximal mixing) can be inferred from the off-diagonal element of the squared-mass matrix Eq. (8): where we have used M A 750 GeV 2 κ λ μ. For κ = O(0.3) and tan β = O(10), we obtain a typical spread of 10 GeV in mass. This is the scenario we will be focusing on in the following.
From now on, we thus assume M A 750 GeV 2 κ λ μ, which causes a strong mixing between H D and H S , so that , both states having masses close to 750 GeV. 7 While our study was already in progress, we became aware of the discussion in [52].

Other phenomenological constraints
After identifying how the NMSSM Higgs spectrum could provide an interpretation of the diphoton excess via the pp → → 2( → γ γ ) topology, we will consider additional constraints and consequences on this scenario.
The SM-like properties of the Higgs state at ∼125 GeV place a major phenomenological limit on the existence of a light pseudoscalar: as a general rule, the channel H SM → A 1 A 1 should remain subdominant -compared to the SM width SM 4 MeV -or it would induce a sizable unconventional decay of the state at ∼125 GeV, which would cause tensions with the results of the LHC Run-I. In our case, with a very light A 1 being experimentally indistinguishable from a photon, the limits on H SM → A 1 A 1 are even more severe in order to avoid a large apparent decay H SM → γ γ , again in contradiction with the Run-I results. The corresponding width can be estimated as The One should observe that, typically, λ, κ = O(0.1), 2 < tan β O (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20), and S 13 , P d 1. As a result, the two 'dangerous' terms in Eq. (14) are those of the last line. Considering Eq. (2) and S 13 , P d 1, the condition on g H SM A 1 A 1 implies where we have applied the mass condition M A 750 GeV 2 κ λ μ for the final step. A first alternative would thus consist in choosing λ 3 · 10 −2 , which will also lead to κ 0.1: then, however, the decay H S → A 1 A 1 no longer competes with the fermionic decays of H D and the mixing among H S and H D is reduced. Our scenario would thus be invalidated. We will thus prefer to keep λ 0.1 and satisfy the limit from [H SM → A 1 A 1 ] with the condition κ λ sin 2β ∼ 0.5. It should be noted that this condition is actually less ad hoc than it looks at first glance: indeed, given the mass condition M A 750 GeV 2 κ λ μ, A λ − M A 2 1 − 2κ λ sin 2β , so that, together with the requirement of a small A κ , the assumption κ λ sin 2β ∼ 0.5 places us naturally in the R-symmetry limit of the NMSSM.
On the other hand, the perturbativity of the couplings up to the GUT scale approximately implies λ 2 + κ 2 0.5. Using the condition κ λ sin 2β ∼ 0.5 then places an upper bound on λ sin 2β 2 1+4 sin 2 2β . Moreover, our scenario requires that the width [H S → A 1 A 1 ] remains competitive in view of the fermionic decays of H D , [H D → tt/bb]: considering the production cross sections of H D , the efficient branching ratio BR[A 1 → γ γ ] 0.99, mediated by the pion, and the magnitude of the diphoton excess at ATLAS and CMS, the condition on κ can be lowered to κ 0.1 -instead of 0.25 as we discussed above -and translates into a lower bound of λ ∼ 2κ sin 2β 0.2 sin 2β.
A further implication of κ λ sin 2β ∼ 0.5 together with the mass condition M A 750 GeV 2 κ λ μ reads μ ∼ M A sin 2β. An immediate consequence is that μ ≤ 750 GeV, and even μ ≤ 375 GeV as soon as tan β 3.7, so that the higgsino states will typically be light, i.e. for tan β 3.7 H 2 and H 3 can have a non-negligible decay rate to higgsinos. This will dilute somewhat more the H 2,3 → A 1 A 1 rates, although the typical widths are of order ∼(15 GeV)λ 2 1 − λ 2 κ 2 3/2 . While nothing forbids that the gauginos also intervene at a low mass, we will assume in the following that they are heavier. A neutral higgsino is thus the lightest supersymmetric particle (LSP), which is consistent with cosmological limits on the dark matter density: the thermal higgsino relics are typically below the observed value of the relic density [12]. However, scenarios with low thermal relic density can be consistent with the measured abundance [53][54][55].
However, μ cannot be too small, since this would contradict the unsuccessful results of LEP in chargino searches [56]. Using the limit μ 100 GeV then provides a bound on tan β, in which the choice κ λ sin 2β ∼ 0.5 can be conciliated with the mass requirement M A 750 GeV 2 κ λ μ: tan β 15.
In [37], limits from ϒ decays play a determining role in constraining the parameter space: the Wilson formula [57] gives BR[ϒ(1S) → γ A 1 ] ∼ 1 · 10 −4 (P d tan β) 2 , which would typically lead to a ∼ 10 −6 -10 −5 effect in ϒ(1S) radiative decays (and similar values for ϒ(2, 3S)). However, we are not aware of experimental limits applying to a pseudoscalar with mass close to that of the pion. Searches in radiative ϒ decays typically ignore mass scales below 0.2 GeV, except for invisible final states, which are not relevant in our scenario. Similar contributions of the A 1 to the radiative decays of J/ψ are suppressed, due to the smaller charm mass and the tan −2 β suppression of the A 1 -charm coupling. Furthermore, as in [37], we find that the impact of A 1 in radiative Z -decays is orders of magnitude below the experimental bounds [47,58].
As [59] pointed out, the presence of a light Higgs pseudoscalar generically leads to tensions in flavor physics. Limits from invisible decays do not apply in our scenario, as A 1 would decay within ∼1 cm at B-factories, but the rare transitions B → K e + e − and K → π e + e − should be considered carefully. Indeed, following [60], one observes that such transitions can be mediated by a light A 1 , as BR[A 1 → e + e − ] is sizable (at the percent level) below the μ + μ − threshold. The strongest limit comes from [61]: The actual bound is in fact much stronger than one can infer from the decay rate alone. The e + e − spectrum is well measured [61], with low background and small theory uncertainty: therefore a peak in the e + e − invariant mass spectrum would be clearly visible. This implies that the A 1 contribution has to be strongly suppressed: even though, in our configuration, BR[A 1 → e + e − ] falls below the percent level due to the large pion-mediated width, the typical magnitude of the effectivebs A 1 andsd A 1 couplings, C A ∼ 10 2 -10 3 GeV 2 and C A ∼ 10 0 -10 1 GeV 2 , would result in an excess -see e.g. Eqs. (21) and (25) in [59]. The conclusion that these flavor-changing processes exclude the considered scenario would be premature, however. First, one should keep in mind that such A 1 -mediated signals may hide in the background of the pion, due to the proximity in mass. Then the actual size of C A and C A depends on the details of the sfermion spectrum and, in particular, these coefficients vanish in the super-GIM limit [60]. The sfermion sector is largely free in what precedes our analysis: its only role so far was to ensure the correct mass for the SM-like state via radiative corrections -this essentially translates into a scalar top spectrum of a few TeV, or very large mixing in the stop sector. We check that C A and C A can be made arbitrarily small for suitable choices of the squark spectra, so that flavor constraints -and not only those involving a A 1 → e + e − decay -can be generally circumvented; see benchmark point P2 below.
More general limits on the spectra, such as B → X s γ or B s → μ + μ − , should also be considered [62], but note that the already large m H ± 750 GeV, the moderate value of tan β, the flexibility of the squark spectra, and the fact that A 1 is off-resonance contrive to place our scenario within 95 % of these flavor constraints.
In Ref. [59] it was argued that important constraints on Higgs-like pseudoscalars can be obtained in beam dump experiments. The most sensitive one is the CHARM search for axions [63]. Using Eq. (3) of Ref. [63] one can estimate the F X parameter to be F X < 10 GeV, assuming that [A 1 → γ γ ]/ [π 0 → γ γ ] > 10 −4 as required by the decay length 0.5 m of A 1 at the LHC. Looking at the exclusion plot in Fig. 4 of [63] one immediately sees that this is way below sensitivity of the experiment. In any case, with this decay length at the LHC, the decay length at CHARM would be of order O(1) cm. After 60 cm the flux suppression would be ∼ 2 60 10 18 . Taking into account that the CHARM detector was 480 m away from the target and that initial A 1 flux was < 10 17 one clearly sees that possibly no pseudoscalars could have reached the detector.
Finally, the anomalous magnetic moment of the muon may also be of relevance: a light A 1 is indeed known to widen the discrepancy between the prediction of the model and the experimental measurement; see [64]. Yet, the moderate value of tan β and the presence of light higgsinos concur to make the supersymmetric corrections to (g − 2) μ the dominant new-physics effect. The overall contribution thus improves the agreement with the BNL measurement as compared to the SM. Placing (g − 2) μ within 2 σ of the experiment remains problematic, however, and can be achieved only in the upper reach of tan β ∼ 15. It should be noted that this observable depends on the details of the slepton masses, such that lighter smuons and sneutrinos would improve the situation. The LHC searches on light neutralinos and sleptons will be discussed in the following section.

Favored parameter space
To summarize this analysis, it appears that most of the parameters in the NMSSM Higgs sector are fixed or bounded in the scenario that we consider: • tan β 15 is constrained by the lower bound on chargino searches μ 100 GeV, as a result of the various correlations; note that tan β = O(10) satisfies the requirements on the fermionic decays of the states at ∼750 GeVwhich should remain moderate; • A κ O(0.1) GeV conditions a light CP-odd singlet; note that, together with the requirement A λ → 0, which, in our scenario, follows the assumptions on κ, λ, μ, and M A , A κ → 0 places us in the approximate R-symmetry limit of the NMSSM, and that A 1 thus appears as the pseudo-Goldstone boson of this R-symmetry.
Moreover, the requirements of a ∼125 GeV mass for the SM-like Higgs state and flavor physics constrain the squark spectra, while (g − 2) μ and slepton searches impact the slepton spectrum. We stress that the singlino and higgsino masses are essentially determined by the choices in the Higgs sector and that light higgsinos (constituting the LSP in the simplest configuration) appear as a trademark of this scenario.
It should be noted that most of the properties of the Higgs sector can be transposed to a simpler singlet extension of the 2HDM, without the complications in the supersymmetric spectrum and with increased number of parameters and degrees of freedom. In the latter setup, one should consider the low-energy constraints more carefully, however, as charged-Higgs contributions to B → X s γ can no longer be balanced by the SUSY loops and the anomalous magnetic moment of the muon would suffer from the negative contributions driven by the loop involving the muon and the light pseudoscalar if the 2HDM is of Type II (which determines the production of the states at ∼750 GeV). It should be kept in mind that the singlet + 2HDM framework receives no deep motivation from the hierarchy problem, DM or gauge unification.
Naturally, certain attractive features of the NMSSM Higgs sector, such as the possibility of a light CP-even singlet, appear as a necessary sacrifice in order to conciliate an interpretation of the ∼750 GeV excess with the parameter space and constraints of the NMSSM. Moreover, it could be argued that the mechanisms which we invoke -from the sizable singlet-doublet mixing at ∼750 GeV, or the condition of a A 1 -π 0 interplay, to the collimated diphoton decays, indistinguishable from a single photon -are quite elaborate. Still, it is remarkable that all the necessary properties to fit the signal can be united in a phenomenologically realistic way within as theoretically simple a model as the NMSSM, without e.g. requiring additional ad hoc matter.

Benchmark points
To investigate the NMSSM parameter space more thoroughly than the derivation at tree level allows, and account for e.g. higher-order corrections or verify various phenomenological constraints, we use the public package NMSSMTools 4.8.2 [65][66][67]. The Higgs spectrum is computed with precision setting 2, i.e. including full one-loop, Yukawadriven two-loop as well as pole corrections [68]. Note that we dismiss the width and branching fractions computed by this code for the light A 1 as they do not implement the effect of hadronic states. We simply tune the mass m A 1 to ∼135 MeV and then invoke a mixing with the pion at the level of θ ∼ 10 −2 -in practice, this might require further adjustment in the choice of m A 1 but this can be achieved with completely negligible consequences for the rest of the spectrum. Additionally, NMSSMTools is interfaced [29] with HiggsBounds 4.2.1 [69][70][71][72] and HiggsSignals 1.4.0 [73] in order to test the properties of the Higgs sector in view of current collider limits. However, we still have to check 'by hand' that the decay H 1 → A 1 A 1 does not induce a large apparent H 1 → γ γ branching ratio. The sparticle decays are obtained with NMSDecay [74,75] and the Higgs production cross sections at the LHC are obtained with SusHi 1.5.0 [76][77][78][79][80][81][82][83][84][85], interfaced with LHAPDF 5.9.1 [86] and using MSTW parton distribution functions (PDFs) at NNLO [87]. At the outcome of this search, phenomenologically realistic points exhibiting the characteristics which we described above are obtained and presented in Table 1.
In Table 1, we provide the input for NMSSMTools as well as relevant masses. The squark soft mass parameters are all chosen degenerate as mQ (for simplicity). So are also the soft masses of the sleptons, mL . We observe that the Higgs mass predictions of NMSSMTools for very heavy sfermions may not be entirely reliable, as a resummation of large log( mQ m t ) may be necessary: such effects are addressed, e.g. in [88] but the details of the correspondence between the parameters and the spectrum are of secondary importance for our conclusions. Note that all the considered points satisfy the phenomenological tests implemented within NMSSMTools -except maybe (g − 2) μ (satisfied for P4 and P9) -and HiggsBounds within 2 σ . We make sure that BR[H 1 → A 1 A 1 ] < 1 · 10 −4 . The fit values to the Higgs measurement at ∼125 GeV obtained with HiggsSignals are all competitive with the SM -i.e. χ 2 < χ 2 SM or |χ 2 /χ 2 SM − 1| 1. Concerning the flavor constraints associated to A 1 → e + e − or, more generally, to a mediation by A 1 , we stress that the proximity in mass of A 1 to the pion would certainly require a more careful analysis on the experimental side, so that the current experimental limits are likely not to apply. Yet, for completeness, we wish to show to which extent A 1 -mediated contributions to rare B and K decays can be reduced in our scenario. We thus undertook the task of tuning the parameters in the sfermion sector in order to suppress the effective flavor-changing couplings for point P2 only: it is quite clear that such a requirement can always be applied independently of the properties of the Higgs states. For this point, the trilinear Higgs-stop coupling A t is adjusted in such a way that the effective flavor-changing A 1 couplings amount to C A 1.3·10 −6 GeV 2 and C A 2.6·10 −8 GeV 2 . These extremely suppressed numbers come at the price of the 7-digit precision in the value of A t . We note that to simultaneously maintain m H 1 125 GeV a significantly heavier squark sector (compared to P1) becomes necessary: this is not unexpected as it is unlikely to combine a maximal stop mixing (which provides a large contribution to m H 1 ) and minimal effective flavor-changing A 1 couplings with the sole handle of A t . 8 The contribution to BR[B 0 → K 0 e + e − ] is then at the level of ∼ 10 −20 and, for BR[K + → π + e + e − ] at the level of ∼ 10 −23 . Such effects are far too small to be measurable experimentally. However, if we consider P1 for comparison, where the sfermion sector was not tailored to accommodate these channels, C A ∼ 130 GeV 2 and C A ∼ 3 GeV 2 , contributing to the branching ratios at the level of 10 −4 and 10 −7 , respectively, i.e. far beyond existing limits. Of course, the precision that we requested in the suppression of C A and C A for point P2 is not really necessary in view of the e + e − channels: C A ∼ 1 GeV 2 and C A ∼ 0.1 GeV 2 would be sufficient and can be achieved with the simpler requirement A t −8.5 TeV for P2. Furthermore, experimental cuts would typically require m e + e − > 140 MeV [61], so that our scenario with m A 1 m π 0 is not affected by these limits in general. However, we wish to stress that any bound on flavor transitions mediated by A 1 could be circumvented in such a fashion. Therefore, we will pay no further attention to these flavor limits on the basis that they strongly depend on the details of the sfermion sector.
We now discuss the phenomenology of the benchmark points. tan β ranges from 4 to 15 in Table 1: this implies a variety of regimes for the points under consideration, as we will see later. In particular, this determines the value of μ, i.e. the higgsino spectrum: the lightest neutralino mass varies between 100 GeV at tan β = 15 (P4) to 370 GeV at tan β = 4 (P7). κ and λ are always of order 0.1 and their ratio also depends on tan β (see previous section). Note that larger values of λ are accessible but tend to result in too efficient a cross section for the diphoton signal, as we shall discuss later. The values of A κ appear with a sizable number of digits: this corresponds to the precision necessary to keep m A 1 within 135 ± 0.5 MeV. M A falls within 10 to 60 GeV of 750 GeV. The squark masses are chosen, together with the trilinear coupling A t , so as to generate a mass close to ∼ 125 GeV The only quantity deserving discussion at this level is the mass-splitting between H 2 and H 3 : it ranges from ∼10 to ∼30 GeV, depending on the values of κ and tan β. As we aimed at a large singlet-doublet mixing close to 50 %, this mass gap is essentially determined by Eq. (12).
Finally, we indicate the magnitude of the doublet component in A 1 , P d , which plays a central role for the characteristics of this state. It follows the approximate rule P d = 0.232 · λ, that one can infer from the tree-level definition, Eq. (2), together with the various conditions on M A , μ, κ, and sin 2β presented in Sect. 2.4.
In Table 2, we indicate several Higgs branching fractions as well as the production cross sections in ggf and bbh at the LHC for the heavy states. First, we check that BR[H 1 → A 1 A 1 ] < 1 · 10 −4 , so that no non-SM diphoton decay of H 1 conflicts with the LHC Run-I results -the total width of H 1 is always SM-like: H 1 4 · 10 −3 GeV. Concerning the heavy CP-even states, their widths fall typically between 1 and 2 GeV and are dominated by the fermionic channels bb and/or tt depending on tan β. BR[H 2,3 → A 1 A 1 ] is typically at 10-30 % and comparable for both states (due to an efficient mixing): we observe that larger values ( 50 %) are accessible for larger κ (λ) but such values tend to overshoot 7 · 10 −6 6 · 10 −5 1 · 10 −5 9 · 10 −6 3 · 10 −6 3 · 10 −7 8 · 10 −6 1 · 10 −6 5 · 10 −6 H2 (GeV) the magnitude of the observed diphoton cross section. We also provide the branching ratios to higgsinos, regrouping all the decays to the lightest (next-to-lightest) neutralino and chargino states (h {χ 0 1 ,χ 0 2 ,χ ± 1 }), as well as to ττ -this will be discussed in Sect. 3.2 in connection with future tests of our scenario.
Regarding the production cross sections of the heavy CPeven Higgs states at the LHC, they are essentially driven by the doublet components of these states. They are shared in roughly equal proportions by the H 2,3 states, as a consequence of the ∼50 % mixing. The cross section in ggf is suppressed as tan −2 β (in agreement with the tan −1 β suppression of the H D -coupling to tops) while that in bbh varies as tan 2 β (in accordance with the tan β enhancement of the H D -coupling to bottoms): summing over both states, the production cross section in ggf at 8 TeV follows the approximate rule σ ggf (0.07 fb) tan 2 β; at 13 TeV, σ ggf 13TeV [H 2 + H 3 ] (600 fb) tan −2 β and σ bbh 13TeV [H 2 + H 3 ] (0.4 fb) tan 2 β. At low tan β ∼ 4-5, the ggf channel thus dominates, while the bbh is more efficient at large tan β ∼ 10-15. Their sum is maximal at large tan β: ∼17 fb at 8 TeV and ∼95 fb at 13 TeV, and it is minimal for tan β ∼ 6: ∼6 fb at 8 TeV and ∼30 fb at 13 TeV. The enhancement factor between 8 and 13 TeV ranges from ∼4.75 at tan β = 4 to 5.65 at tan β = 15. Considering that the 8 TeV data did not show any significant diphoton excess at ∼750 GeV, one would prefer this enhancement factor to be as large as possible, so that it avoids tensions with the limits from Run-I. An enhanced associated production with b-quarks, which occurs at large tan β is thus slightly preferred since the ratio between the 13 TeV and 8 TeV production cross section is the largest for bb initial states. The production cross sections of A 2 at 13 TeV are also given. They follow coarsely the same patterns as their analogs for H 2 /H 3 , with a larger ggf though.
These production cross sections and branching ratios allow us to derive the relevant cross section for the pp → H 2 (H 3 ) → 2(A 1 → γ γ ) process (we assume BR[A 1 → γ γ ] 0.99): this quantity is documented for our points in the last two lines of Table 2, at 8 and 13 TeV. Depending on the characteristics of our points, σ incl 8TeV ∼1-2 fb and σ incl 13TeV ∼5-13 fb, which is the relevant order of magnitude for an interpretation of the diphoton excess. These figures shall be analyzed with further detail in the following section.

Analysing the diphoton signal with current data
We consider resonant H 2 (H 3 ) production, where X stands for the rest of the event. Due to the large boost of A 1 , the two photons of the A 1 decay will be very collimated and thus the opening angle between both photons in the electronic calorimeter will be well below the angular resolution of electromagnetic calorimeters [89,90]. Therefore, our final state resembles a resonant diphoton final state.
The ATLAS conference note [91] studied the signature of a 125 GeV Higgs boson decaying to four photons at √ s = 7 TeV with the full data set. The search estimates the efficiency of photon-pair identification as a single photon at about 85-95 % for photons with p T ≈ 100 GeV and mass m A 1 = 100-200 MeV. The efficiency heavily depends on the mass of the heavy resonance, as can be seen in Fig. 3, which shows the η := |η(γ 1 ) − η(γ 2 )| distribution of the two photons from the A 1 → γ γ decay where η denotes the pseudorapidity. The dark (blue) curve shows the result for a Higgs boson mass of 125 GeV, while the light (beige) curve for 750 GeV. As expected the opening angle of the photons from the 750 GeV resonance is much smaller compared to the 125 GeV case. However, it is difficult to determine the exact efficiency without performing a full detector simulation. Hence, we will choose = 90 % in the remainder of the paper. In any case, our conclusions will not considerably change if we assume a higher efficiency as the one from the ATLAS study [91]. This means that about 80 % of all four-photon events within the fiducial region will be classified as diphoton events. This choice is further supported by the analysis in Ref. [92].   Table 3. The resulting signal efficiency varies between 20 and 60 % depending on the signal region, the experiment and the center-of-mass energy.
Using the above setup we calculate the expected number of events in the signal regions centered at 750 GeV for each parameter point P1-P9 in four 8 TeV searches and in two 13 TeV searches. The results are collected in Table 4. For reference, in column two and three we provide the observed number of events above the SM background ("sig.") and the observed S95 exclusion limits calculated using CheckMATE [99,102]. Our benchmark points offer a range of cross sections for the desired signal, from 5.8 to ; see also the discussion in [5]. The remaining points fulfill all constraints. In a model independent χ 2 fit we estimate that the best-fit cross section at 13 TeV using ATLAS and CMS results is 8.3 fb [5], with points P2 and P7 being closest in value, cf. Table 2. A similar analysis for the 8 TeV data yields 0.5 fb, while a combination of all the available data gives a range of cross sections 4.9-5.7 fb at 13 TeV. The exact result depends on the details of the production mechanism, but our benchmark points cover well the desired range. The best fit to all data is provided by point P5.
In Fig. 4, we show the diphoton invariant mass distribution of our diphoton signal for two different bin sizes. We consider benchmark point P6 for illustration. The distribution with the large bin size of 40 GeV corresponds to the experimental bin size of the ATLAS study [1], as shown in the left panel. The experimental photon energy resolution of about 5-10 % would allow for a higher precision [103] but due to the small statistical sample, both experiments have to choose a rather large bin size. One can clearly see that our benchmark point with two scalars cannot be distinguished from a wide resonance with the current data. For comparison we have included into this plot the original data from ATLAS after subtracting the expected background. One can see that the events predicted for our P6 benchmark provide a good reproduction of the experimental shape. We also display in the right panel of Fig. 4 the invariant mass distribution with a 5 GeV binning. While currently the experimental resolution in m γ γ exceeds 10 GeV, one can speculate that further improvements during the current LHC run will be made. With the accuracy of ∼5 GeV and an increased luminosity, the broad excess, provided it is real, might be resolved as two narrow resonances [8].
Our estimate of the number of signal events in Table 4 has a theoretical uncertainty. The choice of the parton distribution function, missing higher-order calculations and the details of the parton shower as well as the fragmentation induce an uncertainty on the signal rate. However, our finalstate configuration is relatively simple and thus the details of the MC tuning will not significantly alter our results. The size of the uncertainty from the PDFs can be sizable. The variation between the different PDFs changes the hadronic cross section. In addition, the scale dependence of the production cross section on the renormalization and factorization scale affects the signal rate. We estimate the sum of these effects to be of the order of 10 % [104]. Furthermore, the normalization can heavily depend on the value of the bottom Yukawa coupling [104]. Finally, we did not model the detector response of identifying a photon from the pseudoscalar decay into the diphoton state but rather choose a flat efficiency factor of = 90 %. Here, one can assume a conservative uncertainty of about 20 % on the signal rate due to the uncertainty in the photon identification. We conclude that the total uncertainty is of the order of O(20) %, plus an additional uncertainty from the definition of the bottom quark mass.

Future directions
So far we have assumed that our signal in Eq. (17) mimics the diphoton signal since the two collimated photons of the light-pseudoscalar decay are indistinguishable from an isolated photon. However, if the four-photon final state was discriminated from the diphoton signature, it would be a strong hint at our scenario. References [105][106][107] considered photon jets (two or more collimated photons) at hadron colliders. In particular, Ref. [107] discussed the possibility of photon conversion into e + e − pairs and its discriminating power between photon jets and isolated photons. For a photon jet, the prob- ability of photon conversion is higher than for a single photon, and Ref. [107] showed that already several tens of events are sufficient to discriminate between both hypotheses and a few hundred events allow for a 5σ discrimination assuming prompt photons. However, their conclusions assume a pseudoscalar mass of 1 GeV and the results are very sensitive to this parameter. For long-lived pseudoscalars, the discriminating power is reduced since photon conversion cannot start before the pseudoscalar decay. As a consequence, the discriminating power becomes worse for increasing lifetimes.
Apart from the diphoton signal, which is the main motivation of the current study, the NMSSM parameter points discussed here also have additional distinctive features closely related to the diphoton signal. We shortly discuss these collider signatures.
As discussed in the previous section, the light pseudoscalar, A 1 , has a small branching fraction of 1 % for decays to electron pairs. Because of its short life-time it would typically decay promptly to a highly collimated e + e − pair, so-called "electron jet". Such electron jets, prompt or displaced, were searched for by the LHC experiments as they can appear in many different models of new physics. In our case, two signatures can appear: two highp T electron jets or one electron jet and an energetic photon. Note that even though we have suppressed the branching ratio of the SM-like Higgs boson to pseudoscalars such a signal is also possible apart from the decays of the heavy Higgs states. For the 125 GeV Higgs an associated production of pp → hW (→ ν) (where h denotes the SM-like Higgs boson) was studied by ATLAS [108]. The obtained limit is weak and together with the already mentioned suppression of H 1 → A 1 A 1 this is an unlikely discovery channel. The searches for the direct production of the scalar decaying to two electron jets could provide further constraints, but the limits have been obtained only for the light SM-like Higgs boson [109,110]. Nevertheless this signature might become interesting in the current run once the diphoton signal is firmly confirmed. Corresponding studies directly applicable to the heavy Higgs particles have also been performed [111,112]. It is interesting to note that for a hypothetical scenario with a 1 TeV scalar resonance, the sensitivity of the CMS search [111] is in the fb range. While the discussed 8 TeV searches lack the sensitivity to constrain our scenario now, they clearly offer interesting prospects for observing electron decay modes of A 1 (possibly accompanied by the photon jet from the opposite decay chain) at the increased center-of-mass energy and the high luminosity run of the LHC.
Our scenario can also be probed via the "classic signature" for additional heavy neutral Higgs bosons, pp → → τ + τ − , where the limits are set in the m -tan β space. Within the MSSM, assuming the additional Higgs bosons at a mass around ∼750 GeV, the (expected) limits on tan β are around ∼35 based on Run I data [113][114][115] (see also [116]). In our NMSSM scenario we have three Higgs bosons with a mass around 750 GeV contributing to this search channel, H 2 , H 3 , and A 2 , where the overall number of τ + τ − events is roughly 25 % lower than in the MSSM, mainly due to the decay of H 2,3 → A 1 A 1 . Consequently, a similar, but slightly higher limit on tan β can be set in our NMSSM scenario. With increasing luminosity this limit could roughly improve to tan β ∼ 5-10 at the LHC after collecting 300-3000/fb of integrated luminosity (see also [117]). Therefore, the proposed scenario could eventually lead to an observable signal in the τ + τ − searches for heavy Higgs bosons at the LHC, depending on the details of the scenario (value of tan β, masses of electroweak particles etc.).
Another prediction that arises for parameter points considered in this study are light higgsinos. With the masses of 100-300 GeV they are well within the kinematic reach of the LHC. However, the small mass differences, O(10 GeV), within the light higgsino sector hinder their observation at the LHC. If all the non-higgsino SUSY particles are sufficiently far in mass (points P1-P4, P8, P9), the decay of the second neutralino,χ 0 2 proceeds almost exclusively via the light pseudoscalar A 1 . With the following significant branching ratio to a soft γ γ pair the observation in the soft di-and trilepton searches [118,119] becomes practically impossible. The radiative production at a high-energy e + e − collider remains a valid possibility though [120,121].
Finally, light smuons are required in order to obtain phenomenologically viable muon anomalous magnetic moment and to counteract the effects of a very light pseudoscalar. For our parameter points we fix slepton masses at 300-400 GeV. While this is close to the existing simplified model limits, see e.g. [122], in our case due to the significant branching ratio BR(˜ L →χ ± 1 ν) 50 % these constraints are significantly relaxed. Nevertheless, if the slepton and higgsino spectra are favorable, the observation in the current LHC run is plausible.

Conclusions
We have proposed an NMSSM scenario that can explain the excess in the diphoton spectrum at 750 GeV recently observed by ATLAS and CMS. In our scenario the heavy neutral (and charged) Higgs bosons have a mass around ∼750 GeV, while one light CP-odd Higgs boson has a mass around the mass of the pion, ∼135 MeV. The 750 GeV excess is generated by the production of the heavy neutral CP-even Higgs bosons, which subsequently decay to two light pseudoscalars. Each of these pseudoscalars then decays, mainly via the mixing with the π 0 , to a collimated photon pair that appears as a single photon in the electromagnetic calorimeter. The mass gap between heavy Higgses of O (20) GeV mimics a large width of the 750 GeV peak. Furthermore, the production of the heavy neutral CP-even Higgs bosons may contain a large component of bb initial state, thus ameliorating a possible tension with 8 TeV data compared to other production modes. The main virtue of our scenario is that all necessary properties to fit the signal can be united in a phenomenologically realistic way within as theoretically simple a model as the NMSSM, without e.g. requiring additional ad hoc matter.
We derived the NMSSM parameter space in which this scenario can be realized. It is characterized by a heavy Higgs boson mass scale, M A , around 750 GeV. The Yukawa-like couplings λ and κ are found to satisfy 0.4 tan β 1+tan 2 β λ 2 √ 2 tan β √ 1+18 tan 2 β+tan 4 β and κ λ 2 sin 2β . The μ parameter is given by μ ∼ M A sin 2β, or 2 κ λ μ 750 GeV. We furthermore find 5 tan β 15, A κ O(0.1) GeV, and A λ v. Due to this choice of parameters the two neutral heavy CPeven Higgs bosons are strongly mixed doublet/singlet states and the light CP-odd Higgs boson can have a mass around m π . The light CP-even Higgs boson with SM-like properties can have a mass around ∼125 GeV, mainly by choosing the scalar top parameters accordingly. The parameter choice furthermore forbids a large decay rate of the SM-like Higgs boson to the two light CP-odd states, which would be in contradiction with the LHC measurements.
In order to validate our scenario we have chosen nine benchmark points, all satisfying the above constraints, but with a strong variation within the allowed intervals. Using state-of-the-art tools, including higher-order corrections, these points have been analyzed to reproduce the observed "excess" in the diphoton search at the LHC Run-II, including detector simulation and efficiencies. We have furthermore checked explicitly that these points fulfill all other experimental constraints. These include LHC Higgs (and SUSY) searches, Higgs boson rate measurements, as well as flavor observables and electroweak precision data. We have shown explicitly that the two collimated photon pairs would be seen as a single photon each, applying the same settings as in the ATLAS/CMS analyses.
Finally, we have analyzed how our scenario can be probed in the upcoming continued LHC Run-II. Possibly striking features are the absence of any other relevant decay mode, such as the decay to massive gauge bosons, as well as an increased rate of photon conversion to electron jets with respect to the "simple" diphoton decay mode, or the distinction of a photon pair from a single photon. Furthermore, the heavy neutral Higgs bosons should be visible in the conventional τ + τ − searches at high luminosity. Other characteristic features of our scenario are relatively light higgsinos and possibly sleptons that can be probed at the LHC Run-II. Using these characteristics, our scenario should be distinguishable from most other physics scenarios that have been proposed to explain the LHC diphoton "excess".