Probing the two light Higgs scenario in the NMSSM with a low-mass pseudoscalar

In this article we propose a simultaneous collider search strategy for a pair of scalar bosons in the NMSSM through the decays of a very light pseudoscalar. The massive scalar has a mass around 126 GeV while the lighter one can have a mass in the vicinity of 98 GeV (thus explaining an apparent LEP excess) or be much lighter. The successive decay of this scalar pair into two light pseudoscalars, followed by leptonic pseudoscalar decays, produces clean multi-lepton final states with small or no missing energy. Furthermore, this analysis offers an alternate leptonic probe for the 126 GeV scalar that can be comparable with the ZZ* search channel. We emphasize that a dedicated experimental search for multi-lepton final states can be an useful probe for this scenario and, in general, for the NMSSM Higgs sector. We illustrate our analysis with two representative benchmark points and show how the LHC configuration with 8 TeV center-of-mass energy and 25 inverse femtobarns of integrated luminosity can start testing this scenario, providing a good determination of the light pseudoscalar mass and a relatively good estimation of the lightest scalar mass.


Introduction
The experimental evidence for a scalar boson with a mass around 126 GeV [1,2] has been received with enthusiasm by the particle physics community, since it could correspond to the missing piece of the standard model (SM), the Higgs boson. Although most of the properties of this new particle (in particular the branching ratios of various decay modes) are compatible with a SM-like Higgs, an apparent di-photon excess [3,4], has been considered as a possible hint for physics beyond the SM, such as supersymmetry (SUSY). However, this excess remains yet to be confirmed [5].
It has been pointed out that in order to reproduce the mass of the Higgs boson in the Minimal Supersymmetric Standard Model (MSSM) one generally needs a very heavy spectrum, raising concerns about the naturalness of this construction. This problem can be alleviated in extended Supersymmetric constructions. In particular, the Next-to-Minimal Supersymmetric Standard Model (NMSSM) can accommodate a 126 GeV scalar [6][7][8][9][10][11], while maintaining the idea of naturalness [11][12][13][14]. The advantage of the NMSSM resides in the presence of an extra singlet field, which provides an additional contribution to the tree-level Higgs mass [15][16][17]. Thus, contrary to what happens in the MSSM, loop corrections can be smaller and the correct Higgs mass is achieved for a lighter supersymmetric spectrum. This extra singlet superfieldŜ is introduced to promote the bilinear µ-term µĤ dĤu in the MSSM to a trilinear coupling λŜĤ dĤu . After electroweak symmetry breaking (EWSB) takes place, the vacuum expectation value (VEV) of the singlet field triggers an effective µ parameter, λv s , which is naturally of the order of the electroweak scale, thus curing the µ-problem of the MSSM [18]. On top of this, the singlet-doublet mixing of the resulting scalar states is useful to produce the observed di-photon excess [6,19].
The presence of light scalar and pseudoscalar Higges can lead to characteristic Higgs-to-Higgs cascade decays [20,21]. Although these cascades are generally suppressed, a substantial enhancement can occur when particles are light and therefore induce distinctive signatures. 1 For example, if the lightest pseudoscalar, denoted as a 0 1 , has a mass m a 0 1 > ∼ 2m b , the process h 0 1 → 2a 0 1 → 2bb with 4b final states is dominant [22,[36][37][38][39][40]. This seems a good channel for collider study, 2 however, in hadron colliders like the LHC this signature is hidden by a huge QCD background. A viable alternative is looking for 4τ [20,22,25,31,35,39,41,42] or 2µ2τ [42] final states, for which the problem of large QCD backgrounds is somehow ameliorated. However, this approach suffers from the poor tau-identification efficiency, dependent on the transverse momentum [43]. Besides, since tau decays are always accompanied by neutrinos, there is no sharp peak of a 0 1 invariant mass. For these reasons, as discussed in the literature, 4µ [39,44,45], 4γ [21,33] or 2γ + 2-jets [46] search channels with m a 0 1 < 2m τ turn out to be the ideal way to probe NMSSM Higgses with very light pseudoscalars.
Recently, the possibility of having two light scalar Higgses in NMSSM has received a lot of attention. One of these Higgses, h 0 2 , would correspond to the scalar observed by ATLAS and CMS with a mass around 126 GeV and the other one, h 0 1 , with a mass around 98 GeV [10], would be consistent with the small excess in the LEP search for e + e − → Zh, h → bb [47]. 3 Alternative scenarios with an even lighter h 0 1 have also been considered earlier [22,50]. Since the recent LHC observation of the 126 GeV scalar, which is well compatible with a SM-like Higgs boson, many of these studies have become extremely constrained, although a small window for new physics is still open if the di-photon excess is confirmed. On the other hand, scenarios with very light pseudoscalars are also affected by stringent experimental bounds. This is the case for m a 0 1 < 2m b (which leads to leptonic final states 4µ, 2µ2τ , and 4τ ) [26,34,42,51], with tight constraints in the case m a 0 1 < 2m τ for CPeven Higgs masses in the range of 86 -150 GeV by the recent CMS analysis [52], especially in the case of singlet-like states. The latter is a consequence of the fact that BR(h 0 1 /h 0 2 → a 0 1 a 0 1 ) ∼ 1 and thus σ(pp → h 0 1 /h 0 2 → 2a 0 1 → 4µ) is sizable and would have already been observed. 4 In this work we investigate potential detection channels for scenarios involving very light pseudoscalar particles and two light scalar Higgses in the NMSSM. For concreteness, we consider a SM-like scalar Higgs in the range 124 GeV < m h 0 2 < 128 GeV (consistent with LHC findings) and a lighter h 0 1 in the mass range 96 -100 GeV (to account for the LEP excess) or lighter than 86 GeV (to avoid the current CMS limit). In order to avoid the above mentioned constraints that affect the lightest pseudoscalar, a 0 1 is assumed to have a mass in the range 2m τ m a 0 1 < 2m b , since the µ + µ − final state is then a sub-leading leptonic mode with a branching ratio of the order of 10 −2 and σ(pp → h 0 1 → 2a 0 1 → 4µ) is small enough to elude the CMS limit of the 4µ search. We point out that the SM-like Higgs decay modes h 0 2 → 2a 0 1 → 4µ, 2µ2τ , and 4τ have branching fractions comparable to those for the same final states from h 0 2 → ZZ * and therefore provide alternative four-lepton 3 A two light Higgs but without a low-mass pseudoscalar has been discussed in the MSSM [48], as well as in the NMSSM [49]. 4 One way to avoid this is to consider the case when the production cross section for a singlet-like h 0 1 /h 0 2 is very suppressed. signatures accompanied with zero or a small missing energy at the LHC. We also note that the long Higgs-to-Higgs cascades of h 0 which constitutes a clean and distinctive multi-lepton signal at the LHC. We design a set of cuts that isolates the signal events from the background, study the distributions of the kinematic variable M T2 [53], and determine to what extent the masses of the pseudoscalar and scalar Higgs particles can be reconstructed. In our analysis we include constraints from B-physics [54][55][56][57], the muon anomalous magnetic moment [58][59][60][61], and the LHC result for h 0 2 → γγ. In addition, we impose the WMAP bound on the lightest neutralino relic abundance [62] and we ensure that its spin-dependent scattering cross section is consistent with XENON100 data [63].
The paper is organized as follows. We start with a brief description of the model and determine the relevant signatures in Sec. 2. The choice of model parameters is justified in Sec. 3 where we also describe the details of the analysis. Sec. 4 is dedicated to the analysis of the numerical results. We also introduce some discriminating variables to suppress SM backgrounds and we investigate the statistical significance of the proposed signals. Finally we summarize our conclusions in Sec. 5.

Higgs signals with a low-mass pseudoscalar in the NMSSM
The superpotential of the Z 3 -symmetric NMSSM is given by where W ′ is the MSSM superpotential without the µ-term,Ŝ is a superfield singlet under the SM gauge group, and a, b = 1, 2 are SU(2) L indices. Similarly, with L ′ soft representing the MSSM soft-terms without the B µ -term, the Lagrangian containing the soft SUSY-breaking terms in a supergravity framework is given by The last trilinear term in Eq. (1) is essential to avoid an unacceptable axion associated with the breaking of a global U(1) symmetry [16]. After EWSB, this term also generates an effective Majorana mass 2κv s for the singlino-like neutralino. The second term of Eq. (1) not only generates an effective µ parameter, µ eff = λv s , but also provides an extra tree-level contribution to the lightest doublet-like Higgs boson mass. The complete expression is now [15][16][17], where tan β = v u /v d is the ratio of up and down-type Higgs VEVs. Thus even maintaining perturbativity (λ 0.7), the extra contribution to the tree-level Higgs mass can be sizable for small values of tan β. The last term of Eq. (2) is important to understand the phenomenology of singlet-like scalars. In particular, when the lightest CP-even scalar h 0 1 as well as the lightest CP-odd scalar a 0 1 are predominantly singlet-like, the decay width for h 0 1 → a 0 1 a 0 1 is proportional to A κ κ. As already mentioned in the previous section, the Higgs sector of the NMSSM is richer than that of the MSSM, featuring three CP-even scalar states, (h 0 1 , h 0 2 , h 0 3 ) and two pseudoscalars, (a 0 1 , a 0 2 ). This is of particular interest for studying Higgs cascade decays, which can now be long and often result in multi-particle final states with no missing energy. The decay chains become more interesting in the presence of light scalar and pseudoscalar Higgses, which can be in agreement with present constraints if their singlet composition is dominant. In particular, a viable hierarchy for the scalar mass eigenstates in the NMSSM consists of having a light singlet-like h 0 1 , a SM-like Higgs h 0 2 with a mass around 126 GeV and a much heavier h 0 3 which is also doublet-like. This structure with two light scalar Higgses might produce distinctive features in collider searches, as compared to the signals that come from the MSSM Higgs sector, and therefore we concentrate on this possibility.
Although very light singlet-like pseudoscalars can give rise to potentially characteristic signatures, one must be careful with recent experimental constraints. In particular, for a scalar Higgs mass in the mass range of 86 -150 GeV, the upper has been found to be 3 -5 fb, depending on the pseudoscalar mass for 1 GeV < m a 0 1 < 2m τ in the 7 TeV CMS search [52]. Thus, if we want to study the 98 + 126 GeV scenario for the scalar Higgs bosons, in order to avoid this constraint, we will consider a light pseudoscalar with a mass in the range 2m τ m a 0 1 < 2m b . We call this Scenario I. Another way to avoid the CMS bound is to consider a lighter scalar Higgs with m h 0 1 < 86 GeV. Although this would relax the LHC limit on the pseudoscalar mass, in fact 2m τ < m a 0 1 < 2m b is still very constrained from the LEP limit on 4τ [64] and multi-lepton/jet final states [36,65] in the case when m h 0 1 m h 0 2 /2. We have included these constraints in our analysis and checked that viable points can still be obtained within this mass range for the lightest pseudoscalar. We call this Scenario II.
For clearness we summarise here the properties of the two kind of scenarios considered in our work We have also checked that these benchmark scenarios are in agreement with ATLAS and CMS direct pseudoscalar searches, pp → a 0 1 → µ + µ − [66], due to the singlet nature of a 0 1 and its dominant decay branching fraction to τ . We are now ready to specify the potential signatures that we will study for these scenarios. We consider the following decay modes, where ℓ = e, µ, τ . We here assume that h 0 1 and h 0 2 have been produced through gluon-fusion and consider inclusive decay modes of the tau lepton. The missing energy / E T is associated with neutrinos coming from tau-decays and is generally small. It is important to emphasize that direct-pseudoscalar decays to a muon-pair will eventually yield a clean four-muon final state with vanishing missing energy and the only source of electrons is leptonic τ -decay, since BR(a 0 1 → e + e − ) ∼ 0. The second decay mode can act as an alternative source of four-lepton final state apart from h 0 2 → ZZ * → 4ℓ. In the chosen corner of parameter space, we found that the predicted BR(h 0 2 → a 0 1 a 0 1 ) is comparable with that of BR(h 0 2 → ZZ * ), which still matches well the SM prediction. We thus encourage our experimental colleagues to search for possible excesses in 4µ, 4τ or 2µ2τ channels by relaxing the Z-boson resonance condition. In the remainder of the paper we carry out a dedicated analysis including the discussion of SM backgrounds and useful collider variables to search for the four-lepton collider signals of these decay modes.
Concerning possible SUSY backgrounds, these generally lead to large missing energy and can therefore be distinguished with a cut on / E T . There can also be sources of four-lepton final states with a small or no missing energy, provided that one considers R-parity violation [67] in the NMSSM [68] or µνSSM [69]. It is well known that R-parity conserving SUSY models are characterized by large missing energy due to the production of a heavy neutralino LSP, which is stable and escapes the detector, whereas this is not the case in R-parity violating scenarios. Models with broken R-parity and non-minimal superfield content, on the other hand, exhibit moderate or large displaced vertices with non-prompt jets or leptons in the final states (see, for example, Refs. [70,71]) with more complex Higgs cascade decays [71] and occasional correlations with neutrino physics [72] following [73], which can be used for discrimination. We also note that four leptons can be produced in h 0 3 decays, but with a very reduced production cross section.
Notice that one distinguishing feature of Scenario II is a three-step cascade decay of the Higgs boson h 0 2 → 2h 0 1 → 4a 0 1 → multi-lepton/jet. We have checked for several points in the parameter space that, even when m h 0 can be as large as ∼ 10 −3 and thus multi-lepton final states can be experimentally accessible. Due to the large number of leptons and/or jets in the final state, the main backgrounds are expected to be SUSY processes (such as cascade decays mediated by neutralinos and/or charginos), rather than SM ones. The analysis of this decay signal will be carried out in a separate work [74] since the search strategy is very different than for the other two.

Experimental constraints and the choice of parameters
The recent discovery of a 126 GeV scalar along with the reported di-photon excess set stringent limitations on the NMSSM parameter space. In addition, considering bounds from low-energy observables in flavour physics, the supersymmetric contribution to the muon anomalous magnetic moment a SUSY µ , and the relic abundance of the dark matter, which in our case would correspond to the lightest neutralino, narrow down significantly the viable regions.
In our analysis we impose compatibility with all these constraints. The numerical results are obtained with nmssmtools 3.2.1 [75], which calculates the mass spectrum and provides predictions for low-energy observables, as well as computing the decay widths of Higgses [76] and sparticles [77]. The masses of Higgs bosons include full two-loop contributions. We also consider the lightest neutralino as a dark matter candidate and include the WMAP upper constraint on its relic abundance [62] as well as the upper bound on the spin-independent neutralino-nucleon scattering cross section, σ SI , from XENON100 [63]. The neutralino relic density and its scattering cross section are computed through an interface with micromegas [78]. We have also modified the code nmssmtools to include three-body Higgs decay in the case of the NMSSM.
We have calculated the theoretical predictions for BR(B 0 s → µ + µ − ) and compared them with the experimental value. In the case of the MSSM, the SUSY contributions to this observable [79] are likely to exceed the recent LHCb measure- [56], especially for light pseudoscalars in the large tan β regime. In the case of the NMSSM [80,81], the corresponding Wilson coefficient receives contributions from both pseudoscalar Higgs bosons through their doublet components. Given that the light pseudoscalar is a purely singlet-like field, it leads to a negligible effect on this observable, and since a 0 2 is heavy enough and tan β is small, the experimental constraint is easily fulfilled. On the other hand, in general, BR(b → sγ) provides a strong constraint on SUSY models. On top of the usual MSSM terms, NMSSM-specific contributions arise from the extended Higgs and the neutralino sectors, which come into effect at the two-loop level [80]. This observable has been shown to lead to stringent constraints on the NMSSM parameter space [30] and this is indeed the case in our current analysis. Here, we consider the experimental result [57] and include the theoretical error of the calculation in the SM [82] −1600, 1000, 1800 −1600, 1000, 1300 The SUSY contribution to the muon anomalous magnetic moment, a SUSY µ was also computed. The observed discrepancy between the experimental value [58] and the SM predictions using e + e − data favours positive contributions from new physics in the range 10.1 × 10 −10 < a SUSY µ < 42.1 × 10 −10 at the 2σ confidence level [59][60][61], combining experimental and theoretical errors in quadrature. If tau data is employed, this discrepancy is smaller, 2.9 × 10 −10 < a SUSY µ < 36.1 × 10 −10 [60].
We have carried out a simple scan to sample the NMSSM parameter space, searching for regions in the parameter space where Scenarios I and II could be realised. Out of the viable regions we have selected two benchmark points, BP1 and BP2 with parameters given in Table 1, which constitute representative examples of Scenario I and Scenario II, respectively. The squarks and gluino masses for the first two generation are chosen to be heavy enough to be in agreement with current LHC SUSY searches [83]. We have relaxed gaugino mass unification in order to have more freedom in the neutralino composition, which is important in order to fix its relic abundance and scattering cross section. Also, trilinear parameters are non-universal and chosen in such a way that a SUSY µ is maximized and BR(b → sγ) is in agreement with the experimental value. The resulting mass spectra are shown in Table 2.
The compositions of Higgses, neutralinos, and charginos for our benchmark points are shown in Table 3. Since the µ parameter is chosen to be small and λ ∼ 0.3, the singlet VEV is around 430 GeV. Thus with the small κ value, the singlino mass, 2κv s for the benchmark points BP1 and BP2 is approximately 72 and 100 GeV, respectively. Notice that this value is close to that of the µ parameter and that both of them are smaller than the bino and wino mass terms. This results in an interesting hierarchical neutralino spectrum, especially for the three lightest neutralino states. The lightest neutralino is a singlino-Higgsino state (where the singlino component slightly dominates). There is an orthogonal eigenstate which is also singlino-Higgsino   third neutralino and in BP2 is the second, and close in mass to this one we can find an almost pure Higgsino ( χ 0 2 for BP1 and χ 0 3 for BP2). The heavier states are bino and wino-like.
The resulting predictions for low-energy observables are shown in Table 4. Although the values of a SUSY µ are in the 2σ range of the result from tau data, there is some tension with the one from e + e − data. We have tried reducing this discrepancy by decreasing the masses in the slepton sector. Regarding the neutralino relic abundance, it is well known that a Higgsino-like neutralino generally entails a large annihilation cross section and consequently a small relic abundance, but in our case, the presence of a sizable (slightly dominant as we said above) singlino component is welcome to make it compatible with the upper constraint on the dark matter relic density. Similarly, the elastic scattering cross section for neutralinos with a large Higgsino component is large and can be in conflict with current experimental  σ(gg→h SM →γγ) > 0.8 [10,11], within the 2σ range of the ATLAS and the CMS result [1][2][3][4]. For BR(b → sγ), the theoretical error is added to the experimental one in quadrature.
bounds. In our case, the recent constraints from XENON100 are very severe and exclude a large portion of the parameter space we analyzed. The theoretical predictions for the spin-independent neutralino-nucleon scattering cross section in both BP1 and BP2 are below the current exclusion line which, for a mass in the range 60 -100 GeV, is of the order of < ∼ 10 −9 pb. Notice in this sense that increasing the singlino component through the decrease, e.g., of the κ parameter may be useful to avoid this problem but this alternative might be questionable since it leads to a reduction in the di-photon production from the Higgs decay.

Collider analysis
Let us finally present the collider studies for the 4ℓ + / E T signal of the benchmarks shown in Tables 1 and 2. As discussed in Sec. 2, the benchmark points can be classified in terms of the Higgs masses, namely, In these benchmark points, BR(h 0 2 → a 0 1 a 0 1 ) is comparable to or larger than BR(h 0 2 → ZZ * ), and BR(h 0 1 → a 0 1 a 0 1 ) ∼ 1. The light pseudoscalar boson a 0 1 decays mainly into the di-tau final state with a branching fraction of ∼ 79 (92) % in Scenario I (II). Given that the pseudoscalar mass is larger than 2m τ , the di-muon decay mode is negligible in Scenario II, while it remains sub-leading for Scenario I. Thus, restating the expression (4), the leading final state will be with ℓ = e, µ, or τ h . Here, electrons and muons are coming from leptonic tau decays, while τ h denotes the tau jet originated from hadronic tau decays. Occasional muons can come from a 0 1 → µ + µ − decay process, which is sub-leading or negligible in the chosen benchmarks. The source of the missing energy is associated with neutrinos resulting from tau decays. In order to determine the feasibility of the signal at the LHC run with the 8 TeV center-of-mass energy ( √ s), which has recently finished its operation, we consider inclusive search channels. In other words, our analysis differs from past studies which usually consider one or two exclusive channels such as 4µ or 2µ + 2τ h with the 14 TeV beam condition.
We have generated Monte Carlo (MC) event samples of the Higgs signals for a proton-proton collision at √ s = 8 TeV using Herwig++ 2.6.1 [86] with the CTEQ6L1 parton distribution function (PDF) [87]. h 0 2 → ZZ * → 4ℓ processes are also generated since the final states are similar as in Eq. (5). The generated event samples have been scaled to 25 fb −1 of integrated luminosity, which corresponds to the LHC data accumulated in the year 2012. To simplify the analysis, we here consider leading-order cross sections for all the processes. The production cross section and the decay width of the SM Higgs boson are calculated with higlu [88], then the NMSSM Higgs cross sections have been obtained by corrections according to the total decay widths and BR(h 0 1, 2 → gg) calculated with nmssmtools. The generator-level events are further processed with the fast detector simulation program Delphes 2.0.3 [89] using a modified CMS detector card. In the detector simulation, jets are reconstructed by using the anti-k t algorithm [90] with the radius parameter of 0.5, and the b-tagging efficiency is set to be 70%. We assumed that the mis-tagging rate for c-jet is 10%, and that for the gluon and light-flavor jets is 1%, both of which are taken into account by Delphes. A candidate jet must satisfy p T > 25 GeV for its transverse momentum and |η| < 2.5 for its pseudorapidity. Isolated electrons and muons are required to have p T > 8 and 6 GeV, respectively, and |η| < 2.4. A tau jet is accepted only when it has p T > 15 GeV. In order to increase the purity of the leptonic signal, any charged lepton, including the tau jet, lying within a distance ∆R ℓj < 0.4 from a candidate jet is discarded in the analysis, where ∆R ab is defined as ( The dominant SM backgrounds include Drell-Yan (DY), bb, cc, semi-or dileptonically decaying tt, W W/Z, and ZZ/γ processes. Direct J/ψ or Υ productions can in principle contribute, but they are found to be almost negligible. To estimate the backgrounds, bb and cc event samples are generated with Pythia 6.4 [91] using CTEQ6L1 PDF in the various p T bins, and the remaining dominant background processes are generated by Herwig++ at the matrix-element level. The parton showering and the hadronization are performed by Herwig++ for all the back- ground processes. Then, the MC samples have been fed into Delphes using the same detector card used for the Higgs signals. In order to suppress the backgrounds while keeping the significance of the Higgs signals of interest, the following basic event selection cuts are imposed.
• At least two pairs of oppositely-charged leptons including the tau jet, • no b-tagged jet.
When selecting the leptons, the priority is given to the hard electrons and isolated muons since the reconstruction efficiency of the tau jet is relatively poor. No cut is applied on the missing energy, despite the presence of neutrinos in the final state from tau decays. In fact, the neutrino momenta are approximately collinear with their parent tau momentum, leading to a partial cancellation among them and resulting in a quite small missing energy for signal events. To illustrate this, the signal distribution of the missing energy is compared to some dominant background distributions in the left panel of Fig. 1. As already mentioned in the previous section, multi-lepton and missing energy signatures of SUSY cascade decay processes (not included here) can also be serious backgrounds to the Higgs signal. However, the missing energy in these processes is expected to be much larger and one might attempt to use an upper cut on / E T . In any case, as we will see below, all these backgrounds can be efficiently suppressed by employing other cut variables.
For the events that passed the basic selection cuts, the combinatorial ambiguity of pairing the leptons is resolved by calculating all the possible M T2 defined as in Ref. [53], where M T of the neutrino system, which are determined by a numerical minimization over all possible splittings that maintain the missing energy condition constructed in the event. The M T2 variable was introduced for a decay topology like a pair production of the squark that decays into a quark and the LSP. The decay topology given in expression (5) is essentially different from that because there are at least four neutrinos in each event. However, the M T2 is applicable to the type of Higgs signal events considered here since the pseudoscalar bosons are emitted back-to-back in the rest frame of the scalar Higgs boson. Thus their decay products are nearly collinear and the neutrinos sharing the same parent pseudoscalar boson can be considered as one invisible particle. In the construction of the M T2 , there are two possible ways of pairing the leptons because two pairs of oppositely-charged leptons are initially selected. We choose the pairing which gives the smaller value of M T2 . This method can be justified by the fact that the M T2 value would be bounded from above by the parent particle mass m a 0 1 if the right pairing is chosen, whereas there is no such a restriction if the pairing was wrong [93]. Since we here consider a light pseudoscalar boson, an upper cut on the M T2 < 25 GeV is further implemented. The M T2 distribution is shown in Fig. 2 for both benchmark points. We can appreciate a very clear peak and the endpoint structure around the pseudoscalar mass for both scenarios.
After choosing the pairing of the leptons, the transverse momentum p ℓ + ℓ − T ≡ |p ℓ + T + p ℓ − T | and the angular separation ∆R ℓ + ℓ − are calculated for each lepton pair. It is likely that the leptons coming from signal events are fairly boosted in the nearlycollinear direction. This can be observed in the right panel of Fig. 1, and provides a clear criterion for discriminating the signal and background. Motivated by this, we include selection cuts on the p ℓ + ℓ − T and ∆R ℓ + ℓ − as follows, An important observation is that after applying the above cuts, the h 0 2 → ZZ * → 2ℓ + 2ℓ − background is found to be almost negligible since the collinearity of the leptons is not valid for such kind of events. On top of that, the invariant masses of the oppositely-charged lepton pairs are required to be m ℓ + ℓ − < 20 GeV to ensure that the lepton pairs come from the light pseudoscalar. Collectively, a low ∆R ℓ + ℓ − together with a low m ℓ + ℓ − cut can discriminate the studied signal from the h 0 2 → ZZ * → 2ℓ + 2ℓ − process.
We also calculate the cluster transverse mass defined as where V denotes the four-lepton system [94]. The endpoint position of the M T distribution corresponds to the parent particle mass of the pseudoscalar bosons, i.e., m h 0 1 or m h 0 2 if there exists an event with vanishing invariant mass of the neutrinos and all the particle tracks are on the transverse plane. However, one cannot distinguish the decay event of h 0 1 from that of h 0 2 as they lead to the same final states. The M T distributions are shown in Fig. 3. For both scenarios, the signal distributions are largely populated around/below m h 0 1 , whereas no clear peak structure is observed around m h 0 2 . The main reason is the different production cross sections of h 0 1 and h 0 2 . At the generator level we obtain σ gg→h 0 1 →2a 0 1 →4ℓ /σ gg→h 0 2 →2a 0 1 →4ℓ ∼ 7, and this ratio is still as large as ∼ 4 in the detector-level data when applying the basic selection cuts to the triggered events in the case of the BP1 5 . For the same reason, the M T distribution of the BP2 has a clear peak around the m h 0 1 value as well. Although the contribution from the small h 0 2 decay events smears out the edge of the M T distribution, one can still estimate the light Higgs boson mass scale by the peak position.
The M T variable can also be used to define the signal region. We consider the events with 35 < M T < 145 GeV, which corresponds to the region where two light Higgs bosons can be discovered. This definition of the signal region might be useful if there exists a new physics process whose final state contains several leptons and massive invisible particles. Even though the new processes can contribute to the background of the Higgs signal, they can be readily isolated since the M T values of such processes are generically large.
For events passing all the selection cuts, we also construct the invariant mass of all visible leptons, m 2ℓ + 2ℓ − . The peak position of the invariant mass distribution is generically different from m h 0 1 and m h 0 2 due to invisible neutrinos from the tau decays. We impose a cut, requiring the invariant mass to be larger than 40 GeV, which accommodates to the possible spread of the invariant mass distribution. Another possible collider variable is the invariant mass reconstructed by the M T2 -assisted on-shell approximation of the invisible momenta [95], which has been found to be useful for reconstructing the resonance mass peak in the similar decay topology [96]. However, by definition, we would need to sacrifice a certain signal statistics by imposing a strong M T2 cut to attain a good accuracy. We expect that this variable would be useful with higher luminosity.
After imposing all the event selection cuts, we find that all the considered SM backgrounds turn out to be completely negligible. In order to estimate the signal significance, we show how the number of events of the Higgs signals and the backgrounds change under each event selection cut in Table 5. The cuts on p ℓ + ℓ − T and ∆R ℓ + ℓ − are especially efficient in discriminating the signal from the otherwise dominant DY and bb backgrounds, resulting in a good statistical significance for the benchmark point BP1 and a poorer result for BP2. We defineμ as the statistical sig-  nificance calculated with the Poisson probabilities since the number of background events is small when applying all the cuts. Indeed, the cut efficiency of each scenario is rather different. One reason seems to be the dissimilar efficiency of the object reconstructions. In BP1 the decay products of the pseudoscalar boson are relatively adjacent as it is likely that they are more boosted than in the case of BP2, and this can affect the correct identification of the final-state particles.
Notice finally that since there is no source of hard jets except for initial state radiation in the signal process, one can further impose a jet-veto condition that discards events containing hard jets and can be employed for suppressing rare processes like tt+ jets or W W + jets, which have not been included in this study. The (small) effect of this cut with p T > 50 GeV on the signal, is listed in the last row for use in studies that consider the other possible backgrounds.
To summarise, our simulation and analysis establish that the current 8 TeV data from the LHC, with dedicated selection cuts, can be used to test parts of the parameter space of non-standard Higgs scenarios with very light pseudoscalars.

Conclusions
In the light of the recent detection at the LHC of a scalar boson with a mass of 126 GeV and the possible hint for an excess in LEP that could correspond to a second, lighter scalar we investigate the phenomenology and detectability of viable models in the NMSSM. We consider two scenarios that contain two light Higgses, in which the mass of the lighter scalar is either in the range of 96 -100 GeV or lighter than half of the SM-like Higgs mass. Both of them are analyzed in the presence of a very light pseudoscalar boson, 2m τ < ∼ m a 0 1 < 2m b , a range which is yet to be explored further at the LHC. The Higgs phenomenology of these scenarios is extremely rich, involving cascade Higgs decays that can lead up to eight leptons in the final state, with or without missing energy. We perform a reconstruction of these channels for two representative benchmark points, determining the optimal cuts that allow us to single out the signal from the background. We show that the pseudoscalar mass can be successfully determined from the M T2 distribution of four visible leptons and the missing energy. We also investigate the reconstruction of the light scalar Higgs boson from the M T of the four visible leptons and the missing energy. The Higgs boson mass can be estimated by the peak position of the M T distribution. Our analysis suggests that the experimental search for inclusive decay modes using the current 8 TeV LHC data and dedicated cuts can be used to test parts of the parameter space on these possible two-Higgs scenarios with very light pseudoscalars.