Review of LHC experimental results on low mass bosons in multi Higgs models

A number of searches at the LHC looking for low mass ($2m_{\mu} - 62\ \mathrm{GeV}$) bosons in $\sqrt{s} = 8\ \mathrm{TeV}$ data have recently been published. We summarise the most pertinent ones, and look at how their limits affect a variety of supersymmetric and non-supersymmetric models which can give rise to such light bosons: the 2HDM (Types I and II), the NMSSM, and the nMSSM.


Introduction and Motivation
Since the discovery of a Higgs boson in July 2012 by the ATLAS and CMS collaborations at the Large Hadron Collider (LHC) [1,2] innumerable analyses have been performed in order to ascertain its nature. While its profile is largely consistent with the predictions of the Standard Model (SM), there remains the possibility that this object belongs to a Beyond the SM (BSM) scenario in which a SM-like Higgs state is realised in specific configurations of the corresponding parameter space. Since the necessity of BSM physics is evident from both the theoretical (hierarchy problem, absence of coupling unification, etc.) and experimental (neutrino masses, dark matter, etc.) point of view, it is of paramount importance to investigate whether it is possible to access it through Higgs analyses.
A possibility is clearly to improve the precision of the measurements of the discovered SM-like objects as, sooner or later, statistically significant deviations from the SM predictions may well appear. It should be emphasised, however, that accessing BSM physics indirectly, i.e., through the study of SM-like production and decay channels of the 125 GeV Higgs boson, may not be the most efficient way of isolating the underlying BSM scenario.
An alternative procedure is the following one. Whichever the BSM scenario encompassing the discovered SM-like Higgs state, this obviously includes an extended Higgs sector, with respect to the SM, hence a Higgs mechanism of Electroweak Symmetry Breaking (EWSB) giving rise to more physical Higgs states than just the single one of the SM. Crucially, other than with SM states, all these emerging Higgs boson (both neutral and charged, both scalar and pseudoscalar), can interact with each other. For example, the heavier Higgs states can decay into the lighter ones and in these chains the 125 GeV Higgs boson could, if appearing, either be the initiator or else the end product of the various possible decay patterns. Needless to say, to isolate one or more of the latter would be a direct evidence of a non-SM Higgs sector, hence of the existence of BSM physics. Furthermore, the study of the additional Higgs states would certainly gain one much more understanding of the underlying scenario than what can be extracted from the aforementioned analyses of the SM-like Higgs state.
It is the purpose of this paper to review both the theoretical and experimental status of several BSM scenarios predicting such Higgs cascade decays, in particular, those embedding in their particle spectrum a rather light state, with mass below 60 GeV or so, which would be produced in pairs in the last step of the discussed Higgs cascade decays. From the theoretical side, we will concentrate on the most popular BSM Higgs scenarios in which such a light object is realised, which is typically pseudoscalar in nature. From the experimental side, we will adopt published data obtained by the end of Run 1 of the LHC from either ATLAS and CMS, covering several signatures of such a pair of pseudoscalar Higgs states, including decays into pairs of muons, taus, and bottom quarks.
It is natural to organise the discussion of the possible BSM scenarios behind such a decay phenomenology around the divide of BSMs with and without Supersymmetry (SUSY). In fact, among the possible BSM theories, SUSY remains one of the favourite ones. However, while its minimal realisation, the Minimal Supersymmetric Standard Model (MSSM), has been under close experimental scrutiny lately, through direct searches for both its sparticle and Higgs states, much less effort has gone into testing non-minimal SUSY scenarios. Amongst the latter, a particular role is played by the Next-to-MSSM (NMSSM). Further, a slight variation of the latter, known as the New Minimal Supersymmetric Standard Model (nMSSM), has recently also undergone significant phenomenological scrutiny. All such SUSY scenarios are built upon a Higgs sector which is essentially one particular realisation of a 2-Higgs Doublet Model (2HDM), with (NMSSM, nMSSM) or without (MSSM) an additional Higgs singlet field. Thus, if one abandons the paradigm of SUSY, it is natural the examine generic 2HDMs. In fact, all such extended Higgs models are capable of producing the Higgs cascade decays which are of interest here, apart from the MSSM, which we will then not test. Regarding the others, we will tackle them in turn.
This paper is thus organised as follows. In the next section, we shall review the discussed theoretical models (in separate subsections) while in the following one we will describe the experimental analyses exploiting the mentioned signatures. Our results, obtained by confronting predictions from the former with constraints from the latter, are presented in Section 4. Finally, we conclude in Section 5.

Models
We now briefly review a number of models, all based around the Two Higgs Doublet Model (2HDM). Whilst they differ in their input parameters and number of fields, they all share the ability to produce small mass (pseudo)scalars, with sizeable Higgs-to-Higgs couplings.
A scan over parameter space was performed for each model, targetting scenarios with low mass a 1 or h 1 . Scans were subjected to many existing experimental constraints. SM Higgs searches and measurements can be used to place indirect limits on our models. All scans used HiggsBounds 4.3.1 [3][4][5][6][7] to implement current Higgs exclusion limits. HiggsSignals 1.4.0 [8] was used to apply measured Higgs signal rate constraints in a variety of channels. This was run in peak-centered mode, with a Gaussian probability distribution, and requiring the overall p > 0.05. For the NMSSM and nMSSM scans, we also consider the individual requirements on ZZ/γγ/bb Higgs signal rates which can result in different exclusion regions. This will be discussed further in Section 4.
We also consider non-Higgs constraints, including those on flavour variables, the anomalous muon magnetic moment a µ , and dark matter (DM) relic density Ω DM h 2 . All scans use micrOMEGAs [9] to implement the latter constraint. We apply a "relaxed" set of constraints, requiring points to pass all constraints but allowing any ∆a µ > 0 and Ω DM h 2 < 0.131, and ignoring constraints on R(D), R(D * ) [10,11]. This allows for future developments and changes in those calculations, whilst still accommodating some BSM contribution. This does not significantly modify the results in Section 4.

Type I and II Two Higgs Doublet Model (2HDM)
The Two Higgs Doublet Model (2HDM) represents one of the most economical extension of the SM Higgs sector, providing a simple, yet comprehensive, framework for studying extended patterns of EW symmetry breaking. In the 2HDM a second complex Higgs doublet with the same quantum numbers of the SM one is added to the SM Higgs sector. The scalar spectrum of the 2HDM is thus enlarged to include two CP even states, denoted as h and H (with m h < m H ), a CP odd state, A, and a charged Higgs, H ± . In general, the role of the SM like Higgs boson can be played by either h or H.
Denoting the two Higgs doublets as Φ 1,2 , the most generic scalar potential of the 2HDM that respects a Z 2 symmetry distinguishing Φ 1 and Φ 2 can be expressed as [12] The imposition of a Z 2 symmetry, together with the assignment to the right handed SM quarks of a defined Z 2 quantum number, is necessary so as to avoid Higgs mediated flavour changing neutral currents (FCNC). Note that the Z 2 breaking term m 2 12 is generally tolerated, since it breaks the Z 2 symmetry softly, i.e. the symmetry is restored in the UV [13].
In the potential of eq. (2.1) the parameters λ 1−4 , m 2 11 and m 2 22 are real numbers, while m 12 and λ 5 are in principle allowed to be complex valued numbers. However, complex parameters that cannot be made real through a suitable transformation give rise to CP violation in the Higgs sector. Since we are not interesting in the study of these effects, in the following we will consider all the parameters of the potential to be real numbers.
Starting from the scalar potential of eq. (2.1), various 2HDM realisations can then be formulated according on how the SM fermions couple to the two Higgs doublets. In particular we will focus in our analysis on the so called Type I and Type II 2HDMs. In Type I 2HDM all the SM fermions, up and down type quarks and down type leptons, couple to only one doublet while in Type II down type quarks and leptons couple to one doublet and up type quarks to the other doublet.
In order to scan the 2HDM parameter space we have used the package 2HMDC [14] with input parameters defined in the mass basis. In this basis the free model parameters are the physical masses of the four scalar states (m h , m H , m A , m H ± ), the ratio of the two doublets vacuum expectation values (tan β = v 2 /v 2 ), m 2 12 , and sin(β − α), with α the mixing angle between the two scalar states. The parameter ranges used for the scan are indicated in Table 1. The 2HMDC package imposes basic theoretical constraints, such as stability of the potential, tree level unitarity, and consistency with the S, T, and U EW parameters. Finally superiso [15] was used to check compatibility with current flavour constraints. However failing points were not explicitly excluded to increase the overall scan efficiency.

Next-to-Minimal Supersymmetric Standard Model (NMSSM)
The Next-to-Minimal Supersymmetric Standard Model (NMSSM) [16] is a simple extension of the MSSM, which adds a singlet S to its superpotential. Originally proposed to solve the µ-problem of the MSSM, the NMSSM has gained renewed interest as additional tree-level contributions to 10 -10 5 GeV 2 | cos(β − α)| 0.9 -1 the Higgs mass alleviates the need for large loop contributions to achieve its measured value, thus possibly allowing a more natural sparticle spectrum [17][18][19][20][21][22].
The inclusion of a new singlet scalar naturally also leads to more physical scalar particles: one scalar and one pseudoscalar will be added giving in total three scalars (h 1,2,3 ), two pseudoscalars (a 1,2 ), and the usual charged Higgs h ± . A novel feature is that the discovered Higgs can be assigned to either h 1 or h 2 . The latter possibility was found to be excluded in the MSSM by [23,24] due to a combination of flavour observables and LHC searches for scalars decaying to τ τ pairs, though one might add that more recently [25] claims there still is a very constrained possibility that the heavier scalar is the discovered one in the phenomenological MSSM.
The inclusion of the extra singlet superfield results in a modified superpotential, where λ and κ are dimensionless coupling constants, and we have assumed a Z 3 invariant model. The rest of the superpotential is formed from the usual Yukawa terms for quarks and leptons as in the MSSM. Further, one needs to add the corresponding soft supersymmetry breaking terms in the scalar potential, where m S , A λ and A κ are dimensionful mass and trilinear parameters, and one also has the other usual MSSM soft SUSY breaking terms. As the masses of the singlet dominated scalar and pseudoscalar are essentially free parameters, it opens the possibility for them to be very light. If the singlet component of a 1 is large enough, then such light particles can easily escape all exclusion limits from earlier searches. We briefly consider m a1 as a function of selected input parameters, showing the results in Fig. 1. Scan details are explained below. Relaxed constraints have been applied, apart from those on Higgs signal rates. Each horizontal bin is normalised such that the largest bin in each row has contents = 1. This allows one to see which value(s) of input parameter are preferred for a given m a1 . There are a few salient features to note. Most strikingly, panel (a) shows that A κ ∼ 0 or slightly negative is highly favoured for a light a 1 scenario. Panel (b) indicates some preference for κ 0.3, with another "hotspot" of points at κ ∼ 0.02 − 0.04. Panel (c) also shows a weak preference for a fairly small λ ∼ 0.15.
Whilst a scalar with mass ∼ 125 GeV is easily achievable in the NMSSM, it is useful to momentarily review its dependence on the model input parameters. A scalar with mass 125±3 GeV is achievable over the parameter range scanned. Fig. 2 shows the dependence of m h1 on selected   Each horizontal bin is normalised such that the largest bin in each row has contents = 1. Relaxed constraints have been applied, apart from those on Higgs signal rates.
parameters where there are noticeable trends. Relaxed constraints have been applied, apart from those on Higgs signal rates. In particular, A t (left panel) sets an upper limit on m h1 through its effect on the stop mixing which in turn effects the loop contributions to the Higgs mass. Additionally, smaller values of λ (central panel) tend to push m h1 to larger values. It may seem surprising that smaller λ allows larger m h1 , while the NMSSM specific contribution to m h1 is proportional to λ. But in our case all large λ are already excluded by the signal rate constraints and m h1 only shows a clear growth with λ for λ>0.4, below that one also has to remember that λ affects the mixing of the scalars and thus can have a more complicated impact on m h1 . We also see that smaller values of κ∼ 0.1 − 0.3 (right panel) are preferred in order to satisfy signal rate constraints for h 1 .
There have been numerous studies of light pseudoscalars in the NMSSM and their discovery prospects, see, e.g., [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43] but the present study is the first attempt to investigate the impact on the NMSSM parameter space from LHC searches for light pseudoscalars. For our analysis, we have performed scans for both the Z 3 -invariant NMSSM (hereafter referred to as just the NMSSM), and a GUT inspired NMSSM. In the latter, one has a common parameter for all scalar masses (m 0 ), a common parameter for all trilinear parameters except A λ and A κ (A 0 ), and a typical GUT relation between the gaugino masses (M 2 = M 1 /2 = M 3 /3 = m 1/2 at the EW scale). The singlet pseudoscalar mass parameter, M p , is used as an input parameter in the GUT scan instead of A κ , requiring input parameters to be specified at the EW scale to be effective. The parameter ranges for the NMSSM scan are described in Table 2, while the ranges in the GUT inspired scan are given in Table 3; here two scans were made, one (reduced range) focusing on the region with large λ and small tan β to optimise the NMSSM specific contribution to the Higgs mass, and one broader (extended range) to ensure no possibility was missed.
1000 GeV A e/µ/τ 2500 GeV Table 2: NMSSM parameters and their ranges used for the scans. All parameters are specified at the SUSY scale.
Parameter Extended range Reduced range Table 3: Parameter ranges used in the GUT inspired NMSSM scans. All parameters are specified at the EWK scale.
In order to use HiggsSignals with the output from NMSSMTools, we add a DMASS block to the SLHA file to represent theoretical uncertainties on the h 125 mass. This is set to 2 GeV for both h 1 and h 2 . Additionally, HiggsSignals was modified to ensure that either h 1 or h 2 was correctly assigned to h 125 by increasing assignmentrange_massobs to 2.0 in usefulbits_HS.f90.

New Minimal Supersymmetric Standard Model (nMSSM)
In the previous section we have described the properties of the Z 3 invariant NMSSM. However, a general 2HDM+S superpotential might not posses this accidental symmetry. A different realisation, called the new Minimal Supersymmetric Standard Model (nMSSM), possesses instead a discrete Rsymmetry that forbids a cubic singlet term in the superpotential but allows for tadpole terms. While the field content of the nMSSM is the same as that of the Z 3 invariant NMSSM, the phenomenology can be quite different due to the different superpotential and soft SUSY breaking terms.
The first striking feature of the nMSSM is the absence of a mass term for the pure singlino, whose mass can be raised up to ∼ 75 GeV only via mixing effects. The singlino is thus naturally light and the LSP, which generally contains a large singlino component, can have a mass lighter than ∼ 5 GeV, leading to a quite different phenomenology for the nMSSM in both collider and DM searches.
The Higgs sector of the nMSSM superpotential reads [16] (in contrast to eq. 2.2) to which the usual Yukawa terms are added. The corresponding soft SUSY breaking terms are very similar to eq. 2.3, but removing the κ 3 A κ S 3 term and introducing a tadpole term: where ξ F and ξ S are O(M 2 SU SY ) and O(M 3 SU SY ) terms which avoid domains walls and stability problems of the nMSSM (see [16]).
Our reinterpretation of the constraints arising from low mass 8 TeV scalar searches will be based on the results presented in a recent paper [50] that reviews the status of the nMSSM after the first run of the LHC and highlights the prospects for this model for the 13 TeV run of the CERN machine. Referring to [50] for more details, we summarise here the major details of the parameter scan and of the constraints imposed. NMSSMTools has been used to scan over the following parameters: all defined at the GUT scale except tan β, defined at M Z , and λ and µ, both defined at the SUSY scale. We impose the following universal soft terms conditions at the GUT scale: Regions of the parameter space where sparticles are out of the LHC reach have been discarded, thus only focusing on regions with interesting prospects at present and future colliders. Constraints on direct sparticle searches at LEP, Tevatron, and the LHC have been implemented via the SModelS [51,52] and MadAnalysis5 [53][54][55] packages.
In [50], three regions compatible with the aforementioned combination of theoretical, cosmological and collider constraints were identified. In two of them the LSP has a mass of ∼ 45 GeV and ∼ 70 GeV respectively, while a third region features a light LSP, m LSP < 5 GeV. This is the only region with a light spin 0 state, a 1 , in the mass range of interest for this paper. In particular one has m a1 ∼ 2mχ0 1 , which ensures an efficient annihilation in the early Universe and thus provides a relic abundance compatible 1 with the value measured by the Planck collaboration [56]. Within this region, there are two different subregions, denoted as 1A and 1B. Region 1A is characterised by a small m 0 and M 1/2 , both below 1 TeV, whilst region 1B has a small M 1/2 (< 500 GeV) and large m 0 (> 4 TeV). Their full parameter ranges are reported in Table 4.
Unlike the NMSSM, in both these regions the role of the SM Higgs boson is played by h 2 , with h 1 having a mass between 35 and 70 GeV. As previously mentioned, a 1 is the lightest of the Higgs states which has a dominant singlino component, while the remaining heavier Higgs are decoupled. In particular region 1B features an extremely light gluino, with mg 1.2 TeV, and is almost nearly excluded by run 1 searches. LHC results for stop and slepton searches also strongly constrains region 1A, via , which are light in this part of the parameter space where m 0 is small.

New Experimental Analyses
There are several recent experimental analyses searching for light bosons which may impinge on the parameter space of the aforementioned 2HDM and NMSSM/nMSSM scenarios. We provide an overview of the ones most relevant to this investigation, categorised by their final state. Note that while we refer to a 1 , it should be understood that this can refer to a generic light boson, a 1 or h 1 . For scenarios where m a1 << m h , a common theme is that of "boosted" topologies, where the a 1 is significantly boosted, and therefore its decay products are highly collimated [57]. The separation is of the order ∆R ∼ 2m a1 /p a T ∼ 4m a1 /m h , where we have assumed that each a 1 has a transverse momentum p a T ∼ m h /2. For m a1 ∼ 8 GeV, we therefore expect ∆R ∼ 0.3. Analyses must therefore take care to ensure standard isolation criteria do not inadvertently quash any potential signal. At larger m a1 , the a 1 is no longer highly boosted, and there is good separation between its decay products. Standard reconstruction techniques can therefore be used. The intermediate region, m a1 ∼ 15-20 GeV, proves the most challenging since the decay objects are neither neatly collimated, nor well separated.

Adapting Experimental Limits
One can adapt the limit from a search for one final state to place a limit on another, given a relationship between the corresponding final states. The channel widths are given in [58]. Since all leptons and down-type quarks couple to the same doublet in the models under consideration, there is no tan β dependence and the conversion is simple. For µµ → τ τ : is the velocity factor. For bb → τ τ : where the radiative corrections are where N f is the number of active light quarks;ᾱ s is the running strong coupling constant;m q is the running quark mass in the MS scheme; and α is the QED coupling constant. The running parameters are evaluated at scale µ = m a1 using [57,[59][60][61][62].

4τ
For the mass region 2m τ − 2m b , BR(a 1 → τ τ ) is expected to dominate in a Type II scenario with tan β 2. Ditau (or pairs of ditau) final states are therefore a natural search channel. However due to the nature of the tau decay, it can be a difficult object to fully reconstruct in a boosted regime. Taus can decay into 1, 3, or 5 charged particles ("prongs") along with one or more neutral particles, including neutrinos. The 1-prong and 3-prong decays modes make up ∼ 85% and ∼ 15%, respectively, of all tau decays. The multi-particle nature of the decay reduces the visible energy, making passing trigger thresholds and reconstruction more difficult then, e.g. , a 1 → µµ.
The CMS collaboration has published two analyses that search for 4τ final states arising from pairs of low-mass boson decays [63,64]. Whilst both look for h 125 → 2a 1 → 4τ , and cover similar m a1 ranges, they utilise different analysis strategies to identify the boosted tau pairs. Both analyses capitalise on the excellent muon reconstruction and low fake rates, and require two muons in an event.
The approach taken in [63] (CMS HIG-14-019) targets the tau 1-prong and muon decay modes. Ditau pairs are selected by looking for a well-isolated muon with only one nearby track with p T > 2.5 GeV. This forms a µ-track pair, and events are required to have 2 such pairs that are well separated. Backgrounds are almost entirely from QCD heavy-flavour decays, since Drell-Yan, tt, and diboson events are rejected by a same-sign requirements on the two muons. The µ-track invariant mass, m µ-trk , is used as the discriminating variable. A background template is formed from a QCD-rich sideband region, and fitted to the data along with signal template from MC to extract the size of any potential signal. Upper limits on the total σ × BR range from 10.3 pb at m a1 = 5 GeV down to 4.5 pb at m a1 = 8 GeV.
A complementary approach is taken in [64] (CMS HIG-14-022). This analysis targets both the gluon fusion and WH production modes. To target the boosted ditau pair, the standard tau reconstruction is modified. The tau reconstruction is seeded by anti-k T (with a 0.5 cone radius) [65] jet candidates. Candidate jets must have at least one muon constituent, which is removed before passing the remaining jet constituents to the tau reconstruction algorithm. This tau must have p T > 20 GeV and also pass isolation criteria. Events are required to have at least one such muontau pair. There is also an additional muon requirement, which must be well separated from the muon-tau pair. This is designed to be sensitive to W (µν)H production, or a muon from the other ditau pair in the gluon fusion and VBF production modes. The analysis uses the µ-τ invariant mass to define a signal region, only considering events with m µ-τ > 4 GeV. Upper limits on the total σ × BR range from ∼ 500 pb at m a1 = 5 GeV to 3.5 pb at m a1 = 11 GeV.
Both analyses are less powerful at smaller m a1 as a consequence of using the effective ditau invariant mass as the discriminating variable. Background events are characterised by small invariant mass, and thus there is a much larger overlap with a smaller m a1 signal, thereby reducing its discriminating power. In the case of HIG-14-022, the lack of any information below 4 GeV has a severe impact on the limit at small masses. Additionally, the use of the visible ditau invariant mass means that there is no longer a clean, sharp peak on a continuous background, reducing the sensitivity of the searches compared to a fully reconstructible final state e.g. µµ.

2τ 2µ
This final state is a compromise between the large but less clean τ τ final state, and the much cleaner but rarer µµ final state. CMS and ATLAS have both published results looking for a 2τ 2µ final state produced by light bosons [66,67]. Both analyses look for resonances in the dimuon invariant mass distribution, and are triggered by an asymmetric dimuon requirement with similar p T thresholds. The CMS analysis targets a mass range m a1 = [20, 62.5] GeV, whilst the ATLAS result covers a range m a1 = [3.5, 50] GeV, optimising for m a1 = 5 GeV. The two analyses are therefore complementary.
Since the CMS analysis targets much larger values of m a1 , the dimuon and ditau pairs will not be heavily boosted. Therefore the standard hadronic tau reconstruction algorithm and isolation requirements can be used. All four objects are required to be well-separated, and events with additional isolated leptons or b-tagged jets are vetoed. Requirements on the 4-body invariant mass and dimu-ditau mass difference are used to further enhance background rejection. Both the reducible background (from jets faking leptons), and the irreducible background (from ZZ → 4 ), are modelled by Bernstein polynomials. An upper limit on the 4τ cross-section is set, ranging from ∼ 2 pb at m a1 ∼ 20 GeV to ∼ 0.8 pb at m a1 ∼ 60 GeV.
In contrast, since the ATLAS analysis optimised for a much smaller mass, the kinematic and topological regime changes. The dimuon and ditau pairs will now be heavily boosted, and akin to the CMS 4τ analysis the ditau selection criteria avoids the use of a standard tau reconstruction algorithm, instead opting for a µ/e + tracks requirement. The dimuon requirements include an isolation requirement, which is modified to remove the other muon. This improves sensitivity at low m a1 at the expense of reduced sensitivity at higher m a1 . Due to the mass range, the background estimation must now take into account various quarkonia resonances, as well as contributions from a continuum Drell-Yan background at smaller m a1 , and tt at large m a1 . The final upper limit on the 4τ cross-section extends down to < 1 pb for m a1 ∼ 4 GeV, but worsens at higher m a1 , where it only reaches ∼ 20 − 30 pb. Since the selection criteria are not adapted for larger m a1 , this is to be expected.
Interestingly, the ATLAS limit is better at smaller m a1 despite the increase from the Drell-Yan background at smaller m µµ . This is due to an increased signal efficiency. The lighter a 1 receives a larger boost and therefore has a larger p T on average, ensuring that more muons and tracks pass the trigger and selection requirements. Whilst the same is also true in the 4τ analyses, in those analyses the increase in signal efficiency is not sufficient to overcome the propinquity for background to lie at lower invariant masses.

4µ
The region m a1 < 2m µ sees a large increase in BR(a 1 → µµ). Whilst not as large as BR(a 1 → ss, gg), the dimuon final state is very clean with small systematic uncertainties. Note that the other non-coloured final state, γγ, is still several orders of magnitude smaller than µµ. CMS has searched for a 4µ final state [68], targetting the pair production of very light (pseudo)scalars m a1 = [0.25, 3.55] GeV, each decaying to a pair of muons. This analysis searches for two dimuon systems, with invariant masses compatible within detector resolution. The muon pairing criteria takes into account situations in which the two muons are nearly parallel. To reduce backgrounds from heavy-flavour decays, a modified muon isolation requirement is used, in which the other muon in the pair is excluded from the isolation sum. The upper limits on the equivalent total 4τ crosssection is ∼ 0.7 − 0.9 fb.

2b2µ
Focussing on higher masses, once the 2m b threshold has been surpassed then this now becomes the dominant decay channel in the models under consideration (assuming tan β 2 the for Type II models). However a 4b search would have to overcome significant QCD backgrounds 2 . Instead, requiring one a 1 to decay to µµ would allow one to use m µµ as a powerful signal/background discriminant, improving search sensitivity. CMS has performed a search for h → 2a 1 → 2b2µ (HIG-14-041) [70], covering a mass range 25 − 65 GeV. In this mass range the a 1 is no longer boosted, and one can therefore utilise standard particle reconstructions algorithms. This analysis required events to have two isolated muons, along with two b-tagged jets, with the 4-body invariant mass close to 125 GeV. Signal and background functional templates are fit to the m µµ distribution in data, where the background is dominantly Z/γ + jets. An upper limit is set, which is equivalent to a limit on the total 4τ cross-section from 40 fb to 100 pb, assuming the relationships given in Section 3.1. It should be noted that unlike other analyses, this limit is fairly constant with respect to m a1 .

Results
We now analyse how these new constraints affect the model parameter space by first considering the factors that influence the total cross-section, using the NMSSM as an example. The total production cross-section predicted by a given model, σ × BR, is decomposed as follows: • g 2 ggh is the squared reduced ggh coupling, with respect to the SM value (1 in the SM by definition) • BR(h → 2a 1 ) is the branching ratio of h to 2a 1 • BR(a 1 → 2X) is the branching ratio of a 1 to 2X where X = τ, µ, · · · • f is a combinatorics factor: 1 if the final states X and Y are identical, 2 otherwise.
Note that we only consider gluon-gluon fusion production, since it is the dominant production mechanism. There are several scenarios that involves light boson pair-production that we must consider: if h 1 = h 125 , then we could have We now consider the squared reduced gluon-gluon-Higgs coupling, g 2 ggh , which a priori is not constrained by the model. Instead, it is heavily constrained by current experimental results. If h i is assigned to be h 125 , then Higgs coupling measurements mean it must be SM-like, i.e. g 2 gghi ∼ 1. If however it is not h 125 , then current exclusion limits mean its production must be suppressed, i.e. have a small g 2 gghi . The ggh 1 squared reduced coupling g 2 ggh1 is shown in Fig. 3 as a function of several input parameters. Blue points indicate models where h 1 = h 125 , whilst orange diamonds are models where h 2 = h 125 . Relaxed constraints have been applied, along with those on HiggsSignals and HiggsBounds. We note that g 2 ggh1 is far larger in models where h 1 = h 125 compared with models where h 2 = h 125 . Additionally, in the former scenarios g 2 ggh1 is easily able to reach 1 across the whole range of parameters scanned, in the latter it is confined to certain region of parameter space: particularly small κ, and large A λ , with moderately sized λ. Generally, it is somewhat favoured to have h 1 = h 125 . g 2 ggh2 follows a similar pattern: when h 2 = h 125 the reduced coupling can reach 1, whilst it is much smaller when h 1 = h 125 . However, the former scenario is now confined to those aforementioned regions of parameter space: small κ, and large A λ , with moderately sized λ.
The Higgs-to-Higgs branching ratio can also take on a range of values, and is again only limited by current Higgs measurements. Fig. 4 shows heatmaps of BR(h 1 → a 1 a 1 ) against several model input parameters for points where h 1 = h 125 and m a1 < 60 GeV. No Higgs coupling constraints have been applied from either HiggsSignals or NMSSMTools, but all other constraints have been applied. Each plot is normalised such that each horizontal bin is scaled so that the largest bin in each row has contents 1. This allows us to determine the sensitivity of a given BR value against a model parameter. Without any Higgs signal constraints, the BR can take on any value. We can also see clear features that show significant dependence of BR(h 1 → a 1 a 1 ) on these parameters, particularly κ, λ, and A κ ; but also some slight dependence on tan β. The dependence on κ and λ can be understood due to the presence of λ 2 and κλ terms in the relevant coupling. Also A κ appears in that coupling, while the effect of tan β is more indirect as it changes the relative importance of the λ 2 and κλ terms.
Adding in current Higgs coupling constraints requires a SM-like scenario for the SM decay channels and therefore a small BR(h 125 → BSM), with the most recent combined fits from CMS and ATLAS constraining BR(h 125 → BSM) < 0.34 at 2σ [72]. A small BR therefore primarily relies on a small κ 0.3 − 0.4, a small λ 0.2 − 0.3, and a negligible or slightly negative A κ . There is also a preference for large tan β ∼ 10 − 25, and large A λ ∼ 3 TeV. Note that we have not considered h 125 → Za 1 decays, since their BR are typically 10 −8 .
Since we are interested in the product of the reduced coupling and BR, it is useful to plots their correlations. The gluon-gluon higgs reduced coupling g 2 ggh is shown in Fig. 5 plotted against BR(h → 2a 1 ) for all the above assignments. Overlaid are contours of constant g 2 ggh × BR(h i → a 1 a 1 ). Two version of this plot have been made: one (Fig 5a) for points passing the HiggsSignals and HiggsBounds constraints, ignoring the NMSSMTools χ 2 constraints; and and therefore potentially the largest total σ × BR. These points have a very SM-like ggh coupling as a result of meeting visible ZZ/γγ/bb signal rates, and are limited entirely by the experimental constraints on BR(h → a 1 a 1 ). Points where the heavier h in the decay chain is not the 125-like object have the opposite trend. Given the lack of any other observed Higgs boson, these must have a small ggh coupling, but are free to have sizeable BR(h → a 1 a 1 ). However their overall product is typically smaller, 0.05.
A noticeable difference between the two plots is the allowed BR(h → a 1 a 1 ), particularly in the h 1 = h 125 scenario where HiggsSignals + HiggsBounds allows BR 0.5, whilst NMSSMTools constraints this more severely to BR 0.2. Note that the aforementioned combined result from CMS and ATLAS falls halfway between these two values. This is due to the differences between the programs: the experimental results they choose to use, and the manner in which they ap-   HiggsSignals, for models passing each NMSSMTools χ 2 constraint individually, and models passing all NMSSMTools χ 2 constraints. NMSSMTools performs a best-fit to each of the ZZ/γγ/bb final states as described in [73], and compares the model compatibility by calculating χ 2 for each final state. Therefore if at least one of those fails, the point will be rejected. We find that the ZZ χ 2 constraint places the strongest constraint on BR(h → a 1 a 1 ). However, NMSSMTools does not use information from other channels, such as τ τ . HiggsSignals in contrast uses information from a much larger set of analyses (85 in version 1.4.0), and performs a global χ 2 fit. Therefore, one can have a large rate in a certain channel if it is compensated by a low rate in another channel.
The last piece of the eq. 4.1, BR(a 1 → 2X), is shown in Fig. 7. For each final state, at a given a 1 , there is little variance over the range of input parameters. This is a consequence of all decays depending on the same Yukawa couplings; the total width of the a 1 can vary depending on its mixing and the value of tan β, but in the branching ratios all this is factored out and we are left with functions of a few parameters in the Higgs sector that are already fixed by phenomenology. The few points that deviate from the lines in Fig. 7 can be understood from the occasional presence of other channels, e.g. a 1 → γγ that is sometimes enhanced by large chargino loops.
Since the branching ratios are dependent on the Yukawa couplings, we see that it is the heaviest decay products that dominate and this manifests in the boundaries at (∼ 3.5, 10.5 GeV) where heavier final states (τ τ, bb) become kinematically viable, above which they become the favoured decay channels. Note that for τ τ this threshold happens at 2m τ as expected, while for bb the threshold is set to twice the B meson mass, which is somewhat larger than twice the b-quark mass (in principle there could be decays including mesons with b quarks also just below this limit, but the calculation of such channels is very challenging and not included in NMSSMTools).
One striking feature of Fig. 7 is the behavior of BR(a 1 → gg), which in the mass window 3.5  -10.5 GeV is dominated by the b-quark loop. The contribution from this loop increases rapidly until m a1 reaches ∼ 9 GeV at which point the quarks in the loop become real, after which it slowly decreases (due to increasing virtuality of the quarks). This threshold does not coincide with the onset of the bb channel since the loop behaviour is governed by the b-quark pole mass, and not the B meson mass that governs the threshold. This behavior is replicated in BR(a 1 → cc) due to this channel being dominated by a 1 → gg * → gcc, where g * is a virtual gluon. The kink in the BR(a 1 → gg) line at 9 GeV is also mirrored in the other branching ratios since a decreasing width to gluons will result in an increasing branching ratio for all other final states. We further know that the width of a channel typically increases quickly with the mass of the mother particle just above its kinematic threshold, then increases slower when the phase space factors become less dominant. This explains why, for example, BR(a 1 → τ τ ) increases in the region from 2m τ to around 6 GeV, and in turn explains the decrease in BR(a 1 → µµ) and BR(a 1 → ss) in the same region.
Below 2m τ , ss is the dominant decay channel due to its relatively large mass, as well as colour factors that favour quarks over µµ. BR(a 1 → ss) decreases due to the increasing gluon final state, while BR(a 1 → µµ) stays constant as the tendency to decrease due to increasing BR(a 1 → gg) is compensated by the fast increase in width due to being close to threshold. There are also QCD effects giving quark channels a flatter curve close to threshold as compared to leptons; this is why ss decreases while µµ remains constant. This is also why BR(a 1 → τ τ ) is increasing slightly above 10 GeV; BR(a 1 → bb) increases somewhat slower than BR(a 1 → τ τ ) despite being closer to threshold.
From the above studies, we expect total 4τ cross-sections up to ∼ 19.3 × 0.2 × 0.9 2 3 pb if one applies the NMSSMTools Higgs signal rate constraints, or even up to ∼ 8 pb if one uses the HiggsSignals constraints. The experimental Higgs signal rate measurements are therefore the limiting factor in determining the total cross-section due to their impact on BR(h → 2a 1 ), and not any particular model feature. We now combine all these pieces together, and plot the total crosssection as a function of m a1 . We start by considering the 4τ final state in the NMSSM. This is shown in Fig. 8, where σ × BR(gg → h i → 2a 1 → 4τ ) has been plotted against m a1 , for masses greater than 2m τ , with different assignments for h 125 , and different Higgs signal rate requirements applied.
There are very few h 2 → 2h 1 → 4τ points, and these have not been shown due to significantly smaller cross-sections. Points are required to pass the "relaxed" set of constraints, where we also require all other NMSSMTools constraints but allow any ∆a µ > 0, Ω DM h 2 < 0.131, and ignore limits on R(D) and R(D * ). Requiring lower bounds on ∆a µ and Ω DM h 2 does not change the overall result, and only reduces the overall number of points. Overlaid are the observed exclusion limits from relevant searches. One can see a wide variety of predicted cross-sections compatible with current experimental constraints, ranging from < 1 fb up to 8 pb. As previously mentioned, models with h i = h 1 (of which many have h 1 = h 125 ) generally have a larger cross-section than those with h i = h 2 . The large decrease in cross-section for masses m a1 > 2m b ∼ 10.5 GeV is due to the decrease in BR(a 1 → τ τ ) as the bb final state becomes kinematically available. The ATLAS 2τ 2µ analysis is more powerful for masses 4 -10 GeV, especially at smaller masses, and is therefore complementary to the 4τ analyses which lose sensitivity at smaller masses. This analysis can exclude a significant number of points of h 1 → 2a 1 , excluding cross-sections as small as 1 -2 pb, even taking into account the more restrictive Higgs signal rate constraints from NMSSMTools. However, it is not yet sensitive enough to probe the alternate scenario where h 2 → 2a 1 . The 4τ analyses start to intrude on the model space, although only if one assumes the more relaxed rate constraint from HiggsSignals. These excluded points are typically those where h i = h 125 , as such configurations often give a larger cross-section as shown in Fig. 5.
A minor detail seen in Fig. 8a is that the rates can go slightly higher when h 1 = h 125 . This is a somewhat complicated effect from the structure of the parameter space; first, if λ is large it is difficult to achieve acceptable SM signal rates for h 125 if m a1 < m h /2, mostly because BR(h 125 → a 1 a 1 ) tends to increase with λ, but also due to interplay with λ affecting the mixing of the h 125 . Furthermore, if λ is not too large we can only have h 2 = h 125 if κ is also small (the singlet scalar mass goes as κs and since λs cannot be too small κ > λ means the singlet scalar is heavy). The coupling h 125 a 1 a 1 has a term proportional to λκ which can saturate BR(h 125 → a 1 a 1 ) with respect to the h 125 signal rates if κ is large. Hence the rate shown in Fig. 8a can reach its maximum for h 1 = h 125 but struggles to do so for h 2 = h 125 .
Expanding our mass range up to 60 GeV means the limits from the CMS 2τ 2µ and 2b2µ analyses can also be included. This is shown in Fig. 8b, where one can see that the latter analysis is powerful enough to start to probe phase space (if one uses the more optimistic constraints from HiggsSignals). The 2τ 2µ analysis is not yet able to probe NMSSM phase space. However it could offer some sensitivity if one were instead dealing with a model where a 1 → τ τ was enhanced over a 1 → bb, for example a Type III (IV) 2HDM with large (small) tan β. Crucially direct searches have similar, often better, sensitivity measuring BR(h 125 → BSM) than limits from indirect searches, assuming BSM is solely a 1 a 1 .
One additional point to note in this Figure is the lack of points with m a1 ∼ 4 − 4.5 GeV and m a1 ∼ 5 − 5.5 GeV. These masses are heavily suppressed due to flavour constraints: the former mass range is excluded by BR(B → X s µµ), whilst the latter range is excluded by B s,d → µµ.
One might also consider the cross-section as predicted in the GUT-constrained NMSSM, shown in Fig. 9 for the 4τ final state. This shows a very similar result to that in Fig. 8b, whilst the limiting factor remains that on Higgs-to-Higgs decays. This bound is in general easy to satisfy, and hence the upper limit on possible rates in the channels we have studied is essentially independent of model details such as GUT scale unification.
Let us also consider the 4µ final state. σ × BR(gg → h i → 2a 1 → 4µ) in shown in Fig. 10 as a function of m a1 . Colours and shape assignments are the same as for the 4τ figure. The relevant experimental limits now include the CMS 4µ search. This probes cross-sections down to 1 fb, and therefore excludes many model points. There are almost no points below m a1 < 2.5 GeV. Points with m a1 < 1 GeV are rejected on grounds that their decay widths are difficult to calculate accurately due to hadronisation effects and QCD effects, while points with 1 < m a1 < 2.5 GeV are rejected by constraints on B → X s µµ.
Since the total cross-section is driven by the limit on BR(h 1 → 2a 1 ), which in turn has a strong dependence on several input parameters (Fig. 4), one can look at the impact of these new limits on possible model parameter values. κ, λ, and A κ are of particular interest. Histograms of the distributions are shown in Fig. 11, where they have been divided into points surviving all new constraints (blue) and failing any of the new constraints (red) for models with m a1 < 10.5 GeV. Also shown is the ratio of failing to surviving points for each bin, and the global fraction of points failing. κ and λ show a clear trend that higher values are more likely to be excluded, which is expected as the h i a 1 a 1 coupling depends on λ 2 and λκ. Additionally, small positive values of A κ also show a similar trend. Although these new constraints do not place a hard limit on values of these parameters, as such limits improve over time they will point towards models with smaller values. If experimental limits can exclude cross-sections down to 100s of fb then Fig. 12 shows that these parameters may be far more constrained, particularly κ due to the "knee" shape of its distribution, and A κ due its "wedge" shaped distribution. Constraining λ to smaller values is of particular interest since the tree-level Higgs mass has an additional contribution ∝ λ compared to the MSSM, one of the strengths of the NMSSM with respect to the MSSM. If this extra contribution is small, then a larger (and potentially more uncomfortable) degree of fine-tuning is required to achieve a mass of 125 GeV.
We now consider the results of scans for the other models. Although one might assume the nMSSM would give similar results to the NMSSM, the spectra of cross-sections and masses as shown in Fig. 13 is very different. From the Figure we observe that the a 1 mass (which we recall has a large singlino component) is constrained to be in a small mass window, between ∼5 and ∼11 GeV. As mentioned, this is due to the fact that the DM candidate (the lightest neutralinoχ 0 1 , which is almost a pure singlino) has a mass around 5 GeV and its relic abundance is fixed via annihilation through the lightest pseudoscalar a 1 . This constrains m a1 to be near the resonant peak, with 0 < m a1 − 2mχ0 1 < 1 GeV [50]. The singlino nature of the lightest pseudoscalar and the small mass ofχ 0 1 also makes BR(a 1 →χ 0 1χ 0 1 ) to be the dominant decay channel for the lightest pseudoscalar, therefore causing a reduction of the a 1 → τ τ rates and hence of the 4τ cross sections. Finally, we mention that in Fig. 13 we have included all points surviving the scan of [50]. However the values of m 0 and M 1/2 have a strong impact on the particle spectrum of the model. In particular, the region with small M 1/2 features a light gluino which is on the edge of the exclusion from 8 TeV searches that will soon be tested by the current run of the LHC. A similar consideration can be made for the region with small m 0 , that features light scalar superpartners (especially stops and sleptons). In this respect, the results of Fig. 13 has to be intended as to show only the current reach of light scalar searches in a different supersymmetric scenario, thus neglecting information arising from other LHC searches.
Lastly, we return to the more general Type I and II 2HDMs. Shown in Fig. 14 is the result of those scans for the 4τ final state. Both possible assignments for h 125 are shown. The Type II models predict significantly larger cross-sections (∼ 7 -8 pb) than in the Type I (∼ 1 pb), due to a different tan β dependence of the light pseudoscalar couplings, which favours higher BR(a 1 → τ τ ) in type II with respect to type I. In the Type I model, there are also fewer points with H = h 125 , whereas in the Type II model there is no such favouritism. In the Type II model, the cross-section range is similar to that in the NMSSM, since the limiting factor is the experimental constraint on BR(h/H → AA). Overall, we find that in both configurations of 2HDM Yukawas, current searches targeting light (pseudo)scalars are starting to scratch the edge of the predicted models cross sections thus making the LHC Run 2 a crucial probe also for these scenarios.  Observed exclusion limits ( s = 8 TeV)    Observed exclusion limits ( s = 8 TeV )

Conclusion
In summary, following the end of Run 1 at the LHC, we have assessed the status of direct searches for a light neutral Higgs boson in popular BSM scenarios with two Higgs doublets in both non-SUSY (2HDMs Type-I and II) and SUSY (NMSSM and nMSSM) frameworks, the latter also including an additional Higgs singlet field. The ability to extract signals of such a particle state would not only be a proof of a non-SM Higgs sector but also a circumstantial evidence of either a non-minimal SUSY (as such a signal is not available in the MSSM) or a non-SUSY scenario. The mass region concerned is up to 50 GeV or so. In such a range, the accessible decays, depending on the actual value of the light Higgs boson mass, are µ + µ − , τ + τ − and bb. The topologies searched for exploit a cascade chain wherein such a light Higgs state is produced in pairs from the decay of another Higgs state, where the latter could be the SM-like Higgs boson discovered in 2012 at the LHC or not. Hence, final state topologies are a combination of two amongst the aforementioned two-particle decays. Those pursued experimentally during Run 1, covering the discussed mass interval, were 4τ , 2τ 2µ, 4µ and 2b2µ. We exploited public results produced by ATLAS and CMS for these final states in order to set limits on the parameter space of all four scenarios considered, 2HDMs Type-I and II plus NMSSM and nMSSM. In doing so we have employed different numerical tools implementing these theoretical scenarios and/or corresponding experimental constraints, so as to enable us to distinguish genuine physics differences in the scope afforded by the various channels from artifacts due to the different degrees of accuracy in the model implementation.
Needless to say, the yield of these channels is not currently available in public tools, nor is the dedicated recasting procedure from one signature to another and onto a particular theoretical model that we have pursued here, so that our study represents an advancement in relation to current phenomenological knowledge, as the latter primarily rely on the study of SM-like signatures of additional Higgs states. Specifically, we have established that combinations of such signatures exclude substantial regions of the 2HDM Type-II (typically for masses below 10 GeV) but not in Type-I, which remains essentially untouched. As for the NMSSM and nMSSM, again, only one of these two scenarios is currently probed over significant portions of its parameter space (NMSSM), over the same mass range, while the other (nMSSM) is largely unaffected. Furthermore, the experimental searches considered do not make any assumption on the nature (whether scalar or pseudoscalar) of the light Higgs states, hence our results are applicable to whichever Higgs-to-two-Higgs decay pattern. We finally remark that all available experimental constraints were implemented, stemming from collider searches, both past (from LEP/SLC and Tevatron) and current (from LHC Run 1 and 2) ones, as well as from flavour and DM probes.
An obvious outlook of our work is to extend our analysis to forthcoming LHC Run 2 results for these and similar topologies, wherein we expect a substantially increased experimental sensitivity to the theoretical models considered.