# A 750 GeV diphoton signal from a very light pseudoscalar in the NMSSM

Ulrich Ellwanger, Cyril Hugonie

► To cite this version:

Ulrich Ellwanger, Cyril Hugonie. A 750 GeV diphoton signal from a very light pseudoscalar in the NMSSM. Journal of High Energy Physics, 2016, 05, pp.114. 10.1007/JHEP05(2016)114. hal-01996238

HAL Id: hal-01996238

https://hal.science/hal-01996238

Submitted on 7 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RECEIVED: March 24, 2016

REVISED: May 2, 2016

ACCEPTED: May 11, 2016

PUBLISHED: May 19, 2016

# A 750 GeV diphoton signal from a very light pseudoscalar in the NMSSM

Ulrich Ellwanger$^{a,b}$ and Cyril Hugonie$^c$

$^a$Laboratoire de Physique Théorique, UMR 8627, CNRS, Université de Paris-Sud, Univ. Paris-Saclay, Bât. 210, 91405 Orsay, France

$^b$School of Physics and Astronomy, University of Southampton, Highfield, Southampton SO17 1BJ, U.K.

$^c$LUPM, UMR 5299, CNRS, Université de Montpellier, Place Eugène Bataillon, 34095 Montpellier, France

E-mail: ellwanger@th.u-psud.fr, cyril.hugonie@umontpellier.fr

**ABSTRACT:** The excess of events in the diphoton final state near 750 GeV observed by ATLAS and CMS can be explained within the NMSSM near the *R*-symmetry limit. Both scalars beyond the Standard Model Higgs boson have masses near 750 GeV, mix strongly, and share sizeable production cross sections in association with b-quarks as well as branching fractions into a pair of very light pseudoscalars. Pseudoscalars with a mass of ~210 MeV decay into collimated diphotons, whereas pseudoscalars with a mass of ~500–550 MeV can decay either into collimated diphotons or into three $\pi^0$ resulting in collimated photon jets. Various such scenarios are discussed; the dominant constraints on the latter scenario originate from bounds on radiative $\Upsilon$ decays, but they allow for a signal cross section up to 6.7 fb times the acceptance for collimated multiphotons to pass as a single photon.

**KEYWORDS:** Supersymmetry Phenomenology

ARXIV ePRINT: 1602.03344
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Parameter regions with diphoton-like events at 750 GeV in the NMSSM</td><td>3</td></tr><tr><td>3</td><td>Viable scenarios with a light NMSSM pseudoscalar</td><td>7</td></tr><tr><td>3.1</td><td>$M_{A_1}$ near 210 MeV</td><td>8</td></tr><tr><td>3.2</td><td>$M_{A_1}$ at 510–550 MeV</td><td>11</td></tr><tr><td>3.3</td><td>The spread in rapidity of multiphotons</td><td>15</td></tr><tr><td>4</td><td>Summary and conclusions</td><td>16</td></tr></table>

## 1 Introduction

In December 2015 the ATLAS and CMS collaborations have reported excesses in the search for resonances decaying into pairs of photons for diphoton invariant masses around 750 GeV [1, 2]. In ATLAS, excesses appeared in the two $M_{γγ}$ bins 710–750 GeV (14 events vs. 6.3 expected) and 750–790 GeV (9 events vs. 5.0 expected), with a local significance of $3.9 \sigma$ (assuming a large width of ~45 GeV; $3.6 \sigma$ in the narrow width approximation). In CMS, excesses appear in the $M_{γγ}$ bin 750–770 GeV for photons in the EBEB category (5 events vs. 1.9 expected) and EBEE category (6 events vs. 3.5 expected), but less in the bin 730–750 GeV (4 events vs. 2.1 expected for photons in the EBEB category, 1 event vs. 4.0 expected for photons in the EBEE category, considered as less sensitive). The local significance of the excesses is $2.6 \sigma$ for CMS in the narrow width approximation.

The global significances of the signals of $O(2-3\sigma)$ are not overwhelming and compatible with statistical fluctuations. Still, the fact that the region of invariant diphoton masses is very similar for ATLAS and CMS has stirred quite some excitement resulting in a huge number of possible explanations. (The number of proposed models exceeds the number of observed signal events.)

Fits to the combined data should, in principle, also consider the informations from diphoton searches at 8 TeV [3, 4] where a mild excess was observed by CMS. However, the extrapolation of signal cross sections from 8 to 13 TeV depends on the assumed production mechanism [5–9]. Assuming the production of a resonance around 750 GeV by gluon fusion (ggF), combined fits to the signal cross sections at 13 TeV are in the range 2–10 fb [5–7, 9], with slightly better fits and a larger signal cross section assuming a larger width of 30–45 GeV [5, 6, 9].

It is notoriously difficult to construct a consistent model for such a resonance “X”: its production channel in proton proton collisions is typically assumed to be ggF through loops of colored particles. If these are the quarks of the Standard Model (SM), X would decay into them leaving little branching fraction for X-decays into $γγ$, which has to be generated by loop diagrams as well.
Accordingly simple two Higgs doublet (or MSSM) extensions of the Standard Model, which could contain a resonance X near 750 GeV [10–23], require additional scalars or vector-like fermions whose loops generate the coupling of X to gluons and/or $\gamma\gamma$ (unless R-parity is broken [24, 25]). Large Yukawa couplings are required for a sufficiently large cross section, which risk to generate new hierarchy problems/Landau singularities (unless compositeness is invoked). Also in the Next-to-Minimal supersymmetric extension of the Standard Model (NMSSM) it has been argued [26–28] that additional vector-like quark superfields have to be introduced. In [29] a two-step decay cascade involving the two pseudoscalars of the NMSSM with masses of about 750 GeV and 850 GeV has been proposed which requires, however, to tune the corresponding mixing angle close to 0.

A different approach towards an explanation of the diphoton events is to consider that a single photon in the detector can represent a collimated bunch of photons (typically two of them) which originate from a single very light state, for instance a light pseudoscalar $A$ [5, 30–37]. Then the observed processes correspond to an initial resonance $X$ decaying into a pair $AA$, where $M_A$ must be well below 1 GeV for the resulting photons to be sufficiently collimated (see below). This scenario opens the possibility to explain the diphoton events in different models which can accomodate resonances $X$ and a light pseudoscalar $A$. In this paper we show that the simple $Z_3$-invariant NMSSM belongs to this class of models. (This has also been observed in [38].)

In the NMSSM (see [39, 40] for reviews), two CP-even Higgs states beyond the Standard Model-like Higgs (subsequently denoted as $H_{SM}$) can play the role of a resonance $X$. In terms of weak eigenstates, a singlet-like state $S$ can have a large coupling to a pair of mostly singlet-like pseudoscalars $A_1$, originating from a cubic singlet self coupling $\kappa$ in the superpotential (see below). However, a coupling to quarks or gluons inside protons has to be induced by a mixing of $S$ with one of the two SU(2) doublet-like Higgs states. If this state is $H_{SM}$, the mixing reduces the couplings of $H_{SM}$ to SM particles (notably $W^\pm$ and $Z$) and is severely constrained [9] by the measured signal rates by ATLAS and CMS [41]. An alternative is that $S$ mixes strongly with the other “MSSM”-like CP-even state $H$. Then the physical eigenstates — preferably both of them with masses near 750 GeV — can profit from an enhancement of the couplings of $H$ to $b$-quarks by $\tan\beta$, leading to sufficiently large signal cross sections into the $A_1A_1$ (and hence diphoton) final state via associated production with $b$-quarks. Given the diphoton mass resolution of the detectors and the slightly preferred large width of the excess it is clear that two (narrow) CP-even states near 750 GeV, mixtures of $H$ and $S$, can also provide a good fit to the data. (A similar scenario has been discussed in [42].) For one of the benchmark points presented below (BP1) the signal originates, however, from one CP-even state only, the other one being significantly heavier.

A light pseudoscalar can appear in the NMSSM in the form of a pseudo-Goldstone boson (PGB). A priori two global symmetries can lead to such PGBs: first, a Peccei-Quinn symmetry emerges in the limit $\kappa \to 0$ [40, 43–45]. However, $\kappa \neq 0$ is required for the couplings of the heavy Higgs states to $A_1A_1$. Second, the scalar potential of the NMSSM is invariant under an $R$-symmetry [40, 46–48] if the soft supersymmetry breaking trilinear couplings $A_\lambda$ and $A_\kappa$ vanish, leading to a PGB due to its spontaneous breakdown
by the phenomenologically required vacuum expectation values. We find indeed, that the interesting part of the parameter space of the NMSSM corresponds to small values of $A_\lambda$ and $A_\kappa$. However, since the R-symmetry is broken by radiative corrections to the scalar potential involving the necessarily non-vanishing gaugino masses and trilinear couplings $A_t$ and $A_b$, it helps only partially to explain a very light pseudoscalar $A_1$. Still, it represents a “go-theorem” showing that a standard supersymmetric extension of the SM — without additional vector-like quarks and/or leptons — could explain the observed diphoton excess.

Different assumptions on the mass of $A_1$ can be made. For one set of scenarios we assume $M_{A_1} \sim 210$ MeV, just below the $2\mu$ threshold. These scenarios lead to visibly displaced vertices from the $A_1 \to \gamma\gamma$ decays. For a large value of the NMSSM trilinear coupling $\kappa \sim 1.65$, the signal can originate from a single Higgs state near 750 GeV. For smaller values of $\kappa$, the signal can originate from two Higgs states with masses near 750 GeV. For another set of scenarios we assume $M_{A_1} \sim 510-550$ MeV, not far from the $\eta$ mass. For $M_{A_1}$ near 550 MeV, $A_1$ mixes with the $\eta$ meson and inherits its decays into $\gamma\gamma$ and $3\pi^0$; the latter lead to photon-jets. The average separation in rapidity of the diphotons and the two leading photons from $3\pi^0$ will be studied. For $M_{A_1}$ near 510 MeV, constraints from searches for radiative $\Upsilon(1S)$ decays into $\gamma + \eta$ by CLEO [49] are alleviated, but estimates of the $A_1$ decay widths are more uncertain. But in both cases the $A_1$ life time is short enough avoiding macroscopically displaced vertices, and two Higgs states near 750 GeV can generate a signal.

In the next section we describe with the help of analytic approximations to the mass matrices (including only the dominant radiative corrections) which region in the parameter space of the NMSSM can generate the diphoton events. In section 3 we discuss various constraints from low energy physics on light pseudoscalars, and discuss separately the different scenarios. Benchmark points are presented with the help of the public Fortran code NMSSMTools [50, 51]. For the different $A_1$ masses we study the average separation in rapidity of the diphotons and the two leading photons from $3\pi^0$, which allows to estimate the corresponding acceptances. In the final section 4 we summarize and discuss possible alternative signatures, which could help to distinguish different scenarios if the excess survives the next runs of the LHC.

## 2 Parameter regions with diphoton-like events at 750 GeV in the NMSSM

We consider the CP-conserving $Z_3$-invariant NMSSM. The superpotential of the Higgs sector reads in terms of hatted superfields

$$W_{\text{Higgs}} = \lambda \hat{S} \hat{H}_u \cdot \hat{H}_d + \frac{\kappa^3}{3} \hat{S}^3. \qquad (2.1)$$

Once the real component of the singlet superfield $\hat{S}$ develops a vacuum expectation value (vev) $s$, the first term in $W_{\text{Higgs}}$ generates an effective $\mu$ term

$$\mu = \lambda s. \qquad (2.2)$$
The soft SUSY-breaking terms consist of mass terms for the gaugino, Higgs and sfermion fields

$$ -\mathcal{L}_{\frac{1}{2}} = \frac{1}{2} \left[ M_1 \tilde{B} \tilde{B} + M_2 \sum_{a=1}^{3} \tilde{W}^a \tilde{W}_a + M_3 \sum_{a=1}^{8} \tilde{G}^a \tilde{G}_a \right] + \text{h.c.}, \\ -\mathcal{L}_0 = m_{H_u}^2 |H_u|^2 + m_{H_d}^2 |H_d|^2 + m_S^2 |S|^2 + m_Q^2 |Q|^2 + m_T^2 |T_R|^2 \\ + m_B^2 |B_R^2| + m_L^2 |L|^2 + m_\tau^2 |\tau_R^2|, \qquad (2.3) $$

as well as trilinear interactions between the sfermion and the Higgs fields, including the singlet field

$$ -\mathcal{L}_{\text{tril}} = \left( h_t A_t Q \cdot H_u T_R^c + h_b A_b H_d \cdot Q B_R^c + h_\tau A_\tau H_d \cdot L \tau_R^c \right. \\ \left. + \lambda A_\lambda H_u \cdot H_d S + \frac{1}{3} \kappa A_\kappa S^3 \right) + \text{h.c.} . \qquad (2.4) $$

The tree level scalar potential can be found in [40], from which the 3 × 3 mass matrices in the CP-even and CP-odd sectors can be obtained. Once the soft Higgs masses are expressed in terms of $M_Z$, $\tan \beta$ and $s$ using the minimization equations of the potential, the mass matrices depend on the six parameters

$$ \lambda, \quad \kappa, \quad \tan \beta = \frac{v_u}{v_d}, \quad \mu, \quad A_\lambda \quad \text{and} \quad A_\kappa. \qquad (2.5) $$

Initially, the CP-even mass matrix $M_S^2$ is obtained in the basis of the real components ($H_{d,r}$, $H_{u,r}$, $S_r$) of the complex scalars ($H_d$, $H_u$, $S$) after expanding around the vevs $v_d$, $v_u$ and $s$. It is convenient, however, to rotate $M_S^2$ by an angle $\beta$ in the doublet sector into $M_S'^2$ in the basis $H'_{SM}$, $H'$, $S_r$:

$$ M_S'^2 = R(\beta) M_S^2 R^T(\beta), \qquad R(\beta) = \begin{pmatrix} \cos \beta & \sin \beta & 0 \\ \sin \beta & -\cos \beta & 0 \\ 0 & 0 & 1 \end{pmatrix}. \qquad (2.6) $$

The advantage of this basis is that only the component $H'_{SM}$ of the Higgs doublets acquires a vev $v$ and that, for typical parameter choices, it is nearly diagonal: $H'_{SM}$ has SM-like couplings to fermions and electroweak gauge bosons, the heavy doublet field $H'$ is the CP-even partner of the MSSM-like CP-odd state $A_{MSSM}$, while $S_r$ remains a pure singlet. The mass matrix $M_S'^2$ in the basis ($H'_{SM}$, $H'$, $S_r$) has the elements

$$ \begin{aligned} M_{S,11}'^2 &= M_Z^2 \cos^2 2\beta + \lambda^2 v^2 \sin^2 2\beta + \sin^2 \beta \Delta_{\text{rad}}, \\ M_{S,12}'^2 &= \sin 2\beta \left( \cos 2\beta (M_Z^2 - \lambda^2 v^2) - \frac{1}{2} \Delta_{\text{rad}} \right), \\ M_{S,13}'^2 &= \lambda v (2\mu - (A_\lambda + 2\kappa s) \sin 2\beta), \\ M_{S,22}'^2 &= M_A^2 + (M_Z^2 - \lambda^2 v^2) \sin^2 2\beta + \cos^2 \beta \Delta_{\text{rad}}, \\ M_{S,23}'^2 &= \lambda v (A_\lambda + 2\kappa s) \cos 2\beta, \\ M_{S,33}'^2 &= \lambda A_\lambda \frac{v^2}{2s} \sin 2\beta + \kappa s (A_\kappa + 4\kappa s), \end{aligned} \qquad (2.7) $$
where $v^2 = 2M_Z^2/(g_1^2 + g_2^2) \sim (174 \text{ GeV})^2$ and

$$M_A^2 = \frac{2\mu}{\sin 2\beta} (A_\lambda + \kappa s) \qquad (2.8)$$

is the mass squared of the MSSM-like CP-odd state $A_{\text{MSSM}}$. $\Delta_{\text{rad}}$ denotes the dominant radiative corrections due to top/stop loops,

$$\Delta_{\text{rad}} = \frac{3m_t^4}{4\pi^2 v^2} \left( \ln \left( \frac{m_{ST}^2}{m_t^2} \right) + \frac{X_t^2}{m_{ST}^2} \left( 1 - \frac{X_t^2}{12m_{ST}^2} \right) \right) \qquad (2.9)$$

where $m_{ST}^2 = m_Q m_T$ and $X_t = A_t - \mu / \tan \beta$.

As discussed in the introduction, we intend to describe the diphoton signal at $\sim 750 \text{ GeV}$ by a mixture of the two states $H'$ and $S_r$. Then, both diagonal matrix elements $M_{S,22}'$ and $M_{S,33}'$ should have values close to $(750 \text{ GeV})^2$. Furthermore we will be interested in the R-symmetry limit $A_\lambda, A_\kappa \to 0$. This implies the relations (for $\tan^2 \beta \gg 1$)

$$M_{S,22}' \sim M_A^2 \sim \frac{2\mu\kappa s}{\sin 2\beta} \sim \frac{\kappa}{\lambda}\mu^2 \tan \beta \sim (750 \text{ GeV})^2 \qquad (2.10)$$

and

$$M_{S,33}' \sim (2\kappa s)^2 \equiv 4 \left(\frac{\kappa}{\lambda}\right)^2 \mu^2 \sim (750 \text{ GeV})^2. \qquad (2.11)$$

The matrix element inducing $H' - S_r$ mixing is given by

$$M_{S,23}' \sim 2\kappa v \mu, \qquad (2.12)$$

and the matrix element inducing $H'_{\text{SM}} - S_r$ mixing by

$$M_{S,13}' \sim 2\lambda v \mu. \qquad (2.13)$$

Next we turn to the CP-odd sector. The $3 \times 3$ CP-odd mass matrix contains always a Goldstone boson which will be eaten by the Z boson. The remaining CP-odd states are a singlet $A_S$, and the “MSSM”-like SU(2)-doublet $A_{\text{MSSM}}$. In the basis ($A_{\text{MSSM}}, A_S$), in the R-symmetry limit $A_\lambda, A_\kappa \to 0$, the CP-odd mass matrix is given by

$$M_A^2 = \frac{2\kappa\mu}{\sin 2\beta} \begin{pmatrix} s & -v \sin 2\beta \\ -v \sin 2\beta & \frac{v^2}{s} \sin^2 2\beta \end{pmatrix}. \qquad (2.14)$$

Obviously $M_A^2$ has a vanishing eigenvalue $M_{A_1} = 0$, and is diagonalised by an angle $\alpha$ with (for $\tan^2 \beta \gg 1$)

$$\sin \alpha \approx \frac{2v}{s \tan \beta}. \qquad (2.15)$$

An important quantity will be the (reduced) coupling $X_d$ of $A_1$ to down quarks and leptons, which is obtained through the mixing of $A_S$ with $A_{\text{MSSM}}$. Since the reduced coupling of the MSSM-like state $A_{\text{MSSM}}$ is given by $\tan \beta$, one obtains

$$X_d \sim \sin \alpha \tan \beta \sim \frac{2v}{s} \equiv \frac{2\lambda v}{\mu}. \qquad (2.16)$$
Radiative corrections to the tree level potential and hence to the CP-odd mass matrix include terms proportional to the electroweak gaugino masses $M_1$ and $M_2$, and terms proportional to the soft SUSY breaking trilinear couplings $A_t$ and $A_b$. These corrections break the R-symmetry present for $A_\lambda$, $A_\kappa \to 0$, which is expected since $A_\lambda$, $A_\kappa = 0$ is not invariant under scale transformations. Hence, depending on the scale where $A_\lambda$, $A_\kappa = 0$ is assumed to hold, $A_1$ is a pseudo-Goldstone boson with a mass of typically a few GeV. For $A_\kappa$ small, but $\neq 0$ one can obtain $M_{A_1} \sim 210$ MeV or $M_{A_1} \sim 510-550$ MeV as it will be assumed in the next section.

Finally we note that, for the parameter region considered below, the dominant contribution to the coupling of $A_1$ to scalars originates from the quartic coupling $\sim \kappa^2|S|^4 \to 2\kappa^2 S_r^2 A_1^2$. After shifting $S_r$ by its vev $s$ one obtains

$$g_{SA_1A_1} \sim \sqrt{2\kappa^2}s. \qquad (2.17)$$

Next we observe that eqs. (2.11) and (2.16) allow to express $\kappa$ in terms of $X_d$: from (2.11) one finds

$$750 \text{ GeV} \sim 2\kappa s = 2\kappa \frac{\mu}{\lambda} = \frac{4\kappa v}{X_d} \qquad (2.18)$$

where (2.16) was used in the last step. Inserting $v \sim 174$ GeV one obtains

$$\kappa \sim 1.1 X_d. \qquad (2.19)$$

In the next section, for the scenarios with $M_{A_1} \sim 510-550$ MeV, we will obtain upper bounds on $X_d$ from upper bounds for the $BR(\Upsilon(1S) \to \gamma\eta)$ from CLEO [49]. These will thus imply upper bounds on $\kappa$ according to (2.19). On the other hand a large signal rate, generated by a mixture of the states $H'$ and $S$ decaying into $A_1A_1$, requires $g_{SA_1A_1}$ to be as large as possible. Accordingly $X_d$ and $\kappa$ should saturate corresponding upper bounds.

If the 750 GeV signal is generated by a superposition of signals of two nearby physical states formed by the $H' - S_r$ system, their mass splitting should not be too large, preferably of $\mathcal{O}(20 \text{ GeV})$. Then the matrix element $M_{S,23}^2$ given in (2.12) should be as small as possible. With $\kappa$ already determined, this implies $\mu$ as small as possible, preferably close to the lower bound $\sim 100$ GeV from the LEP lower bound on higgsino-like charginos. Then (2.16) requires that $\lambda$ is relatively small. (Simultaneously, this avoids a strong push-down effect on the mass of the SM-like Higgs boson from $H_{SM} - S_r$ mixing, which is induced by the matrix element $M_{S,13}^2$ given in (2.13).) Finally the condition (2.10) on $M_A^2$ fixes $\tan \beta \approx 15$.

The remaining NMSSM parameters in (2.5) are $A_\lambda$ and $A_\kappa$. Both R-symmetry breaking parameters have an impact on the mass of the pseudo-Goldstone boson $A_1$. We find that one can choose small values of $A_\lambda$ and $A_\kappa$ such that $M_{A_1}$ assumes the desired value; due to radiative corrections to the scalar potential the precise value of $A_\kappa$ depends on the other R-symmetry breaking parameters $M_1$, $M_2$, $A_t$ and $A_b$. Herewith all NMSSM parameters are nearly uniquely determined.
### 3 Viable scenarios with a light NMSSM pseudoscalar

As discussed in the introduction we will study scenarios with different values of the mass of a light pseudoscalar, denoted subsequently by $M_{A_1}$. Constraints on such a light NMSSM pseudoscalar with a mass below ~1 GeV have been discussed previously in [52–59]. Strong constraints originate from the mediation of FCNCs. Assuming minimal flavour violation, flavour violating couplings of $A_1$ still originate from SUSY loops involving stops, sbottoms and charginos and depend on the corresponding masses and trilinear couplings like $A_t$. These contribute notably to B-physics observables like $B_s \to \mu^+\mu^-$, $\Delta M_d$ and $\Delta M_s$. We have implemented the computation of these and many more B-physics observables and some K-physics observables in the code NMSSMTools [50, 51] following the update in [59] and checked that, for the scenario presented here, the constraints are satisfied due to the mostly singlet-like nature of $A_1$ and the relatively heavy SUSY spectrum.

For $M_{A_1}$ near 210 MeV, additional strong constraints originate from rare flavour changing processes $K^\pm \to \pi^\pm e^+ e^-$. (In [53] it has been argued that the corresponding constraints exclude scenarios with $M_{A_1} \lesssim 210$ MeV, where the branching fraction of $A_1$ into $e^+e^-$ is sizeable.) We have verified the assertion in [38] that, for suitable choices of soft SUSY breaking parameters, the coupling $C_A$ responsible for these processes (see [52]) can be arbitrarily small.<fnref n="1" /> Light pseudoscalars have been searched for in radiative $\Upsilon(1S)$ decays by CLEO in [60]; these are also verified by NMSSMTools 4.9.0 and satisfied by the benchmark points given below.

Due to the mostly singlet-like nature of $A_1$, its contributions to the muon anomalous magnetic moment are negligibly small. However, for $\tan\beta \sim 15$ and assuming relatively light slepton masses of 300 GeV, the scenarios below can reduce the discrepancy between the measured value and the Standard Model to an acceptable $2\sigma$ level.

Further constraints stem from possible $A_1$ production in $Z$ and $H_{SM}$ decays. The relevance of bounds on light pseudoscalars (or axion-like particles) from searches for $Z \to \gamma\gamma$ at LEP (where a photon can correspond to a bunch of collimated photons) has been investigated in [58]. These bounds constrain the loop-induced coupling $g_{ZA\gamma}$. This coupling is also constrained by the upper bound on $BR(Z \to \eta\gamma) < 5.1 \times 10^{-5}$ [61]. We have checked that in our cases this coupling is about four orders of magnitude below the bounds derived from [58, 61]. Searches for $H_{SM} \to A_1A_1 \to 4\gamma$ have been undertaken by ATLAS using $4.9 \text{ fb}^{-1}$ of integrated luminosity at 7 TeV c.m. energy in [62] for $M_{A_1} < 400$ MeV. One can assume that the corresponding upper bound on $BR(H_{SM} \to A_1A_1) \times BR(A_1 \to \gamma\gamma)^2 \lesssim 6.6 \times 10^{-3}$ applies to our scenario as well, which leads to $BR(H_{SM} \to A_1A_1) \lesssim 1.7 \times 10^{-2}$. If $A_1$ imitates a single photon, bounds on $BR(H_{SM} \to \gamma\gamma)$ should be respected. In our scenarios we require $BR(H_{SM} \to A_1A_1) \lesssim 5 \times 10^{-4}$, hence these constraints are well satisfied. Notably this small branching fraction has no impact on the measured signal rates of $H_{SM}$ into the other Standard Model channels, which agree well with the Standard Model predictions.

Additional constraints depending on $M_{A_1}$ will be discussed in the corresponding subsections below.

<footnote n="1">¹We thank F. Domingo for help for this calculation.</footnote>
## 3.1 $M_{A_1}$ near 210 MeV

For a light $A_1$, too light for hadronic final states ($M_{A_1}$ below $3M_\pi$), the possible decays are into $\mu^+\mu^-$, $e^+e^-$ and the loop induced decay into $\gamma\gamma$. The couplings of $A_1$ to Standard Model fermions are obtained via mixing with $A_{\text{MSSM}}$ as discussed in eqs. (2.14) and (2.15) in the previous section, and lead to a reduced coupling of $A_1$ to leptons $\sim X_d \sim \kappa$, see (2.19). These couplings determine also the partial width into $\gamma\gamma$. For a sizeable branching fraction into $\gamma\gamma$, the decay into $\mu^+\mu^-$ must be kinematically forbidden. On the other hand, for $M_{A_1} < 200$ MeV the remaining decays into $e^+e^-$ and $\gamma\gamma$ lead generically to a too small total width implying, for a boosted $A_1$ with an energy of about 375 GeV, a decay length larger than the size of the detectors (unless $A_1$ mixes strongly with $\pi^0$ as discussed in [38]). However, for $M_{A_1}$ very close to $2m_\mu$, the loop contribution of muons to the width $\Gamma(A_1 \to \gamma\gamma)$ reaches a maximum. It is given by (neglecting all other contributions; see, e.g., [63])

$$ \Gamma(A_1 \to \gamma\gamma)|_{\text{muons}} = \frac{G_\mu \alpha_{em}^2 M_{A_1}^3}{128 \sqrt{2} \pi^3} X_d^2 |A_{1/2}^A(\tau)|^2 \quad (3.1) $$

with $\tau = M_{A_1}^2/(4m_\mu^2)$ and, for $\tau \le 1$,

$$ A_{1/2}^A(\tau_f) = 2\tau^{-1} \arcsin^2 \sqrt{\tau}; \quad (3.2) $$

accordingly it increases with $M_{A_1} \to 2m_\mu$ (remaining finite for $M_{A_1} = 2m_\mu$). We find that, for $M_{A_1}$ near or slightly above 210 MeV, the partial width $\Gamma(A_1 \to \gamma\gamma)$ dominated by the muon contribution is large enough to dominate the width $\Gamma(A_1 \to e^+e^-)$ leading to a $BR(A_1 \to \gamma\gamma) \sim 74\%$.

The total $A_1$ width depends then essentially on its reduced coupling to muons $X_d$ related to $\kappa$ via (2.19). First we consider a scenario with a total width of $\sim 1.7 \times 10^{-13}$ GeV, leading to a decay length of $A_1$ for an energy of 375 GeV of about 2 m. Given that the distance of the EM calorimeter cells to the interaction point is larger than 1.3 m for the ATLAS and CMS detectors (depending on the angle $\eta$), one can estimate that somewhat more than 60% of all pseudoscalars decay before the EM calorimeter cells.

This scenario requires $\kappa \gtrsim 1.65$, in which case $\kappa$ runs into a Landau singularity at about 400 TeV where the NMSSM would require a UV completion (e.g. GMSB). Then a single Higgs state near 750 GeV is able to generate a visible signal. (The second Higgs state is heavier near 1 TeV and has a significantly smaller production cross section. A scenario where a single Higgs state near 750 GeV is responsible for the signal and another Higgs state is far below 750 GeV is not possible: then the lighter state would generate a larger signal, which is excluded.) For a large enough production cross section of the state near 750 GeV from its coupling to b-quarks it must have a dominant $H'$ (MSSM-like) component. Still, for a large enough branching fraction into $A_1A_1$, the $H' - S$ mixing angle (2.12) in the heavy scalar Higgs sector must not be too small and, notably, the coupling $g_{SA_1A_1}$ in (2.17) must be large. Both of these conditions are satisfied for $\kappa \sim 1.65$, which is required if a single state should generate a visible signal.

Suitable values of $\lambda$, $\tan\beta$ and $\mu$ for the desired masses and mixings are given by a benchmark point BP1 in table 1. (Since the mass of the second heavy Higgs state is near
<table><caption>**Table 1.** Parameters for the four benchmark points. The soft Susy breaking gaugino masses are $M_1 = 600$ GeV (500 GeV for BP1), $M_2 = 1$ TeV, $M_3 = 3$ TeV, all squarks are assumed degenerate, and all slepton masses are 300 GeV (with vanishing trilinear couplings). More digits (for all parameters) than shown here are necessary in order to reproduce $M_{A_1}$ given in the table 2 below, and in order to obtain $C_A$ sufficiently small for BP1 and BP2, see the text.</caption><thead><tr><th>Couplings and mass parameters</th><th>BP1</th><th>BP2</th><th>BP3</th><th>BP4</th></tr></thead><tbody><tr><td>λ</td><td>0.528</td><td>0.212</td><td>0.0332</td><td>0.0644</td></tr><tr><td>κ</td><td>1.65</td><td>0.75</td><td>0.121</td><td>0.215</td></tr><tr><td>tan β</td><td>9.57</td><td>16.8</td><td>15.5</td><td>14.5</td></tr><tr><td>μ (GeV)</td><td>138.5</td><td>101.1</td><td>102.3</td><td>111.3</td></tr><tr><td>A<sub>λ</sub> (GeV)</td><td>32.2</td><td>15.6</td><td>0.0</td><td>0.0</td></tr><tr><td>A<sub>κ</sub> (GeV)</td><td>1.16</td><td>7.67 × 10<sup>-2</sup></td><td>-4.69 × 10<sup>-4</sup></td><td>-1.49 × 10<sup>-3</sup></td></tr><tr><td>M<sub>squarks</sub> (TeV)</td><td>6</td><td>7.5</td><td>2</td><td>3</td></tr><tr><td>A<sub>t</sub> (TeV)</td><td>-3.48</td><td>-3.95</td><td>3</td><td>3</td></tr></tbody></table>

1 TeV and not near 750 GeV, these values deviate somewhat from the ones obtained in the previous section. Radiative corrections of $O(\kappa^2/4\pi^2 \sim 0.07)$ can require corresponding readjustments of these values.) Since $\tan\beta$ is $\sim 10$, the NMSSM-specific uplift of the Standard Model like Higgs mass at low $\tan\beta$ is not available. Then the Standard Model like Higgs mass of $\sim 125$ GeV requires large radiative corrections as in the MSSM.

As stated above and discussed in [38], the squark masses and $A_t$ can be chosen such that flavour violating couplings of $A_1$ are suppressed. In order to generate simultaneously large enough radiative corrections to the Standard Model like Higgs mass, both parameters have to be relatively large in the multi-TeV range. Possible numerical values are also indicated in table 1. The remaining NMSSM specific parameters $A_\lambda$ and $A_\kappa$ are chosen small, such that the $BR(H_{SM} \to A_1A_1)$ (depending somewhat on $A_\lambda$) is below $5 \times 10^{-4}$, and $M_{A_1}$ sufficiently close to $2m_\mu$ such that the total width of $A_1$ is large enough, i.e. that its decay length $l$ at 375 GeV is small enough: for the BP1 in table 1 with $M_{A_1} \sim 211.3$ MeV we get $\Gamma_{tot}(A_1) \sim 1.74 \times 10^{-13}$ GeV and $l \sim 2$ m, for which we estimate that $1 - e^{-d/l} \sim 63\%$ of all $A_1$ decays take place before the EM calorimeters. ($d$ denotes the average distance to the calorimeter cells of $\sim 2$ m.) For the production cross sections of the Higgs state $H_2$ at 750 GeV we find from SuShi 1.5.0 [64–71] (at NNLO with MMHT2014 PDFs) $\sigma_{ggF}(H_2) \sim 4.8$ fb, $\sigma_{bbH}(H_2) \sim 36.8$ fb, and from NMSSMTools we find $BR(H_2 \to A_1A_1) \sim 0.51$ with a total width of $H_2$ of $\sim 7$ GeV. Together with a $BR(A_1 \to \gamma\gamma) \sim 0.74$ we obtain a signal rate of $\sim 4.6$ fb. This signal rate remains to be multiplied by an acceptance $Acc(\gamma)$ for the diphotons to simulate a single photon in the detector. This issue will be discussed for all scenarios in section 3.3; for the time being the signal rates appear with a factor $Acc(\gamma)$ in table 2.

If we assume a slightly smaller value of $M_{A_1} \sim 210.5$ MeV, $\Gamma_{tot}(A_1)$ decreases to $\sim 1.65 \times 10^{-13}$ GeV leading to $l \sim 2.2$ m, reducing the percentage of decays before the EM calorimeters to $\sim 60\%$ and hence the signal rate by $\sim 10\%$.
<table><caption>**Table 2.** Higgs masses, production cross sections and branching fractions for the 4 benchmark points. For the points BP1 and BP2, the signal cross section takes into account losses from $A_1$ decays beyond 2 m according to a factor $(1 - e^{-2/l})^2$ with $l$ in m. For the points BP3 and BP4 the $BR(A_1 \to \gamma\gamma)$ includes the $BR(A_1 \to 3\pi^0)$. $Acc(\gamma)$ denotes the acceptance for the di- or multiphotons from two pseudoscalars to simulate two single photons in the detector as discussed in section 3.3. $M_{H_1}$ allows for a theoretical error of ~3 GeV. For BP4, $\Gamma_{tot}(A_1)$ and $l(A_1)$ are the unreliable parton level results (at NLO).</caption><thead><tr><th></th><th>BP1</th><th>BP2</th><th>BP3</th><th>BP4</th></tr></thead><tbody><tr><td><i>M</i><sub>H<sub>1</sub></sub> (GeV)</td><td>122.1</td><td>124.3</td><td>123.7</td><td>122.2</td></tr><tr><td><i>M</i><sub>H<sub>2</sub></sub> (GeV)</td><td>750</td><td>730</td><td>744</td><td>740</td></tr><tr><td><i>M</i><sub>H<sub>3</sub></sub> (GeV)</td><td>1003</td><td>762</td><td>750</td><td>750</td></tr><tr><td><i>M</i><sub>A<sub>1</sub></sub> (MeV)</td><td>211.3</td><td>211.1</td><td>548.7</td><td>510.3</td></tr><tr><td><i>M</i><sub>A<sub>2</sub></sub> (GeV)</td><td>763</td><td>747</td><td>748</td><td>745</td></tr><tr><td><i>M</i><sub>H±</sub> (GeV)</td><td>757</td><td>749</td><td>752</td><td>749</td></tr><tr><td>σ<sub>ggF</sub>(<i>H</i><sub>2</sub>) (fb)</td><td>4.8</td><td>2.2</td><td>1.7</td><td>1.9</td></tr><tr><td>σ<sub>bbH</sub>(<i>H</i><sub>2</sub>) (fb)</td><td>36.8</td><td>67.2</td><td>44.7</td><td>44.9</td></tr><tr><td>σ<sub>ggF</sub>(<i>H</i><sub>3</sub>) (fb)</td><td>0.1</td><td>1.8</td><td>2.0</td><td>1.9</td></tr><tr><td>σ<sub>bbH</sub>(<i>H</i><sub>3</sub>) (fb)</td><td>0.2</td><td>52.5</td><td>54.3</td><td>44.3</td></tr><tr><td>BR(<i>H</i><sub>2</sub> → <i>A</i><sub>1</sub><i>A</i><sub>1</sub>)</td><td>0.51</td><td>0.66</td><td>0.082</td><td>0.21</td></tr><tr><td>BR(<i>H</i><sub>3</sub> → <i>A</i><sub>1</sub><i>A</i><sub>1</sub>)</td><td>0.72</td><td>0.53</td><td>0.048</td><td>0.16</td></tr><tr><td>BR(<i>A</i><sub>1</sub> → γγ)</td><td>0.74</td><td>0.73</td><td>0.72</td><td>0.66</td></tr><tr><td>Γ<sub>tot</sub>(<i>A</i><sub>1</sub>) (10<sup>-13</sup> GeV)</td><td>1.74</td><td>0.65</td><td>7500</td><td>19000</td></tr><tr><td><i>l</i>(<i>A</i><sub>1</sub>) at 375 GeV</td><td>2 m</td><td>5.5 m</td><td>0.18 mm</td><td>0.08 mm</td></tr><tr><td>Signal cross section (fb)</td><td>4.6 ×Acc(γ)</td><td>3.7 ×Acc(γ)</td><td>3.4 ×Acc(γ)</td><td>6.7 ×Acc(γ)</td></tr></tbody></table>

Scenarios with smaller values of $\kappa$ are also possible. Then, however, the reduced coupling $X_d$ of $A_1$ to leptons is smaller (see (2.19)), and the total $A_1$ width decreases. Hence the decay length increases, and a smaller fraction of $A_1$'s decay before 2 m. This loss can be compensated for if two states $H_2$ and $H_3$ with large production cross sections and branching fractions into $A_1A_1$ contribute to the signal.

The benchmark point BP2 is of this type, where we take $\kappa = 0.75$, nearly (but not quite) small enough for the absence of a Landau singularity below the GUT scale. For $M_{A_1} \sim 211.1$ MeV the total $A_1$ width is $\sim 6.5 \times 10^{-14}$ GeV, leading to $l \sim 5.5$ m. We estimate that then only $\sim 30\%$ of all $A_1$ decays take place before the EM calorimeter cells. On the other hand, two Higgs states $H_2$ and $H_3$ with masses near 730 GeV and 762 GeV contribute to the signal. Both are strong mixtures of the pure MSSM-like and singlet-like states. For $H_2$, the production cross section is $\sigma_{ggF+bbH}(H_2) \sim 69.4$ fb, and $BR(H_2 \to A_1A_1) \sim 0.66$. For $H_3$, the production cross section is $\sigma_{ggF+bbH}(H_3) \sim 54.3$ fb, and $BR(H_2 \to A_1A_1) \sim 0.53$. Together with a $BR(A_1 \to \gamma\gamma) \sim 0.73$ we obtain a signal rate of $\sim 3.7$ fb times $Acc(\gamma)$, as shown in table 2.
## 3.2 $M_{A_1}$ at 510–550 MeV

The partial widths of a light pseudoscalar in this mass range can be estimated employing two complementary approaches. To begin with one can ask what one would obtain within the parton model, extrapolated into the nonperturbative domain of QCD. First, for a reduced coupling of $A_1$ to leptons $X_d \sim 0.1$ as considered below, the partial width of $A_1$ into muons can still be computed reliably and is

$$ \Gamma(A_1 \to \mu^+\mu^-) \sim 5 \times 10^{-11} \text{ GeV}. \qquad (3.3) $$

The loop induced partial width of $A_1$ into $\gamma\gamma$ is $\sim 4 \times 10^{-15}$ GeV and hence negligibly small. At NLO QCD the partial width of $A_1$ into strange quarks is about $5 \times 10^{-10}$ GeV and the loop induced width into gluons of the same order as the width into $\mu^+\mu^-$. These widths can only be rough estimates, however.

An alternative approach is to consider the case $M_{A_1} \approx 550$ MeV, where one can expect that $A_1$ mixes with the $\eta$ meson with a mass of 547.85 MeV. (The possible rôle of $\eta$ for the decays of a light pseudoscalar has been indicated earlier in [54] without quantitative statements, however.) Mixing with the $\pi^0$ meson of a lighter $A_1$ with $M_{A_1} \approx 135$ MeV has been considered in [38], where Partial Conservation of Axial Currents (PCAC) or the sigma model for light mesons is employed in order to determine the off-diagonal element of the $A_1$-meson mass matrix; the same formalism will be used here for $A_1 - \eta$ mixing for $M_{A_1} \approx 550$ MeV.

First we discuss this latter case, where the results can be considered as more reliable. Only subsequently we turn to the case $M_{A_1} \approx 510$ MeV, motivated by the alleviation of constraints from radiative $\Upsilon(1S)$ decays in this mass range, see below. There, however, estimates of partial widths of $A_1$ are more speculative.

For $M_{A_1} \approx 550$ MeV, the relevant mass matrix of the $A_1 - \eta$ system reads in the basis $(A_1, \eta)$

$$ \frac{1}{2} \begin{pmatrix} M_{A_1}^2 & \delta m_{A_1 \eta}^2 \\ \delta m_{A_1 \eta}^2 & m_{\eta}^2 \end{pmatrix}. \qquad (3.4) $$

For a small mixing angle $\theta$,

$$ \theta \sim \frac{\delta m_{A_1 \eta}^2}{M_{A_1}^2 - m_{\eta}^2} \ll 1, \qquad (3.5) $$

the eigenstate $A'_1$ contains a small $\eta$ component: $A'_1 \sim A_1 + \theta\eta + \dots$. For the partial widths of $A'_1$ one obtains then

$$ \Gamma(A'_1 \to X) \approx \Gamma(A_1 \to X) + \theta^2 \Gamma(\eta \to X). \qquad (3.6) $$

The dominant $\eta$ decays are [61]

$$ \begin{gathered} BR(\eta \to \gamma\gamma) \sim 39\%, \quad BR(\eta \to 3\pi^0) \sim 33\%, \quad BR(\eta \to \pi^+\pi^-\pi^0) \sim 23\%, \\ \Gamma_{\text{tot}}(\eta) \sim 1.3 \times 10^{-6} \text{ GeV}. \end{gathered} \qquad (3.7) $$

Next we require that the $\eta$-induced decays into $\gamma\gamma$ or $3\pi^0$ of the eigenstate $A'_1$ dominate its width into $\mu^+\mu^-$, since we ignore the unreliable widths of $A_1$ into strange quarks or
gluons in this subsection. (Since the latter decays can also generate $\gamma\gamma$ or $3\pi^0$ final states, this assumption is conservative.) This leads to

$$ \theta^2 > \frac{\Gamma(A_1 \to \mu^+\mu^-)}{\Gamma_{\text{tot}}(\eta)}, \quad \theta \gtrsim 6 \times 10^{-3}. \qquad (3.8) $$

In order to estimate the mixing matrix element $\delta m^2_{A_1\eta}$ above we use, following [38], PCAC. There one introduces the SU(3)$_{axial}$ flavour currents $J^\mu_{Ai}$ where $i$ denote the SU(3) generators. Assuming that $\eta$ is a pure octet, $J^\mu_{A8}$ satisfies

$$ \partial_\mu J^\mu_{A8} = f_\pi m_\eta^2 \eta \qquad (3.9) $$

with $f_\pi \sim 93$ MeV. At the quark level one has

$$ \partial_\mu (\bar{s}\gamma^\mu \gamma_5 s) = -\sqrt{\frac{2}{3}} \partial_\mu J^\mu_{A8} + \frac{1}{\sqrt{3}} \partial_\mu J^\mu_{A0} \qquad (3.10) $$

where $J^\mu_{A0}$ is the (anomalous) U(1)$_A$ current whose divergence involves the $\eta'$ meson. Using these relations, one can re-write the coupling of $A_1$ to strange quarks in the Lagrangian (proportional to the corresponding Yukawa coupling $X_d m_s / v$)

$$ \frac{-im_s X_d}{\sqrt{2}v} A_1 \bar{s} \gamma_5 s = -\frac{X_d}{2\sqrt{2}v} A_1 \partial_\mu (\bar{s} \gamma^\mu \gamma_5 s) = \frac{X_d}{2\sqrt{3}v} A_1 \partial_\mu J^\mu_{A8} + \dots = \frac{X_d f_\pi m_\eta^2}{2\sqrt{3}v} \eta A_1 + \dots \qquad (3.11) $$

where we have dropped the terms $\sim \partial_\mu J^\mu_{A0} \sim \eta'$. From (3.11) one can read off

$$ \delta m^2_{A_1\eta} = \frac{X_d f_\pi m_\eta^2}{2\sqrt{3}v}. \qquad (3.12) $$

Then the request (3.8) becomes, again for $X_d \sim 0.1$ and using (3.5),

$$ |M_{A_1} - m_\eta| < 10^{-3} m_\eta \sim 0.5 \text{ MeV}. \qquad (3.13) $$

This estimate can be refined by including mixing with the $\eta'$ meson, the anomalous U(1)$_A$ current $J^\mu_{A0}$ and the loop-induced coupling of $A_1$ to $F_{\mu\nu}\tilde{F}^{\mu\nu}$, where $F_{\mu\nu}$ is the QCD field strength.<fnref n="2" /> The additional contribution to $\delta m^2_{A_1\eta}$ leads to a replacement of the right hand side of (3.13) by $\sim 1$ MeV.

Assuming such a small $A_1 - \eta$ mass difference, the decay length of $A_1$ is below a mm, and its branching fractions are the ones of $\eta$ given in (3.7) above. For $M_{A_1} > M_K$, constraints from rare K decays are no longer relevant. However, since $A_1$ has couplings to b-quarks $\sim X_d$, constraints from the search for the radiative decays $\Upsilon(1S) \to \gamma\eta$ by CLEO in [49] apply (and are more relevant than searches for $\Upsilon(1S) \to \gamma A_1 \to \gamma\mu^+\mu^-$).

The $BR(\Upsilon(1S) \to \gamma A_1)$ can be obtained from the Wilczek formula [72, 73]

$$ \frac{BR(\Upsilon(1S) \to \gamma A_1)}{BR(\Upsilon(1S) \to \mu^+\mu^-)} = \frac{G_F m_b^2 X_d^2}{\sqrt{2}\pi\alpha_{em}} \left(1 - \frac{M_{A_1}^2}{M_{\Upsilon(1S)}^2}\right) \times F, \quad \text{hence} \\ BR(\Upsilon(1S) \to \gamma A_1) \sim 1.03 \times 10^{-4} \times X_d^2 \qquad (3.14) $$

<footnote n="2">²We thank F. Domingo for providing us with his studies of this issue.</footnote>
where $BR(Υ(1S) → μ^+μ^-) ∼ 2.48\%$ and $F$ is a correction factor $∼ 0.5$. The upper bound of CLEO [49] on $BR(Υ(1S) → γη)$ is $1.0 × 10^{-6}$ at the 90% CL level, or $1.3 × 10^{-6}$ at the 95% CL level. Applying this bound to the $BR(Υ(1S) → γA_1)$, (3.14) gives

$$X_d \lesssim 0.11 \qquad (3.15)$$

as used above. From (2.19) one finds that $κ$ must then also be quite small, leading to relatively small branching fractions of the heavy Higgs states $H_2$ and $H_3$ into $A_1A_1$. Hence both of these states should contribute to the signal.

Next we consider the $A_1$ decays induced by its mixing with $η$ where $η$ decays as in (3.7). The decays into $π^0π^+π^-$ give diphotons plus muons, but due to the escaping neutrinos this final state will not allow to reconstruct the masses of the original resonances near 750 GeV. In addition to the $A_1 → γγ$ mode, the $A_1 → 3π^0$ mode leads to photon jets. The compatibility of such photon jets with a single photon signature in the detectors has been discussed in detail in [35]. In particular, due to the enhanced probability for photon conversions into $e^+ + e^-$ in the inner parts of the detectors, such scenarios can be distinguished from single photons (or even diphotons) once more events are available. Adding both modes, about 72% of all $A_1$ decays lead to di- or multi-photons. The resulting signal cross section remains to be multiplied by the acceptance $Acc(γ)$ for the di- or multi-photons to fake a single photon discussed in section 3.3.

The parameters, masses, branching fractions, production and signal cross sections of a corresponding benchmark point BP3 are shown in tables 1 and 2.

If $M_{A_1}$ differs by a few tens of MeV from the $η$ mass it becomes more difficult to estimate its decays; its mixing angle with the on-shell $η$ meson using PCAC as above becomes tiny. Its Yukawa couplings to Standard Model fermions are obtained through mixing with the (heavy) MSSM-like pseudoscalar $A_{\text{MSSM}}$. At the parton level and for $tan β ∼ 10–15$, the relative couplings squared and hence the corresponding widths of $A_1$ are dominantly into $s\bar{s}$ ($≈ 87\%$), into $gg$ via top quark loops ($≈ 5\%$), and into $μ^+μ^-$ ($≈ 8\%$). The hadronic or $γγ$ decays of $A_1$ can then be considered as being mediated by the CP-odd isospin and color singlet interpolating composite fields $s\bar{s}$ and $F_{μν}\tilde{F}^{μν}$. Both are known components of the $η$ wave function in Fock space, and the most reasonable assumption is that their hadronisation (decays into physical hadrons and $γγ$) proceeds again with branching fractions similar to the ones of $η$.

The partial width for the sum of these decays of $A_1$ is less clear, however. It is relevant, since it competes with the width of $A_1$ into $μ^+μ^-$ and determines consequently the branching fraction $BR(A_1 → \text{hadrons or } γγ)$ via the above interpolating fields relative to the $BR(A_1 → μ^+μ^-)$. Since the widths for the above mentioned decays of $η$ into $γγ$ or pions are small (being electromagnetic or suppressed by isospin), one must assume that the widths for the decays of $A_1$ via the above interpolating fields are also smaller than estimated from the couplings squared at the parton level as at the beginning of this section. A quantitative statement is difficult, however, without a nonperturbative evaluation of the relevant matrix elements between physical states.

On the other hand, the sum of the couplings squared of $A_1$ to $s\bar{s}$ or $gg$ and hence the sum of the partial widths of $A_1$ into $s\bar{s}$ or $gg$ (for both of which $η$-like branching fractions are
assumed) is considerably larger than into $\mu^+\mu^-$: at NLO one has $BR(A_1 \to s\bar{s} \text{ or } gg) \sim 0.92$. Hence, reducing the sum of the partial widths of $A_1$ into $s\bar{s}$ or $gg$ by a factor 1/10 leaves us still with a dominant $BR(A_1 \to s\bar{s} \text{ or } gg) \sim 0.55$.

In the scenario where $M_{A_1}$ differs by a few tens of MeV from the $\eta$ mass we will make the assumption that the reduction of the width of the decays $A_1 \to \gamma\gamma$ or hadrons is not too dramatic, i.e. the relevant branching fractions of $A_1$ can be parametrized as

$$BR(A_1 \to \gamma\gamma, 3\pi^0, \pi^0\pi^+\pi^-) \sim F_A \times BR(\eta \to \gamma\gamma, 3\pi^0, \pi^0\pi^+\pi^-) \quad (3.16)$$

where the factor $F_A$ is not too small ($\gtrsim 0.1$).

Let us have another look at the searches by CLEO which were performed separately for the $\eta \to 3\pi^0$, $\eta \to \pi^0\pi^+\pi^-$ and $\eta \to \gamma\gamma$ final states. The windows for the invariant masses were chosen differently for different final states, $M_{3\pi^0} \gtrsim 475$ MeV and $M_{\pi^0\pi^+\pi^-} \gtrsim 515$ MeV. ($M_{\gamma\gamma}$ is fitted to a double Gaussian function centered at $M_\eta$.) No candidates were found in the (background free) $3\pi^0$ and $\pi^0\pi^+\pi^-$ final states, but two events in $\pi^0\pi^+\pi^-$ with $M_{\pi^0\pi^+\pi^-} \sim 510$ MeV just below the $M_{\pi^0\pi^+\pi^-}$ window. Also a mild excess of events for $M_{\gamma\gamma} \sim 510$ MeV is observed. These events are not numerous enough to allow for the claim of a signal, but we conclude that the $\pi^0\pi^+\pi^-$ and $\gamma\gamma$ final states do not lead to stronger upper limits on $BR(\Upsilon(1S) \to \gamma A_1)$ for $M_{A_1} \sim 510$ MeV than the limit on $BR(\Upsilon(1S) \to \gamma\eta)$ from the remaining $3\pi^0$ final state. After translating the 90% CL upper limit from the latter final state into a 95% CL upper limit, we find from only the $3\pi^0$ final state

$$BR(\Upsilon(1S) \to \gamma\eta) \times 0.33 \lesssim 3.8 \times 10^{-6} \quad (3.17)$$

where 0.33 is the $BR(\eta \to 3\pi^0)$. For $A_1$, assuming $\eta$-like decays, this upper bound becomes

$$BR(\Upsilon(1S) \to \gamma A_1) \times 0.33 \times F_A \lesssim 3.8 \times 10^{-6}. \quad (3.18)$$

Combining (3.18) and (3.14) gives

$$X_d \lesssim 0.19/\sqrt{F_A} \quad (3.19)$$

and, from (2.19),

$$\kappa \lesssim 0.21/\sqrt{F_A}. \quad (3.20)$$

Since smaller $F_A$ alleviates the constraint on $\kappa$, we should study its impact on the signal rate. First, the branching fractions of the heavy Higgs states induced by the coupling $g_{SA_1A_1}$ in (2.17) behave roughly like $g_{SA_1A_1}^2 \sim \kappa^4 \lesssim 0.002 F_A^{-2}$, i.e. smaller $F_A$ allows for larger branching fractions. On the other hand, by assumption (see (3.16)) the branching fractions of $A_1$ into $\gamma\gamma$ or $3\pi^0$ are proportional to $F_A$. Hence the factors of $F_A$ cancel approximatively in the final signal rate. (Depending on the other parameters we found, however, that smaller $F_A$ can lead to a decrease of the signal rate if the heavy Higgs branching fractions increase somewhat less than indicated above.)

For a rough estimate we have constructed a benchmark point BP4 with $M_{A_1} \sim 510$ MeV. It has $X_d \simeq 0.206$ and satisfies the CLEO constraints for $F_A < 0.87$. Its parameters are given in table 1. For the branching fractions of $A_1$ into $\gamma\gamma$ or $3\pi^0$ we assume
$F_A \sim 0.87$ which gives, using the corresponding branching fractions of $\eta$, $BR(A_1 \to \gamma\gamma$ or $3\pi^0) \sim 0.87 \times (0.39 + 0.33) \simeq 0.63$. Together with the production cross sections and branching fractions of the Higgs states $H_2$ and $H_3$ in table 2 we obtain finally a signal cross section of $\sim 6.7$ fb times $\text{Acc}(\gamma)$.

### 3.3 The spread in rapidity of multiphotons

The probability for a di- or multiphoton system to fake a single photon depends to a large extent on its angular spread. In the context of a decay of the Standard Model Higgs boson into light pseudoscalars (which decay into diphotons) this has been discussed in some detail in [76]. One has to consider the fineness in rapidity $\eta$ of the strips of the first layer of the electromagnetic calorimeter (EM), ranging from 0.003 to 0.006 (depending on $\eta$) for ATLAS [77]. A relevant quantity for an event is the fraction of the total deposited energy in a single strip, and which fraction is deposited in the adjacent strips [76]. In the case of diphotons, the relevant criterium is then the distribution of $\Delta\eta$ between the two photons. In [76] it is argued (and used in [34]) that only for $\Delta\eta \lesssim 0.0015$ prompt diphotons fake a single photon. (Including converted photons does not seem to modify this estimate [76].) This number remains to be confirmed by the experimental collaborations, however, and will depend on $\eta$ and the transverse energy of the individual events in practice.

In the case of displaced vertices, as for our benchmark points BP1 and BP2, the situation is more involved as discussed in [34]: diphotons from pseudoscalars decaying between the original vertex and the EM are more collimated. Note that the signal cross sections given in table 2 (without the factor $\text{Acc}(\gamma)$) take already into account the loss from pseudoscalars decaying inside or beyond the EM. The interplay between these losses and $\Delta\eta$ for diphotons reaching the EM has been studied in [34], but is beyond the scope of the present paper. We content ourselves with the fact that $\text{Acc}(\gamma)$ for our benchmark points BP1 and BP2 is then definitively larger than $\text{Acc}(\gamma)$ for pseudoscalars of a corresponding mass of $\sim 200$ MeV which decay promptly. To this end we studied the distribution of $\Delta\eta$ for diphotons from promptly decaying pseudoscalars with the help of a simulation based on MadGraph/MadEvent [78] including Pythia 6.4 [79], where the pseudoscalars originate from a 750 GeV Higgs state. The resulting distribution is shown in red in figure 1.

We find that 70% of these diphotons satisfy $\Delta\eta \lesssim 0.0015$ (in rough agreement with [34]), and correspondingly more due to the displaced vertices for our benchmark points BP1 and BP2. A rough estimate for $\text{Acc}(\gamma)$ for the diphotons from both pseudoscalars to fake a single photon would then be $\approx 75\%$. We underline, however, that the actual number of signal events responsible for the observed excesses is small, and statistical fluctuations of quantities like $\Delta\eta$ can be correspondingly large.

Turning to the scenarios BP3 and BP4, displaced vertices are no longer relevant. However, here we expect about equal fractions of both diphoton and 6-photon final states, the latter from $A_1 \to 3\pi^0$. For $M_{A_1} \sim 500$ MeV, the spread in $\Delta\eta$ for the diphotons is obviously larger, as shown in figure reffig:1 in green. Here only $\sim 25\%$ of all diphotons satisfy $\Delta\eta \lesssim 0.0015$.

For the 6-photon final state we note that the angular spread between the photons should be smaller, as the total invariant mass of the system must remain the same and
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Distributions of $\Delta\eta$ of diphotons from promptly decaying pseudoscalars with a mass of 200 MeV (red) and 500 MeV (green). 500 MeV, 2/6 $\gamma$ (blue) denotes $\Delta\eta$ between the two leading among all 6 photons from the decays $A_1 \to 3\pi^0 \to 6\gamma$.</figcaption></figure>

less energy is available for momenta transverse to the principal axis. Since the relevant quantity in the EM calorimeters is the spread of the deposited energy, we concentrate here on $\Delta\eta$ between the two most energetic photons resulting from a single $A_1 \to 3\pi^0$ decay. It is shown in blue in figure 1, and we find that in $\sim 75\%$ of all cases these are closer than $\Delta\eta \sim 0.0015$. This number is suggestive, but it is not clear whether it coincides with the fraction of 6-photon final states faking a single photon; due to the more homogenous distribution of the deposited energy this fraction could even be larger. If we use it as it is, the $\text{Acc}(\gamma)$ for the di- or multiphotons from both pseudoscalars to fake simultaneously a single photon becomes $\text{Acc}(\gamma) \approx 25\%$. Moreover many of 6 photons will convert, and the acceptance of such events remains to be studied by the experimental collaborations. Again, statistical fluctuations can be large as long as the number of signal events is as low as at present.

# 4 Summary and conclusions

We have shown that the excess of events in the diphoton final state near 750 GeV observed by ATLAS and CMS can be explained within a fairly standard supersymmetric extension of the Standard Model, the NMSSM, without invoking new particles like additional vector-like quarks and/or leptons. The signal cross sections are not very large, but may be sufficient to explain the observed excesses.

The corresponding processes differ, however, from what has been proposed in most of the literature up to now: except for the scenario BP1 (with $\kappa \sim 1.65$), two resonances nearby in mass which share the properties of the additional CP-even scalars of the NMSSM
are responsible for the signal cross section. Their components proportional to the MSSM-like scalar $H'$ lead to enhanced couplings to b-quarks implying sizeable production cross sections via associated production with b-quarks, whereas their components proportional to the singlet-like scalar $S$ lead to sizeable branching fractions into two pseudoscalars $A_1A_1$. These scenarios are not in tension with the upper limit from CMS on the diphoton cross section obtained at the run I at 8 TeV [4]; note that the $b\bar{b}H$ cross sections increase somewhat faster with the c.m. energy than $ggF$.

We note that constraints from other decay modes of the scalars with masses of about 750 GeV are satisfied: upper bounds on signal cross sections into other final states — quark pairs, lepton pairs and electroweak gauge bosons — are discussed in [5, 74, 75]. These bounds are obeyed given the relatively large branching fractions into $A_1A_1$ and finally into di- or multi-photon final states of the heavy scalars in our scenario, which do not require excessive production cross sections.

Four different scenarios have been discussed, which differ in the properties and masses of the light pseudoscalars $A_1$ and the heavy Higgs states:

1) For the benchmark points BP1 and BP2, the mass $M_{A_1} \sim 211$ MeV is just below twice the muon mass. Then the branching fraction of $A_1$ into diphotons is large enough for a satisfactory signal rate. In the case of the BP1 with $\kappa \sim 1.65$, a single heavy Higgs state (still a mixture of the MSSM-like and singlet-like states) with a width of $\sim 7$ GeV is sufficient for a signal. In the case of the BP2 with a more modest value for $\kappa \sim 0.75$, two nearby heavy Higgs states, both with a width of $\sim 6$ GeV, are responsible for the signal. For both BP1 and BP2, the dominant constraints from low energy experiments originate from $K$ decays involving loop-induced flavour changing vertices of $A_1$; it must be assumed that these are cancelled by suitable choices of the SUSY breaking parameters.

2) For the benchmark points BP3 and BP4 it is assumed that $A_1$ shares its branching fractions with the $\eta$ meson. In the case of BP3 with $M_{A_1} \sim 549$ MeV this is guaranteed by $A_1 - \eta$ mixing, estimated with the help of the PCAC formalism. In the case of BP4 with $M_{A_1} \sim 510$ MeV, estimates of the $A_1$ partial widths are on less solid ground. We assumed that the non-leptonic decays of $A_1$ proceed via $s\bar{s}$ or $gg \sim F_{\mu\nu}\tilde{F}^{\mu\nu}$ interpolating fields which, in turn, hadronise (decay) again similar to the $\eta$ meson. We showed, however, that sizeable reductions of the corresponding partial widths with respect to the decays of $A_1$ into $s\bar{s}$ or $gg$ by, e.g., a factor 1/10, would not invalidate this scenario. For both BP3 and BP4 two nearby heavy Higgs states with widths of $\sim 1.5-2$ GeV are responsible for the signal. For this range of $M_{A_1}$ the dominant constraints from low energy experiments originate from searches for radiative $\Upsilon$ decays into $\gamma + \eta$ by CLEO. These lead to upper bounds on the coupling of $A_1$ to down-type quarks and leptons and, as we have shown, on $\kappa$.
Interestingly, the four scenarios have different features which allow to distinguish them experimentally also from more “conventional” models:

1) For BP2, BP3 and BP4 the signal originates from two resonances $H_2$ and $H_3$ close in mass, which can imitate a single wide resonance. Of course, small variations of the parameters allow to vary the masses of $H_2$ and $H_3$, the total signal rate, and to reshuffle the individual signal rates of $H_2$ and $H_3$. With more events (and depending on the actual mass difference) the two states could possibly be resolved. A particular feature of BP1 is that the single resonance near 750 GeV responsible for the signal has another large branching fraction of ~25% into $Z + A_1$. With $A_1$ imitating a photon, one obtains signals of the kind $Z+\gamma$ similar to the ones expected if a 750 GeV resonance decays into $\gamma\gamma$ via fermionic loops.

2) The $A_1$ decays differ considerably for the benchmark points. For BP1 and BP2 the decay lengths of $A_1$ are macroscopic leading to measurable displaced vertices if $A_1$ decays inside the calorimeters. BP2 corresponds to an extreme case with a decay length of ~5.5 m, but a very large signal cross section (before reducing it by the number of decays before the EM calorimeters). For both BP1 and BP2, $A_1$ decays into diphotons. However, due to the displaced vertices it will not be straightforward to distinguish them from single photons via the number of converted photons [35]. Moreover, $A_1$ has branching fractions of ~25-30% into $e^+ + e^-$ leading to similar signatures as converted photons. Also the opening angle between the two photons gets reduced for displaced vertices increasing their acceptance as a single photon, see section 3.3. For BP3 and BP4 the $A_1$ decay lengths are short, but $A_1$ decays into diphotons or photon jets from $3\pi^0$. The latter should lead to a very large proportion to “converted photons”; additional potentially relevant observables like EM shower shapes have also been discussed in [35]. Our results for $M_{A_1} \sim 500$ MeV for $\Delta\eta$ between diphotons or the two leading photons from jets from $3\pi^0$ may be useful here. Finally, for BP3 and BP4 $A_1$ has branching fractions into muon pairs (of 6-8% at the parton level) which could be used for alternative signals, once more events are obtained.

Hence, if the excess of events continues, several observables can be used to verify/test/exclude the scenarios discussed here.

Finally we recall that the origin of the light pseudoscalar $A_1$ in the NMSSM is an approximate R-symmetry of the scalar potential, see the small values of $A_\kappa$ and $A_\lambda$ of the benchmark points. For BP1, BP2 and BP3, $M_{A_1}$ has to coincide accidentally with specific values $2m_\mu$ or $m_\eta$. For the BP4 $M_{A_1}$ is actually less constrained (unless one intends to fit the events near $M_{A_1} \sim 510$ MeV observed by CLEO as we did), but in this scenario the $A_1$ decays are less understood theoretically. Although the approximate R-symmetry at the weak (or SUSY) scale is not preserved by radiative corrections, we content ourselves in the present paper with the mere fact that such a scenario would allow to explain the events. Work on an R-symmetric extension of the NMSSM explaining a light pseudo-Goldstone boson naturally is in progress.
Acknowledgments

The authors acknowledge the support of France Grilles for providing cloud computing
resources on the French National Grid Infrastructure and thank J. Jaeckel and M. Spira for
helpful discussions, and notably F. Domingo for communication on possible $A_1-\eta$ mixings.
U.E. acknowledges support from the European Union Initial Training Networks HiggsTools
(PITN-GA-2012-316704), INVISIBLES (PITN-GA-2011-289442), the ERC advanced grant
Higgs@LHC, and from the grant H2020-MSCA-RISE-2014 No. 645722 (NonMinimalHiggs).

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

References

[1] ATLAS collaboration, *Search for resonances decaying to photon pairs in 3.2 fb⁻¹ of pp collisions at √s = 13 TeV with the ATLAS detector*, ATLAS-CONF-2015-081.

[2] CMS collaboration, *Search for new physics in high mass diphoton events in proton-proton collisions at √s = 13 TeV*, CMS-PAS-EXO-15-004.

[3] ATLAS collaboration, *Search for high-mass diphoton resonances in pp collisions at √s = 8 TeV with the ATLAS detector*, Phys. Rev. D 92 (2015) 032004 [arXiv:1504.05511] [INSPIRE].

[4] CMS collaboration, *Search for diphoton resonances in the mass range from 150 to 850 GeV in pp collisions at √s = 8 TeV*, Phys. Lett. B 750 (2015) 494 [arXiv:1506.02301] [INSPIRE].

[5] S. Knapen, T. Melia, M. Papucci and K. Zurek, *Rays of light from the LHC*, Phys. Rev. D 93 (2016) 075020 [arXiv:1512.04928] [INSPIRE].

[6] R. Franceschini et al., *What is the γγ resonance at 750 GeV?*, JHEP 03 (2016) 144 [arXiv:1512.04933] [INSPIRE].

[7] J. Ellis, S.A.R. Ellis, J. Quevillon, V. Sanz and T. You, *On the Interpretation of a Possible ~750 GeV Particle Decaying into γγ*, JHEP 03 (2016) 176 [arXiv:1512.05327] [INSPIRE].

[8] B. Bellazzini, R. Franceschini, F. Sala and J. Serra, *Goldstones in Diphotons*, JHEP 04 (2016) 072 [arXiv:1512.05330] [INSPIRE].

[9] A. Falkowski, O. Slone and T. Volansky, *Phenomenology of a 750 GeV Singlet*, JHEP 02 (2016) 152 [arXiv:1512.05777] [INSPIRE].

[10] A. Angelescu, A. Djouadi and G. Moreau, *Scenarii for interpretations of the LHC diphoton excess: two Higgs doublets and vector-like quarks and leptons*, Phys. Lett. B 756 (2016) 126 [arXiv:1512.04921] [INSPIRE].

[11] D. Bečirevič, E. Bertuzzo, O. Sumensari and R. Zukanovich Funchal, *Can the new resonance at LHC be a CP-Odd Higgs boson?*, Phys. Lett. B 757 (2016) 261 [arXiv:1512.05623] [INSPIRE].

[12] X.-F. Han and L. Wang, *Implication of the 750 GeV diphoton resonance on two-Higgs-doublet model and its extensions with Higgs field*, Phys. Rev. D 93 (2016) 055027 [arXiv:1512.06587] [INSPIRE].
[13] W.-C. Huang, Y.-L.S. Tsai and T.-C. Yuan, *Gauged Two Higgs Doublet Model confronts the LHC 750 GeV di-photon anomaly*, arXiv:1512.07268 [INSPIRE].

[14] S. Moretti and K. Yagyu, *750 GeV diphoton excess and its explanation in two-Higgs-doublet models with a real inert scalar multiplet*, Phys. Rev. D 93 (2016) 055043 [arXiv:1512.07462] [INSPIRE].

[15] M. Badziak, *Interpreting the 750 GeV diphoton excess in minimal extensions of Two-Higgs-Doublet models*, arXiv:1512.07497 [INSPIRE].

[16] L.J. Hall, K. Harigaya and Y. Nomura, *750 GeV Diphotons: Implications for Supersymmetric Unification*, JHEP 03 (2016) 017 [arXiv:1512.07904] [INSPIRE].

[17] N. Bizot, S. Davidson, M. Frigerio and J.L. Kneur, *Two Higgs doublets to explain the excesses* pp → γγ(750 GeV) and h → τ$^{±}$µ$^{∓}$, JHEP 03 (2016) 073 [arXiv:1512.08508] [INSPIRE].

[18] E. Ma, *Diphoton Revelation of the Utilitarian Supersymmetric Standard Model*, arXiv:1512.09159 [INSPIRE].

[19] X.-F. Han, L. Wang, L. Wu, J.M. Yang and M. Zhang, *Explaining 750 GeV diphoton excess from top/bottom partner cascade decay in two-Higgs-doublet model extension*, Phys. Lett. B 756 (2016) 309 [arXiv:1601.00534] [INSPIRE].

[20] A.E.C. Hernández, I.d.M. Varzielas and E. Schumacher, *The 750 GeV diphoton resonance in the light of a 2HDM with $S_3$ flavour symmetry*, arXiv:1601.00661 [INSPIRE].

[21] X.-F. Han, L. Wang and J.M. Yang, *An extension of two-Higgs-doublet model and the excesses of 750 GeV diphoton, muon $g-2$ and $h \to \mu\tau$, Phys. Lett. B 757 (2016) 537* [arXiv:1601.04954] [INSPIRE].

[22] S.F. King and R. Nevzorov, *750 GeV Diphoton Resonance from Singlets in an Exceptional Supersymmetric Standard Model*, JHEP 03 (2016) 139 [arXiv:1601.07242] [INSPIRE].

[23] E. Bertuzzo, P.A.N. Machado and M. Taoso, *Di-Photon excess in the 2HDM: hasting towards the instability and the non-perturbative regime*, arXiv:1601.07508 [INSPIRE].

[24] R. Ding, L. Huang, T. Li and B. Zhu, *Interpreting 750 GeV Diphoton Excess with R-parity Violation Supersymmetry*, arXiv:1512.06560 [INSPIRE].

[25] B.C. Allanach, P.S.B. Dev, S.A. Renner and K. Sakurai, *Di-photon Excess Explained by a Resonant Sneutrino in R-parity Violating Supersymmetry*, arXiv:1512.07645 [INSPIRE].

[26] Y.-L. Tang and S.-h. Zhu, *NMSSM extended with vector-like particles and the diphoton excess on the LHC*, arXiv:1512.08323 [INSPIRE].

[27] F. Wang, W. Wang, L. Wu, J.M. Yang and M. Zhang, *Interpreting 750 GeV diphoton resonance as degenerate Higgs bosons in NMSSM with vector-like particles*, arXiv:1512.08434 [INSPIRE].

[28] W. Chao, *The Diphoton Excess Inspired Electroweak Baryogenesis*, arXiv:1601.04678 [INSPIRE].

[29] M. Badziak, M. Olechowski, S. Pokorski and K. Sakurai, *Interpreting 750 GeV Diphoton Excess in Plain NMSSM*, arXiv:1603.02203 [INSPIRE].

[30] P. Agrawal, J. Fan, B. Heidenreich, M. Reece and M. Strassler, *Experimental Considerations Motivated by the Diphoton Excess at the LHC*, arXiv:1512.05775 [INSPIRE].
[31] J. Chang, K. Cheung and C.-T. Lu, Interpreting the 750 GeV diphoton resonance using photon jets in hidden-valley-like models, Phys. Rev. D 93 (2016) 075013 [arXiv:1512.06671] [INSPIRE].

[32] M. Chala, M. Duerr, F. Kahlhoefer and K. Schmidt-Hoberg, Tricking Landau-Yang: How to obtain the diphoton excess from a vector resonance, Phys. Lett. B 755 (2016) 145 [arXiv:1512.06833] [INSPIRE].

[33] X.-J. Bi et al., A Promising Interpretation of Diphoton Resonance at 750 GeV, arXiv:1512.08497 [INSPIRE].

[34] L. Aparicio, A. Azatov, E. Hardy and A. Romanino, Diphotons from Diazions, arXiv:1602.00949 [INSPIRE].

[35] B. Dasgupta, J. Kopp and P. Schwaller, Photons, Photon Jets and Dark Photons at 750 GeV and Beyond, arXiv:1602.04692 [INSPIRE].

[36] C.-W. Chiang, H. Fukuda, M. Ibe and T.T. Yanagida, 750 GeV diphoton resonance in a visible heavy QCD axion model, arXiv:1602.07909 [INSPIRE].

[37] G. Arcadi, P. Ghosh, Y. Mambrini and M. Pierre, Re-opening dark matter windows compatible with a diphoton excess, arXiv:1603.05601 [INSPIRE].

[38] F. Domingo, S. Heinemeyer, J.S. Kim and K. Rolbiecki, The NMSSM lives: with the 750 GeV diphoton excess, Eur. Phys. J. C 76 (2016) 249 [arXiv:1602.07691] [INSPIRE].

[39] M. Maniatis, The Next-to-Minimal Supersymmetric extension of the Standard Model reviewed, Int. J. Mod. Phys. A 25 (2010) 3505 [arXiv:0906.0777] [INSPIRE].

[40] U. Ellwanger, C. Hugonie and A.M. Teixeira, The Next-to-Minimal Supersymmetric Standard Model, Phys. Rept. 496 (2010) 1 [arXiv:0910.1785] [INSPIRE].

[41] ATLAS, CMS collaborations, Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at $\sqrt{s}$ = 7 and 8 TeV, ATLAS-CONF-2015-044 [CMS-PAS-HIG-15-002].

[42] Q.-H. Cao, Y.-Q. Gong, X. Wang, B. Yan and L.L. Yang, One bump or two peaks: The 750 GeV diphoton excess and dark matter with a complex mediator, Phys. Rev. D 93 (2016) 075034 [arXiv:1601.06374] [INSPIRE].

[43] D.J. Miller, R. Nevzorov and P.M. Zerwas, The Higgs sector of the next-to-minimal supersymmetric standard model, Nucl. Phys. B 681 (2004) 3 [hep-ph/0304049] [INSPIRE].

[44] L.J. Hall and T. Watari, Electroweak supersymmetry with an approximate $U(1)_{PQ}$, Phys. Rev. D 70 (2004) 115001 [hep-ph/0405109] [INSPIRE].

[45] R. Barbieri, L.J. Hall, A.Y. Papaioannou, D. Pappadopulo and V.S. Rychkov, An alternative NMSSM phenomenology with manifest perturbative unification, JHEP 03 (2008) 005 [arXiv:0712.2903] [INSPIRE].

[46] B.A. Dobrescu and K.T. Matchev, Light axion within the next-to-minimal supersymmetric standard model, JHEP 09 (2000) 031 [hep-ph/0008192] [INSPIRE].

[47] R. Dermisek and J.F. Gunion, The NMSSM Close to the R-symmetry Limit and Naturalness in $h \to aa$ Decays for $m_a < 2m_b$, Phys. Rev. D 75 (2007) 075019 [hep-ph/0611142] [INSPIRE].

[48] D.E. Morrissey and A. Pierce, Modified Higgs Boson Phenomenology from Gauge or Gaugino Mediation in the NMSSM, Phys. Rev. D 78 (2008) 075029 [arXiv:0807.2259] [INSPIRE].
[49] CLEO collaboration, S.B. Athar et al., Search for radiative decays of Υ(1S) into η and η', Phys. Rev. D 76 (2007) 072003 [arXiv:0704.3063] [INSPIRE].

[50] U. Ellwanger, J.F. Gunion and C. Hugonie, NMHDECAY: A Fortran code for the Higgs masses, couplings and decay widths in the NMSSM, JHEP 02 (2005) 066 [hep-ph/0406215] [INSPIRE].

[51] U. Ellwanger and C. Hugonie, NMHDECAY 2.0: An updated program for sparticle masses, Higgs masses, couplings and decay widths in the NMSSM, Comput. Phys. Commun. 175 (2006) 290 [hep-ph/0508022] [INSPIRE].

[52] G. Hiller, B physics signals of the lightest CP odd Higgs in the NMSSM at large tan beta, Phys. Rev. D 70 (2004) 034018 [hep-ph/0404220] [INSPIRE].

[53] S. Andreas, O. Lebedev, S. Ramos-Sanchez and A. Ringwald, Constraints on a very light CP-odd Higgs of the NMSSM and other axion-like particles, JHEP 08 (2010) 003 [arXiv:1005.3978] [INSPIRE].

[54] B.A. Dobrescu, G.L. Landsberg and K.T. Matchev, Higgs boson decays to CP odd scalars at the Tevatron and beyond, Phys. Rev. D 63 (2001) 075003 [hep-ph/0005308] [INSPIRE].

[55] F. Domingo and U. Ellwanger, Updated Constraints from B Physics on the MSSM and the NMSSM, JHEP 12 (2007) 090 [arXiv:0710.3714] [INSPIRE].

[56] F. Domingo and U. Ellwanger, Constraints from the Muon g-2 on the Parameter Space of the NMSSM, JHEP 07 (2008) 079 [arXiv:0806.0733] [INSPIRE].

[57] F. Domingo, U. Ellwanger, E. Fullana, C. Hugonie and M.-A. Sanchis-Lozano, Radiative Upsilon decays and a light pseudoscalar Higgs in the NMSSM, JHEP 01 (2009) 061 [arXiv:0810.4736] [INSPIRE].

[58] J. Jaeckel and M. Spannowsky, Probing MeV to 90 GeV axion-like particles with LEP and LHC, Phys. Lett. B 753 (2016) 482 [arXiv:1509.00476] [INSPIRE].

[59] F. Domingo, Update of the flavour-physics constraints in the NMSSM, arXiv:1512.02091 [INSPIRE].

[60] CLEO collaboration, W. Love et al., Search for Very Light CP-Odd Higgs Boson in Radiative Decays of Upsilon(S-1), Phys. Rev. Lett. 101 (2008) 151802 [arXiv:0807.1427] [INSPIRE].

[61] PARTICLE DATA GROUP collaboration, K.A. Olive et al., Review of Particle Physics, Chin. Phys. C 38 (2014) 090001 [INSPIRE].

[62] ATLAS collaboration, Search for a Higgs boson decaying to four photons through light CP-odd scalar coupling using 4.9 fb⁻¹ of 7 TeV pp collision data taken with ATLAS detector at the LHC, ATLAS-CONF-2012-079.

[63] J.R. Ellis, M.K. Gaillard and D.V. Nanopoulos, A Phenomenological Profile of the Higgs Boson, Nucl. Phys. B 106 (1976) 292 [INSPIRE].

[64] R.V. Harlander, S. Liebler and H. Mantler, SusHi: A program for the calculation of Higgs production in gluon fusion and bottom-quark annihilation in the Standard Model and the MSSM, Comput. Phys. Commun. 184 (2013) 1605 [arXiv:1212.3249] [INSPIRE].

[65] R.V. Harlander and W.B. Kilgore, Next-to-next-to-leading order Higgs production at hadron colliders, Phys. Rev. Lett. 88 (2002) 201801 [hep-ph/0201206] [INSPIRE].
[66] R.V. Harlander and W.B. Kilgore, Higgs boson production in bottom quark fusion at next-to-next-to leading order, Phys. Rev. D 68 (2003) 013001 [hep-ph/0304035] [INSPIRE].

[67] G. Degrassi and P. Slavich, NLO QCD bottom corrections to Higgs boson production in the MSSM, JHEP 11 (2010) 044 [arXiv:1007.3465] [INSPIRE].

[68] G. Degrassi, S. Di Vita and P. Slavich, NLO QCD corrections to pseudoscalar Higgs production in the MSSM, JHEP 08 (2011) 128 [arXiv:1107.0914] [INSPIRE].

[69] G. Degrassi, S. Di Vita and P. Slavich, On the NLO QCD Corrections to the Production of the Heaviest Neutral Higgs Scalar in the MSSM, Eur. Phys. J. C 72 (2012) 2032 [arXiv:1204.1016] [INSPIRE].

[70] R. Harlander and P. Kant, Higgs production and decay: Analytic results at next-to-leading order QCD, JHEP 12 (2005) 015 [hep-ph/0509189] [INSPIRE].

[71] S. Liebler, Neutral Higgs production at proton colliders in the CP-conserving NMSSM, Eur. Phys. J. C 75 (2015) 210 [arXiv:1502.07972] [INSPIRE].

[72] F. Wilczek, Problem of Strong p and t Invariance in the Presence of Instantons, Phys. Rev. Lett. 40 (1978) 279 [INSPIRE].

[73] R.N. Hodgkinson, $Υ \to \gamma A_1$ in the NMSSM at large $\tan \beta$, Phys. Lett. B 665 (2008) 219 [arXiv:0802.3197] [INSPIRE].

[74] M. Low, A. Tesi and L.-T. Wang, A pseudoscalar decaying to photon pairs in the early LHC Run 2 data, JHEP 03 (2016) 108 [arXiv:1512.05328] [INSPIRE].

[75] R.S. Gupta, S. Jäger, Y. Kats, G. Perez and E. Stamou, Interpreting a 750 GeV Diphoton Resonance, arXiv:1512.05332 [INSPIRE].

[76] P. Draper and D. McKeen, Diphotons from Tetraphotons in the Decay of a 125 GeV Higgs at the LHC, Phys. Rev. D 85 (2012) 115023 [arXiv:1204.1061] [INSPIRE].

[77] ATLAS collaboration, Measurement of the inclusive isolated prompt photon cross section in pp collisions at $\sqrt{s} = 7$ TeV with the ATLAS detector, Phys. Rev. D 83 (2011) 052005 [arXiv:1012.4389] [INSPIRE].

[78] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, MadGraph 5: Going Beyond, JHEP 06 (2011) 128 [arXiv:1106.0522] [INSPIRE].

[79] T. Sjöstrand, S. Mrenna and P.Z. Skands, PYTHIA 6.4 Physics and Manual, JHEP 05 (2006) 026 [hep-ph/0603175] [INSPIRE].