Possible indications for new Higgs bosons in the reach of the LHC: N2HDM and NMSSM interpretations

In several searches for additional Higgs bosons at the LHC, in particular in a CMS search exploring decays to pairs of top quarks, $t \bar t$, and in an ATLAS search studying tau leptons, $\tau^+\tau^-$, local excesses of about $3\,\sigma$ standard deviations or above have been observed at a mass scale of approximately $ 400$GeV. We investigate to what extent a possible signal in these channels could be accommodated in the Next-to-Two-Higgs-Doublet Model (N2HDM) or the Next-to Minimal Supersymmetric Standard Model (NMSSM). In a second step we analyze whether such a model could be compatible with both a signal at around $400$GeV and $96$GeV, where the latter possibility is motivated by observed excesses in searches for the $b \bar b$ final state at LEP and the di-photon final state at CMS. The analysis for the N2HDM reveals that the observed excesses at $400GeV$ in the $t \bar t$ and $\tau^+\tau^-$ channels point towards different regions of the parameter space, while one such excess and an additional Higgs boson at around $96$GeV could simultaneously be accommodated. In the context of the NMSSM an experimental confirmation of a signal in the $t \bar t$ final state would favour the alignment-without-decoupling limit of the model, where the Higgs boson at $125$GeV could be essentially indistinguishable from the Higgs boson of the standard model. In contrast, a signal in the $\tau^+\tau^-$ channel can only be accommodated outside of this limit, and parts of the investigated parameter space could be probed with Higgs signal-rate measurements at the (HL-)LHC.


Introduction
In the year 2012 the ATLAS and CMS Collaborations have discovered a new particle with a mass of about 125 GeV [1][2][3]. Within the current experimental and theoretical uncertainties the properties of this particle, which in the following is denoted as h 125 , agree with the predictions for the Higgs boson of the Standard Model (SM) of particle physics, but they are also compatible with the interpretation as a Higgs boson in a variety of SM extensions corresponding to different underlying physics. Valid parameter regions in scenarios of physics beyond the SM (BSM) featuring extended Higgs-boson sectors are established taking into account the measurements of Higgs-boson couplings, which are known experimentally to a precision between 10% and 30% [4,5], the existing limits from searches for additional Higgs bosons, as well as other constraints (see the discussion below). Consequently, the question whether the observed scalar boson forms part of an extended Higgs sector is one of the science drivers for the upcoming run of the LHC beginning in 2022 (Run 3), and beyond.
The most obvious possibility for realisations of extended Higgs-boson sectors is that all additional Higgs bosons have masses that are larger than 125 GeV. But also cases where at least one of the additional Higgs bosons is lighter than the one at 125 GeV are phenomenologically viable. Thus, searches for BSM Higgs bosons ranging from very low to very high mass scales are crucial in this context. Among the models with extended Higgs-boson sectors the most studied ones are models with two Higgs doublets (2HDM), possibly with an additional (real or complex) Higgs singlet (N2HDM or 2HDMS), as well as their supersymmetric (SUSY) counter parts, the MSSM and the NMSSM.
The searches at the LHC for BSM Higgs bosons with masses above 125 GeV have shown several excesses in the data recorded between 2015 and 2018 (Run 2), which is a finding that by itself is not unexpected in view of the large number of searches that have been conducted. However, it is remarkable that several searches show an excess of events above the background expectation around the same mass scale of a hypothetical new Higgs boson φ of m φ ≈ 400 GeV. The various excesses each reach the level of about 3 σ local significance or above, but the global significance of any excess individually, taking into account the "look elsewhere effect", stays below 3 σ. The current experimental situation regarding the mentioned excesses at m φ ≈ 400 GeV can be briefly summarized as follows (we use here a notation in terms of the CP-even and -odd eigenstates, H and A, respectively, but interpretations in terms of CP-mixed states would also be possible).
-A → tt: CMS reported a local excess of 3.5 σ in their first year Run 2 data [6]. There is no corresponding Run 2 ATLAS analysis available. φ → τ + τ − : ATLAS reported a local excess of 2.7 σ in their full Run 2 data [7]. The corresponding CMS analysis, using only first year Run II data, does not show an excess, but also has substantially weaker expected sensitivities. -A → Zh 125 : ATLAS reported a local excess of 3.6 σ in their first year Run 2 data at a mass of around 440 GeV in the b-quark associated production channel [8]. The corresponding CMS analysis is consistent with the SM background expectation [9]. An updated full Run 2 ATLAS analysis for the gluon fusion (ggF) production channel also agrees with the SM expectation [10]. This experimental situation (see also the discussion in Ref. [11]) triggers the question whether all or some of these excesses could be fitted simultaneously. In a recent publication [12], in which this question was addressed in a model-independent approach by assuming a CP-odd state at 400 GeV with independently adjustable couplings to the SM particles as the origin of both excesses, it was found that the preferred values of the couplings would correspond to rather contrived scenarios that do not lend themselves to an immediate interpretation within UV-complete models. Within the present paper we take a different viewpoint and analyse the observed excesses in the context of two popular models, namely the N2HDM and the NMSSM. The symmetry properties of these models give rise to correlations between the couplings of the CP-odd Higgs boson to different SM particles, in contrast to the freely adjustable couplings that were considered in Ref. [12].
In our investigation in this paper we will perform a χ 2 analysis taking into account the observed excesses in the pp → A → tt channel at CMS and the pp → φ → τ + τ − channel at ATLAS, since in both cases the observed excesses over the background expectation are not in direct tension with corresponding results (in terms of search channel and integrated luminosity) from the other collaboration. The situation is more ambiguous for the excess in the A → Zh 125 search reported by ATLAS (as we will discuss in more detail in Sect. 2). In view of this situation we will focus on the parameter space that is preferred as a result of the χ 2 analysis where the A → Zh 125 search is not included. However, we investigate whether these parameter regions would also be compatible with a possible signal in the bb → A → Zh 125 channel.
As a second step of our analysis we take into account the possibility that a Higgs boson of an extended Higgs sector could also be lighter than the observed state at 125 GeV. Searches for Higgs bosons below 125 GeV have been performed at LEP [13][14][15], the Tevatron [16] and the LHC [17][18][19]. It is intriguing that also two of those searches show a 2 − 3 σ local excess around the same mass of about 96 GeV, that is not in tension with other search limits: pp → φ → γγ: CMS reported a local excess of about 3 σ in their first year Run 2 data [17], with a similar upward deviation of 2 σ local in the Run 1 data at a comparable mass [20]. The ATLAS results based on the data of the first two years of Run 2 [19] are not sensitive to the excess. e + e − → Z φ → Z bb: LEP reported a local excess of around 2 σ [21]. In previous analyses focussing just on the excesses at about 96 GeV it was shown that type II and type IV of the N2HDM can fit both excesses simultaneously [21]. This was also investigated in combination with a viable dark-matter candidate for the case where a complex singlet instead of a real singlet is considered [22]. It was furthermore demonstrated that SUSY models like the NMSSM [23][24][25] or the µνSSM [26,27] can account for the excesses at a level of roughly 1σ. On the other hand, in the MSSM neither a 400 GeV Higgs boson can be accomodated in view of the existing constraints (see the discussion in Sect. 4), nor the CMS excess at around 96 GeV can be realized [28]. Consequently, in the present paper we focus on the N2HDM and the NMSSM. We address the question, after separately analyzing possible interpretations of the observed excesses at about 400 GeV as described above (for a discussion of the tt excess in a general 2HDM, see Ref. [29]), whether the observed patterns at 400 GeV and 96 GeV can be described simultaneously within the considered models.
Our paper is organized as follows. In Sect. 2 we summarize the experimental results in the various search channels and define the various χ 2 contributions employed in our analysis. Sect. 3 is devoted to possible interpretations of the observed excesses within the N2HDM. We first investigate whether the model can fit the observed excesses at 400 GeV in the pp → A → tt and pp → φ → τ + τ − channels while being in agreement with the relevant theoretical and experimental constraints. We find that the parameter regions that would be preferred by possible signals in the two channels do not overlap with each other. On the other hand, each of the two excesses individually can be described very well by the N2HDM with Yukawa structure of type II. In a second step we then demonstrate that this model can simultaneously also describe both excesses at 96 GeV. In Sect. 4 we extend our analysis to the case of the NMSSM, where the structure of the Higgs sector is more rigid than in the N2HDM. Nevertheless, we demonstrate that also the NMSSM can fit each of the excesses at 400 GeV individually, while complying with the BSM Higgs-boson searches and the signal-rate measurements of the Higgs boson at 125 GeV. In this analysis, we show that a particularly well motivated parameter region to accommodate the tt excess is given by the alignment-without-decoupling limit of the NMSSM [30]. In this limit a light singlet-like scalar is naturally present in the Higgs spectrum. Consequently, as a next step we include the 96 GeV excesses also into the NMSSM analysis and show that only the CMS excess can be accommodated alongside with the 400 GeV excesses. In Sect. 5 these analyses are supplemented with a discussion on the prospects for future experimental tests of the parameter space regions of the N2HDM and the NMSSM that are favored by the experimental excesses. We conclude in Sect. 6. In the Appendix we provide supplementary material on different Yukawa types of the N2HDM and a discussion of the constraints from the signal at 125 GeV in the alignment-without-decoupling limit of the NMSSM.

Experimental situation
Current LHC Run 2 data of proton-proton collisions at a center-of-mass energy of 13 TeV show some tantalizing local excesses in various channels with local significances around the 3 σ confidence level (C.L.) which are summarized in the following.

Search for H, A → tt
A search for additional scalar or pseudoscalar Higgs bosons decaying to a top quark pair was performed by the CMS experiment [6]. The data set analyzed corresponds to an integrated luminosity of 35.9 fb −1 . Final states with one or two charged leptons were considered. The invariant mass of the reconstructed top quark pair system and angular variables sensitive to the spin of the particles decaying into the top quark pair were used to search for signatures of the H or A bosons, taking the interference with the SM tt background into account.
Constraints on the strength modifiers of the Htt(Att) couplings c Htt (c Att ) were derived as a function of the mass and width of the heavy Higgs boson H(A), with the mass of the bosons ranging from 400 to 750 GeV. Here, the modifiers are defined as the coupling strengths relative to the ones of a SM Higgs boson. In this study it is assumed that there is no CP violation in the Higgs sector. As a consequence the H and A Higgs bosons are CP eigenstates and do not mix with each other, which means that there are no inteference effects between H and A. The modelindependent search was performed for each CP state separately, setting the coupling modifier for the other respective CP state to zero. As a result, a moderate signal-like deviation was observed for the hypothesis of a pseudoscalar Higgs boson with the mass m A ≈ 400 GeV and a total relative width of Γ A /m A ≈ 4.5%, with a local significance of 3.5 ± 0.3 standard deviations. This translates to a global significance of 1.9 σ. The excess is also consistent with slightly different masses and widths of a possible signal, however, with less statistical significance. The deviation is driven by  the CP-odd Higgs boson, while the impact from the CP-even Higgs boson is small, because at a mass of ≈ 400 GeV the resonant production cross section of a CP-odd Higgs boson is much larger than that of a CP-even Higgs boson. This can also be inferred, for instance, from Fig. 3 of Ref. [32].
There is no corresponding search at 13 TeV in this channel yet from the ATLAS Collaboration, but there is a search for heavy pseudoscalar and scalar Higgs bosons decaying into a tt pair performed with an integrated luminosity of 20.3 fb −1 at a center-of-mass energy of 8 TeV [33]. No significant deviation from the SM prediction was observed in the tt invariant mass spectrum in final states with an electron or muon, large missing transverse momentum, and at least four jets. However, dileptonic final states and angular variables were not utilized for this analysis, and cross section limits were set only for masses ranging from 500 GeV to 750 GeV. Consequently, a combined analysis of the ATLAS and CMS results is not possible.
The χ 2 -distribution associated to the CMS analysis [6, 31], χ 2 tt is shown in Fig. 1 as a function of the coupling strength c Att for the case of a pseudoscalar Higgs boson with an assumed mass of m A = 400 GeV and different values of the width, Γ A , ranging from 0.1% to 25% of its mass. 1 Width values between 1% and 10% have a χ 2 value close or below 1. The χ 2 is given by twice the negative logarithm of the likelihood function with b i denoting the combined background yield in a given bin i of the invariant tt mass and spin 1 It should be noted that in Fig. 7 of Ref.
[6] there is a typo in the label of the vertical axis, which should read − ln (L(g Att )/L SM ).
analyzing distribution, s A R,i and s A I,i the signal yields in a given bin for the resonant and interference part, respectively, ν the vector of nuisance parameters on which the signal and background yields generally depend, and n i the observed yield in data. The external constraints on the nuisance parameters are taken into account in the likelihood via a product of corresponding probability density functions, G( ν). The χ 2 tt distribution is normalized by L max = L(c Att = 0.94; m A = 400 GeV, Γ A = 4.5% m A , ν max ) such that it vanishes for the most likely choice of mass, width, coupling strength and nuisance parameters, ν max , providing the best description of the data.

Search for φ → τ + τ −
A search for heavy neutral Higgs bosons decaying into two tau leptons was performed by the ATLAS Collaboration utilizing data corresponding to an integrated luminosity of 139 fb −1 [7]. The search covered the mass range 0.2-2.5 TeV of the heavy resonance, decaying into τ + τ − with at least one τ lepton decaying into final states with hadrons. The natural width of the scalar boson is assumed to be negligible compared to the experimental resolution. Upper limits on the production cross section times branching fraction for a Higgs boson φ produced via gluon-gluon fusion (ggF) and b-quark associated production (bbφ) were derived. The limits were calculated from a statistical combination of three different final states, involving one τ lepton decaying into a neutrino and hadrons and one τ lepton decaying into neutrinos and an electron (eτ h ) or into neutrinos and a muon (µτ h ), and involving two τ leptons decaying into a neutrino and hadrons each (τ h τ h ). For ggF a local excess was observed in the data at the 2.2 σ C.L. at m φ = 400 GeV, while for bbφ production a local excess at the 2.7 σ C.L. at m φ = 400 GeV was found. We have performed a combined statistical analysis of the excesses in both production modes assuming m φ = 400 GeV in terms of a two-dimensional likelihood function as provided by ATLAS [7,34]. The point with the highest excess is located at σ(gg → φ) × BR(φ → τ + τ − ) = 20.19 fb and σ(bb → φ) × BR(φ → τ + τ − ) = 38.37 fb.
In our N2HDM and NMSSM analyses all points receive a contribution χ 2 τ + τ − relative to the best-fit point according to the two-dimensional likelihood function. In both models it is in principle possible that more than one Higgs boson can contribute to the observed excess. Thus, we define the signal as the incoherent sum of all Higgs bosons with masses in the range (400 ± 40) GeV, and use the χ 2 τ + τ − values of the sum according to the experimentally measured likelihood function for m φ = 400 GeV. It should be noted, however, that in the N2HDM only the CP-odd Higgs boson with a mass of m A ≈ 400 GeV can give rise to a sizable contribution to the excess. This is due to the fact that for the case where a second (CP-even) Higgs boson h 2 has a similar mass (i.e. 360 GeV ≤ m h 2 ≤ 440 GeV) the applied constraints require that h 2 has to be an almost pure gauge singlet in the parameter space region that is relevant for the τ + τ − excess, and its contributions to the signal cross section are therefore suppressed compared to the ones of the CP-odd state.
The corresponding search of the CMS Collaboration for additional neutral Higgs bosons in the τ + τ − final state was done on a dataset corresponding to an integrated luminosity of 35.9 fb −1 [18]. Model-independent limits were set on the product of the production cross section times the branching fraction for the decay into τ leptons for both the ggF and the bbφ production modes. The search for heavy resonances was performed over the mass range 0.09-3.2 TeV, assuming a narrow width. Final states involving one τ lepton decaying into neutrinos and an electron and one τ lepton decaying into neutrinos and a muon (eµ) were considered as well as eτ h , µτ h and τ h τ h final states. No excess over the SM background was found in this analysis. However, due to the lower amount of data included, the expected and observed cross section limits of the CMS analysis are weaker than the ones from ATLAS utilizing the full Run 2 dataset. Thus, the signal intepretation of the ATLAS excess at around 400 GeV at the best-fit point is compatible with the expected (observed) upper limits from the CMS analysis at the level of 1.1(1.9)σ.

Search for A → Zh 125
The ATLAS Collaboration performed a search for new resonances decaying into a Zh 125 in ννbb and + − bb final states, where ± = e ± or µ ± [8]. The data used corresponds to a total integrated luminosity of 36.1 fb −1 . The search was conducted by examining the reconstructed invariant mass distribution of Zh candidates for evidence of a localised excess in the mass range of 220 GeV up to 5 TeV. The results of the search for the CP-odd scalar boson A in two-Higgs-doublet models were interpreted in terms of constraints on the ggF and b-quark associated (bbA) production cross-sections times the branching fraction of A → Zh 125 and the branching fraction of h 125 → bb.
In the search for bbA production, an excess of events was observed for a resonance mass around 440 GeV, mainly driven by dimuon final states with three or more b-tags required. The local significance of this excess with respect to the background-only hypothesis was estimated to be 3.6 σ, and the global significance, accounting for the look-elsewhere effect, was estimated to be 2.4 σ. In the search for ggF production, agreement between the data and the background-only hypothesis was found in a preliminary update analysing the full Run 2 data set corresponding to an integrated luminosity of 139 fb −1 [10].
The respective analysis performed by the CMS Collaboration for a heavy pseudoscalar boson A decaying to Zh 125 similarly considered final states where the Higgs boson decays to a bottom quark and antiquark, and the Z boson decays either into a pair of electrons, muons, or neutrinos [9]. The analysis was performed using a data sample corresponding to an integrated luminosity of 35.9 fb −1 . Exclusion limits were set in the context of 2-Higgs-doublet models in the A boson mass (m A ) range between 225 and 1000 GeV on the ggF and bbA production cross-sections times the branching fraction of A → Zh 125 and the branching fraction of h 125 → bb. In this search, however, the data was found to be consistent with the background expectations over the whole range of the reconstructed resonance mass m A .
In view of this somewhat ambiguous experimental situation we do not incorporate the excess in the A → Zh 125 channel into our χ 2 analysis. We rather investigate whether the parameter regions in the N2HDM and the NMSSM that are favoured as a result of the other observed excesses discussed in this section could give rise to a possible signal also in the A → Zh 125 channel. For a discussion of the ATLAS excess in the A → Zh 125 channel in the context of the 2HDM see Ref. [35], and in the context of a general NMSSM see Ref. [36].

Searches for Higgs bosons in the low mass region
As discussed above, in the context of searches for additional Higgs bosons it is important to take into account also the mass region below the observed state at 125 GeV. Higgs bosons with relatively low masses can be searched for at the LHC via diphoton resonant searches. CMS found a deviation between the expected and observed exclusion limits at a mass of about 96 GeV in the Run 1 data set, with a local significance of about 2 σ [20]. Remarkably, the updated analysis including the first Run 2 data found an excess of about 3 σ at a comparable mass [17]. Taking into account the accumulated data at 7, 8 and 13 TeV center-of-mass energy, CMS reported a signal interpretation of the excesses, corresponding to a signal strength of Here the signal strength µ CMS is defined, as usual, as the resonant signal cross section assuming a scalar particle around 96 GeV, normalized to the expected signal cross section assuming a SM Higgs boson at the same mass. Similar searches in the diphoton final state have also been performed by ATLAS using 80 fb −1 of the 13 TeV dataset [19]. Here, only a very small excess of events over the SM background has been found around 96 GeV. However, overall the resulting upper limits on the cross sections found by ATLAS are substantially above the corresponding CMS exclusion limits. Consequently, the ATLAS result cannot exclude the signal interpretation of the CMS excess.
The CMS excesses have gained considerable interest in the literature due to the fact that it is consistent with another excess found already at the Large Electron Positron collider (LEP). At LEP, Higgs bosons in this mass range could be searched for also in hadronic final states due to the much lower backgrounds. Here an excess was reported in the search e + e − → ZH → Zbb at a mass of about 98 GeV, with a local significance of 2.3 σ [14]. The signal strength for the excess was found to be [23] µ exp LEP = 0.117 ± 0.057 .
The mass resolution of the LEP excess is rather limited due to the hadronic final state, such that both the LEP and the CMS excesses could have a common origin. Various different scenarios to accommodate both excesses at about 96 GeV have already been discussed in the literature. 2 In particular, possible interpretations comprise type II and type IV of the N2HDM [21] and SUSY models like the NMSSM [23][24][25] and the µνSSM [26,27], where the SUSY models can account for the excesses at a level of roughly 1σ. In the present paper we will first focus on the observed excesses at ≈ 400 GeV and then investigate whether a simultaneous realization of also the excesses at 96 GeV is possible. In order to quantify how well the excesses are fitted for each parameter point in our numerical analysis, we define the χ 2 regarding the excesses at 96 GeV via Here µ CMS and µ LEP are the theoretical predictions. For the SM, for which we have set µ CMS = µ LEP = 0, we find a penalty of χ 2 SM,96 = 13.2 from these contributions.

N2HDM interpretation
In this section we will discuss possible realizations of the described excesses in the Higgs searches within the context of the N2HDM. We start by introducing the model and its parameters, followed by a discussion of the relevant theoretical and experimental constraints. The numerical analysis will be divided into two parts. First, we discuss how the excesses at 400 GeV can be accommodated, 2 A recent overview about the various BSM models that can explain both the CMS and the LEP excesses can be found in Ref. [37]. and whether a simultaeneous realization of both the tt and the τ + τ − excesses is possible. In a second step, we additionally take into account the excesses at 96 GeV in order to address the question whether these can be realized in the parameter region in which one of the excesses at 400 GeV is accommodated.

Model definitions
The N2HDM contains two doublet scalar fields Φ 1,2 and a real scalar singlet field Φ S [38]. In the physical basis, the Higgs particle sector consists of a total of three CP-even Higgs bosons h i , a CP-odd state A and two charged Higgs bosons H ± . In order to avoid tree-level flavorchanging neutral currents (FCNC) a Z 2 symmetry is assumed as in the usual 2HDM, which is only softly broken by the terms −m 2 12 Φ † 1 Φ 2 + Φ † 2 Φ 1 appearing in the Lagrangian. In addition, the scalar potential respects a second Z 2 symmetry, under which only the singlet field Φ S is charged. This symmetry is however broken spontaneously when the singlet field has a non zero vacuum expectation value (vev), Φ S = v S , which we will assume throughout the analysis. The breaking of the electroweak gauge symmetry originates from the vevs of the doublet fields Φ 1, A crucial parameter for the analysis of the experimental excesses is the ratio of the doublet vevs given by The three real neutral components of Φ 1,2,S mix to form the mass eigenstates h 1,2,3 with masses m h 1 < m h 2 < m h 3 . The mixing in the CP-even scalar sector can be described by the three mixing angles α 1,2,3 , defining the matrix At least in principle, the properties of each of the Higgs bosons h i could be such that it can be identified with the SM-like Higgs boson at 125 GeV. The masses of the CP-odd Higgs boson and the charged Higgs bosons are denoted by m A and m H ± , respectively. Assigning consistent charges under the Z 2 symmetry to the fermions yields the typical four types of Yukawa structures that are familiar from the 2HDM. Depending on the type, the couplings of the scalar particles to the fermions will be different. The relevant parameters defining the strength of the couplings normalized to the one of a hypothetical SM Higgs boson of the same mass can be given in terms of tan β and the mixing angles α i of the CP-even sector.
As already mentioned above, in a first step we will analyze whether the N2HDM can realize the excesses at 400 GeV. We will go a step further in the subsequent analysis and investigate whether and how the excesses at 400 GeV and at 96 GeV can be realized simultaneously. Our approach for covering the relevant parameter space will be discussed in detail in App. A.1. Before, we briefly summarize the relevant theoretical and experimental constraints that we take into account in our analysis.

Theoretical and experimental constraints
Numerous theoretical constraints on the N2HDM parameter space have to be taken into account in order to exclude unphysical parameter configurations. To assure the presence of a viable electroweak minimum, described by the vevs v 1 , v 2 and v S as explained in Sect. 3.1, we demand that the scalar potential is bounded from below. The necessary conditions are given in terms of the quartic scalar couplings which define the behavior of the potential for large field values [39]. Furthermore, we verified that the electroweak minimum is either the global minimum of the scalar potential, or, in case other deeper minima exist, that the lifetime of the electroweak minimum is large compared to the age of the observable universe. Furthermore, in order to verify that the perturbative treatment of the model is adequate, we checked that the perturbative unitarity constraints are fulfilled. These are also given in terms of the quartic scalar couplings and set upper limits on their absolute values [40]. All the theoretical constraints mentioned here have been applied using the public code ScannerS [40][41][42][43][44] in combination with the public code EVADE [45,46] for the calculation of the lifetime of metastable electroweak minima.
Besides the mentioned theoretical constraints, we have also included the most relevant experimental constraints on the N2HDM parameter space into our analysis: (i) We perform a χ 2 test regarding the measurements of the signal rates of h 125 using the public code HiggsSignals v.2.6.0 [47][48][49][50]. In the following, we will denote the result of the HiggsSignals test by χ 2 125 , which is constructed taking into account n obs = 107 observables. The value of χ 2 125 for each N2HDM parameter point will be compared with the SM prediction χ 2 SM,125 = 84.42. (ii) We verify for each parameter point that none of the N2HDM scalars is excluded by searches for additional Higgs bosons at the LHC, the Tevatron and at LEP by making use of the public code HiggsBounds v.5.9.0 [51][52][53][54][55][56]. For the considered parameter point this code selects the most sensitive search channel for each Higgs boson based on the expected experimental sensitivity. It then compares the model predictions for the signal rates µ theo. with the observed upper limits at the 95% C.L., µ excl. , and rejects a point if for one of the Higgs bosons the ratio r ≡ µ theo. /µ excl. is larger than one. (iii) We include constraints from electroweak precision observables (EWPO) in terms of the socalled oblique parameters S, T and U [57,58]. Using the implementation of the code ScannerS we exclude points for which the χ 2 value obtained by comparing to the fit result of Ref. [59] deviates by more than 2 σ. The N2HDM model predictions for S, T and U are calculated at the one-loop level following Refs. [60,61]. (iv) Mainly due to the presence of the charged Higgs boson in the N2HDM, constraints from flavor physics are important in some regions of parameter space. Since the extra field of the N2HDM compared to the 2HDM is a gauge singlet, it is sufficient to take over the bounds for most of the flavor physcis observables from the 2HDM [21,44]. We apply the flavor constraints from Ref. [59] as implemented in ScannerS. They yield a lower limit m H ± 550 GeV in both type II and type IV. In type I and type III no such lower limit on m H ± can be established. In all types the flavor constraints give rise to a lower limit on tan β, where the specific value is different in each type, and the strongest constraint is obtained in the type IV (N  The states h a,b,c are automatically arranged according to the mass-ordered notation h 1,2,3 for each point by ScannerS.
Except the χ 2 result of HiggsSignals regarding the signal rates of h 125 , the experimental constraints described in this section have been taken into account in terms of a hard cut, i.e., a parameter point is regarded to be excluded if the considered experimental constraints are not fulfilled for this point. As will be described in Sect. 3.4, we combine χ 2 125 with the χ 2 values for the different observed excesses. In this way we quantify how well the considered models describe both the experimental results for h 125 and the potential signals at 400 GeV, as well as in a second step also the potential signals at 96 GeV.

A Higgs boson at 400 GeV in type II
In order to accommodate the tt excess one needs values of |c Att | 0.5 (see Fig. 1). Since in all four Yukawa types of the N2HDM one finds |c Att | = 1/ tan β, it is easy to understand that beyond the fact that different constraints apply in each type, all four Yukawa types are very similar regarding the tt excess. The situation is more complicated for the τ + τ − excess, where it is required that two conditions are fulfilled: Firstly, in order to accommodate values for the bb associated production cross section of A of the same size or larger than the gluon-fusion production cross section, the condition |c Abb | |c Att | has to be fulfilled. Secondly, in order to obtain sufficiently large values for BR(A → τ + τ − ), the condition |c Aτ + τ − | |c Att | has to be fulfilled. The coupling coefficients c Abb and c Aτ + τ − are different in each type, either proportional to tan β or its inverse. One finds that only the type II Yukawa structure is able to satisfy both conditions (see App. A.1 for a detailed explanation), such that it is the only type that can potentially accommodate the τ + τ − excess. In the following, we will focus our numerical discussion on the N2HDM type II, while in the appendix we will also provide a discussion of the N2HDM type IV.
In this section we start by investigating in which parameter regions of the type II N2HDM the excesses at 400 GeV can be accommodated. 4 We give in Tab. 1 the ranges of the free parameters used here. The choice of the input parameters is the default one of the public code ScannerS [41,44] that was used to generate the parameter points respecting the various theoretical and experimental constraints discussed in Sect. 3.2. Each set of input parameters corresponds to a unique benchmark point under the assumption c h b tt · c h b V V > 0. This condition arises from the constraints on the signal rates of the state h 125 (here h 125 = h b , as defined in Tab. 1). In the following, the CP even scalars will be denoted as h a,b,c when referring to the input parameters as shown in Tab. 1, or alternatively as h 1,2,3 , where the mass ordering m h 1 < m h 2 < m h 3 is implied.
Motivated by the pattern of the observed excesses described in Sect. 2, our scan in the N2HDM is performed such that the observed excesses at 400 GeV are mainly associated with the CPodd Higgs boson A (a CP-even Higgs boson in this mass region could yield only a very small contribution to the signal cross section as a consequence of the constraints and the required large singlet component), and we therefore fix m A = 400 GeV in the scan. The scan range of tan β = [0.5, 12.5] is chosen in view of our qualitative discussion above of the observed excesses and of the bounds from searches for additional Higgs bosons. Moreover, the mass m h b is fixed to m h b = 125.09 GeV[3], and the corresponding couplings are restricted to ranges that allow for the presence of a Higgs boson h b that resembles the properties of the Higgs boson that has been detected at 125 GeV. The lower limit of m H ± ≥ 550 GeV is chosen in view of the constraints from flavor physics observables (see Sect. 3.2). The remaining mass parameters were varied up to values of 1 TeV, except for 10 ≤ v S ≤ 1.5 TeV. The mixing matrix element R a3 , see Eq. (3.1), is scanned over all theoretically possible values, and both possibilities for sign(R b3 ) are taken into account. Due to the large number of free parameters, it is not possible to cover the entire parameter space region defined in Tab. 1 sufficiently well without applying other theoretical prejudices on the model parameters. However, since the properties of the pseudoscalar A are to a large extent determined by the parameters m A and tan β, robust conclusions can be drawn with respect to the realization of the excesses by randomly sampling the parameter space, without further theoretical restrictions, for instance, in terms of priors.
We quantify the results of this scan in terms of the total χ 2 arising from the comparison of the N2HDM predictions with the observed results for the excesses at 400 GeV in the tt and τ + τ − channels as well as with the signal-rate measurements of the Higgs boson at about 125 GeV, such that The evaluation of the individual χ 2 contributions in Eq. (7) has been described in Sect. 2 and Sect. 3.2. In our analysis we consider all points with χ 2 ≤ χ 2 SM as acceptable, where χ 2 SM is evaluated via Eq. (7) assuming no signal contributions to the excesses. For χ 2 tt and χ 2 τ + τ − we remind the reader that the values are defined relative to the best-fit points. One finds χ 2 SM,tt = 13.92 and χ 2 SM,τ + τ − = 9.99 assuming the SM hypothesis. All remaining theoretical and experimental constraints are included as hard cuts, either allowing or excluding a parameter point (see Sect. 3.2). As a consequence of the condition χ 2 ≤ χ 2 SM , the sample of parameter points displayed in our plots below will also contain points with χ 2 125 ≈ χ 2 SM,125 that potentially do not describe either of the excesses with a confidence level below 2 σ. However, since our results allow a point-by-point , respectively. In both plots the best-fit point, labelled as min(χ 2 ), is indicated with a magenta star. The point with the smallest value of χ 2 τ + τ − + χ 2 125 , labelled as min(χ 2 τ + τ − + χ 2 125 ), is indicated with a cyan star.
comparison 5 of the individual χ 2 values, it should be clearly visible which of the displayed the scan points yield a good description of one or both of the excesses, i.e. featuring small values of χ 2 tt and/or χ 2 τ + τ − , while mantaining an acceptable χ 2 125 . As explained in App. A.1, the properties of the A-boson are driven by the dependence of its couplings on tan β. In the left plot of Fig. 2 c Att (= 1/ tan β) is displayed as a function of tan β for different values of Γ A /m A . One can see that the best-fit values for c Att showing the best agreement with the tt excess, which are indicated for different width hypotheses by the horizontal dashed lines, are reached only for small values of tan β 2. As discussed in Sect. 2 (see Fig. 1) the experimental excess is most pronounced for Γ A /m A ≈ 4%. This feature manifests itself in the plot by the fact that the largest displayed values of Γ A /m A correspond to the largest best-fit values of c Att . Overall, the N2HDM of type II yields a good description of the tt excess in the region of small values of tan β, 1.2 tan β 2.5. The lower bound on tan β arises here from the cross section limits on the process pp → tt(A) → tt(tt), as reported by CMS making use of the full Run 2 data set [62]. For values of tan β > ∼ 2.5 we find values of c Att substantially below the experimental best-fit value even for the lowest value of Γ A /m A = 0.5% taken into account in the analysis. Thus, we conclude that the tt excess can be well described in the N2HDM of type II with low tan β, while for tan β > ∼ 2.5 the compatibility with the tt excess is reduced. To further  Figure 3: χ 2 tt (blue) and χ 2 τ + τ − (orange) and χ 2 tt + χ 2 τ + τ − (green) in dependence of tan β. The horizontal blue, orange and green dashed lines indicate the value of χ 2 SM,tt , χ 2 SM,τ + τ − and χ 2 SM,tt + χ 2 SM,τ + τ − , respectively.
illustrate the fit result we show in the middle plot of Fig. 2 the subset of points for which we find 2.0% ≤ Γ A /m A ≤ 3.0% in the plane of m A and c Att , where we also indicate the expected and observed exclusion limits from the experimental analysis assuming a width of Γ A /m A = 2.5% as published in Ref. [6]. One can see that the points stretch over the values of c Att ≈ 0.8 that provide the best description of the tt excess.
In the right plot of Fig. 2 we show the parameter points in the plane of the signal rates regarding the τ + τ − excess [7], where the production cross sections were calculated with the public code SusHi [64,65]. Here, the colors of the points indicate the value of tan β. Our scan results in two distinct parameter regions in the displayed plane. The origin of these regions will be discussed below. We also show the 1σ and 2σ ellipses of χ 2 τ + τ − with the dark blue and the yellow lines. One can see that, as expected following the discussion in App. A.1, there are only points with larger values of tan β > 6 inside the 1σ ellipse. Comparing with the region of tan β in which the tt excess can be well described, we conclude that the N2HDM cannot provide a good simultaneous fit of both excesses.
This feature can be seen more clearly in Fig. 3, where the χ 2 distributions of the points regarding both excesses are displayed, together with their sum. The values of χ 2 tt are shown in blue and the ones of χ 2 τ + τ − in orange for all the parameter points of this scan. In agreement with the observations described above, low values of χ 2 tt 4 can only be found for tan β 2.5. Low values of χ 2 τ + τ − 4 are found for tan β 5.5. Consequently, the sum of both χ 2 values, as indicated by the green points, has the smallest values in both tan β ranges suitable for the tt or the τ + τ − excess (but does not drop below a value of about 10), while larger values are found in between both ranges, at 3.5 tan β 5.5. Therefore, it is not possible to simultaneously accommodate both excesses in the N2HDM.
The two distinct regions found in the right plot of Fig. 2 correspond to the two separate distributions of points that are visible for χ 2 τ + τ − in the region of large tan β. The left one corresponds to the right region of points in Fig. 2, which have larger values of σ(gg → A → τ + τ − ) and do not cross the center of the ellipses. Consequently, χ 2 τ + τ − does not exactly reach zero for this set of , respectively. In both plots the best-fit point, labelled as min(χ 2 ), is indicated with a magenta star. The point with the smallest value of χ 2 points. On the contrary, the left band of points in Fig. 2 exactly crosses the center of the ellipses for values of tan β ≈ 9.5. Therefore, the distribution of points displaying χ 2 τ + τ − near this value of tan β has its minimum at zero. The two different regions have their origin in the presence (or absence) of the decay of A to a lighter Higgs boson h i and a Z boson. If such a decay is possible, it leads to a reduction of BR(A → τ + τ − ).
The impact of such an additional decay mode of the CP-odd Higgs boson A is illustrated in Fig. 4, in which we show the same as in the right plot of Fig  of points that do not have m ha < 300 GeV are in the wrong sign Yukawa coupling regime. For these points the presence of the decay of A into h 125 and a Z boson improves the fit to the τ + τ − excess. However, such a signature has also been probed in direct searches. For the gg production mode (see below for a discussion of the observed excess in the bbA production mode [8]) the limits obtained from the search for gg → A → Z(h 125 ) → Z(bb) recently reported by ATLAS [10] exclude parameter points in the wrong sign Yukawa coupling regime with values of tan β 8.5. This leads to the lower cut on the left branch of points in Fig. 2 (right) and Fig. 4, and also to the lower bound on tan β in the right branch of points in Fig. 3.
It should also be noted that the parameter point with the smallest value of χ 2 τ + τ − + χ 2

125
(indicated by the cyan star) in Fig. 4 lies on the right branch of points. Since this branch is further away from the center of the ellipses, this branch has overall slightly larger values of χ 2 τ + τ − . The fact that we nevertheless find the cyan star in this branch indicates that for these points the compatibility with the measured properties of h 125 can be better (yielding a lower χ 2 125 ) compared to the left branch of points. Thus, parameter points in which the additional decay channels A → Zh 1,2 are relevant (left branch) are associated with mildly larger deviations of the signal rates of h 125 compared to the SM prediction. However, we checked that the values of χ 2 125 can also be very close to the SM value χ 2 SM,125 ≈ 84 in the left branch. As a consequence, the accuracy of the signal rate measurements of h 125 at the (HL)-LHC [66] will not be sufficient to fully discriminate between the points in the two branches.
As already mentioned above, in the N2HDM type II the presence of one of the decays A → h 1,2 Z with sizable branching ratio can have the effect that the τ + τ − excess is described so well that one encounters a vanishing contribution to χ 2 τ + τ − . Regarding the points in the wrong sign Yukawa coupling regime, it is interesting to note that a local 3.6 σ excess at roughly 400 GeV was found by the ATLAS Collaboration in the search bb → A → h 125 Z [8]. The corresponding search for the gg production mode has meanwhile been updated including the full Run 2 dataset, and no excess was found at m A ≈ 400 GeV. As explained above (see also the discussion in Sect. 2.3) the obtained limit gives rise to the lower bound on tan β for the left branch of points in Fig. 2 (right) and Fig. 4 , the blue points in the right plot are in the wrong sign Yukawa coupling regime. One can see that the points in the wrong sign Yukawa coupling regime could indeed contribute to the excess observed by ATLAS in the bb → A production mode, while at the same time being compatible with the limits obtained for the gg production mode. On the other hand, for the points that are not in the wrong sign Yukawa coupling regime in our scan we do not find a sizable contribution to the σ(bb if kinematically open, are suppressed for the case where the Higgs boson in the final state is h 125 .
For the points that are not in the wrong sign Yukawa coupling regime our analysis has revealed that the presence of a second Higgs boson h a with a mass below 300 GeV can improve the description of the τ + τ − excess. It is tempting in this context to entertain the possibility that this additional Higgs boson could have a mass of about 96 GeV and be the origin of the observed excesses in this mass range, see the discussion in Sect. 2.4. Within the context of the N2HDM of type II, the simultaneous realization of the τ + τ − excess and the observed pattern around 96 GeV leads to a very interesting scenario that is very predictive and can serve as a benchmark scenario that can be probed in future experiments. We will investigate this possibility in Sect. 3.5. Before turning to this discussion we will investigate whether the excesses at 96 GeV can also be accommodated in combination with the tt excess. The latter, as discussed in detail in App. A.1, can be realized both in the type II and the type IV N2HDM. A corresponding parameter scan for the type II is described in Sect. 3.4, and for completeness the according scan in type IV can be found in App. A.2.

Higgs bosons at 96 GeV and 400 GeV for low tan β in type II
In this section we present the analysis for the scan in the low tan β regime of the type II N2HDM, which is dedicated to accommodate the tt excess at 400 GeV in combination with the excesses at 96 GeV. The ranges of the free input parameters are given in Tab. 2, where for the scan in the present section the tan β range specified as tan β low is adopted.
We start by briefly commenting on our choice of parameter ranges. 7 The lightest Higgs boson h 1 plays the role of the candidate that is associated with the excesses at 96 GeV. Thus, the mass m h 1 was chosen to be in this range. The SM-like Higgs boson is h 2 with a mass of m h 2 = 125.09 GeV. As before, the tt excess will be described in terms of the CP-odd Higgs boson with a mass of m A = 400 GeV. The scan range of the charged Higgs-boson mass m H ± = [550, 1000] GeV is chosen due to a lower limit of m H ± 550 GeV arising from the flavor constraints. The indirect constraints from the EWPO lead to the restriction that the third CP-even Higgs boson should have a similar mass m h 3 ≈ m H ± , since all remaining Higgs bosons have fixed masses much below the lower limit on m H ± . Thus, m h 3 is scanned in the same mass range as m H ± . The range of tan β = tan β low gives rise to sizable couplings of the CP-odd Higgs boson A to top quarks, which is desired for the description of the tt excess (see the discussion in App. A.1). The ranges for the squared couplings of the SM-like Higgs boson h 2 were taken over from Ref. [21], where it was shown that such values are suitable for explaining the excesses at 96 GeV in the N2HDM. To be precise, these ranges allow a substantial mixing between h 1 and h 2 , such that h 1 can have sizable couplings to quarks and gauge bosons, while the properties of h 2 are still in agreement with the signal-rate measurements of the discovered Higgs boson at around 125 GeV. The remaining parameters in Tab. 2, R 13 , m 12 and v S , were scanned over a wide range, and both possible signs of R 23 were taken into account. We emphasize that the requirement of simultaneously explaining the tt excess and the excesses at 96 GeV practically fixes the entire scalar spectrum of the model. The only remaining free mass parameter is the common mass scale of h 3 and H ± . However, this scale cannot be much larger than m A (which is fixed to 400 GeV) as a consequence of the applied constraints (see also the discussion below).
The results of this scan will be quantified in terms of the total χ 2 given by the contributions arising from the fit to the excesses at 96 GeV, the fit to the tt and the τ + τ − excesses at 400 GeV, and the signal-rates measurements of the Higgs boson at 125 GeV, such that We include here for consistency also the contribution from the τ + τ − excess even though no sizable contribution to the correspoding signal cross sections are expected in the low tan β regime. All points with χ 2 ≤ χ 2 SM are considered to be acceptable, where χ 2 SM is evaluated assuming no signal contributions to the excesses. Their contribution to χ 2 SM amounts to χ 2 SM,96 + χ 2 SM,τ + τ − + χ 2 SM,tt = 37.2. The best-fit point in our scan is determined by the lowest χ 2 value according to Eq. (8).
In Fig. 6 we show the values of the coupling coefficient c Att in dependence of tan β for all the parameter points with χ 2 ≤ χ 2 SM . One can see that the parameter points lie in a narrow range of 1.1 tan β 1.6. The preference for those low values of tan β is related to the fact that for values of tan β 2 no sizable contribution to the tt excess is present. The feature that no points are displayed in our scan for higher values of tan β although our scan has been performed up to tan β = 4 is mainly due to the sampling used in our scan in combination with the fact that many points in this region do not fulfill the condition χ 2 ≤ χ 2 SM . For the allowed points, the size of the couplings c Att and the values for Γ A /m A , indicated by the colors of the points in the left plot of Fig. 6, provide a very good fit to the tt excess. In particular, for values slightly above tan β = 1 we find c Att ≈ 0.8 and Γ A /m A ≈ 2.5%, which gives rise to values of χ 2 tt < 1. In the middle plot of Fig. 6 the colors of the points indicate the value of the charged Higgs boson mass m H ± . Here it is important to note that no values above m H ± ≈ 750 GeV were found. This upper limit arises from theoretical constraints on the allowed values of the quartic couplings under the Best fit cAtt for ΓA/mA = 3.0% Best fit cAtt for ΓA/mA = 2.5% Best fit cAtt for ΓA/mA = 2.0% Best fit cAtt for ΓA/mA = 1.5% Best fit cAtt for ΓA/mA = 1.0% Best fit cAtt for ΓA/mA = 0.5%  The width of the CP-odd Higgs boson A normalized to its mass ranges between roughly 1.5% and 3.5%, as can be seen in the left plot of Fig. 6. These values are substantially larger than the ones that one would obtain within the 2HDM for a CP-odd Higgs boson at 400 GeV. This is due to the fact that in the considered N2HDM scenario A has additional decay channels available into h 1,2 and a Z boson. In the experimental analysis [6] it was found that for values of Γ A /m A ≈ 2.5% the expected 95% C.L. exclusion limit is at roughly c Att ≈ 0.6, while the observed exclusion limit varies (depending on the precise value of Γ A /m A ) between c Att ≈ 0.9 and c Att ≈ 1.1. Thus, our scan essentially covers the range of c Att that is associated with the observed excess, as can be seen in Fig. 6. This is also reflected in the very small values of χ 2 tt , as discussed in the following. We show the parameter points of the scan in the µ CMS -µ LEP plane in Fig. 7. The color coding of the points indicates the values of χ 2 tt and ∆χ 2 125 ≡ χ 2 125 − χ 2 SM,125 in the left and right plot, respectively. One can see that a large fraction of points lies inside the 1σ region of χ 2 96 while at the same time reproducing the excess observed for c Att , since many points have χ 2 tt ≈ 0. Thus, we conlcude that the low tan β region of the type II N2HDM is well suited for accommodating the excesses at 96 GeV in combination with the tt excess at 400 GeV. Moreover, we can observe in the right plot of Fig. 7 that several points have values of ∆χ 2 125 very close to zero. Consequently, the fit to the excesses can be realized without being in tension with the signal rates measurements of the SM Higgs boson. The best fit point with a value of χ 2 = 97.97 is substantially below the SM value χ 2 SM = 121.61. We also performed a scan in the low tan β region of type IV with the goal of explaining the tt excess and the excesses at 96 GeV. We found that also this type can account for the excesses, but with a slightly worse fit result regarding the CMS excess and the properties of h 125 . This scan is presented in App. A.2 for completeness and to allow for a comparison to the type II results.

Higgs bosons at 96 GeV and 400 GeV for large tan β in type II
As was discussed before, the type II N2HDM is the only candidate of the four different Yukawa types that can potentially accommodate the τ + τ − excess at 400 GeV in combination with the two excesses at 96 GeV. The existence of a region of parameter space fulfilling the relevant criteria (see the discussion in App. A.1) depends on the simultaneous enhancement of σ(bb → A) and BR(A → τ + τ − ), such that both the signal rates in the gg production mode and the bb production mode have the adequate size for providing a good description of the data. Here we present a scan in the high tan β regime dedicated to answer the question whether there is a tan β range for which the observed excesses can be accommodated simultaneously. The input parameters are the same as the ones in the low tan β regime (as shown in Tab. 2), but with considerably larger values for tan β in the range 6 ≤ tan β ≤ 12, and the upper limit of c h 2 V V ≤ 0.9 is removed. 8 As before, we refer to the point with the lowest χ 2 value as the best-fit point, where the same definition as in Eq. (8) was used to obtain the total χ 2 value. Moreover, in the results of our scan we only display points with a lower value of the overall χ 2 than the one of the SM, i.e., Type II: Type II: Note that for consistency we again include both the contributions χ 2 τ + τ − and χ 2 tt in the fit. However, for the high tan β values considered here no sizable contribution to the tt excess is expected. The condition χ 2 ≤ χ 2 SM is practically unchanged depending on whether χ 2 tt is taken into account here or not, as can also be seen in Fig. 3, in which a difference of χ 2 SM,tt − χ 2 tt 2 is visible for tan β ≥ 6.
We show the parameter points of the type II scan at high tan β in the plane of the signal cross sections regarding the τ + τ − excess in Fig. 8 and regarding the signal rates of the excesses at 96 GeV in Fig. 9. In the left plot of Fig. 8 one can see that, as anticipated from the couplings of the A boson to fermions in type II, the signal rate for the bbA production mode grows with increasing values of tan β, while σ(gg → A → τ + τ − ) shows less sensitivity to tan β. For values of 6 ≤ tan β ≤ 11 the points lie within the 1σ region regarding χ 2 τ + τ − , which is indicated by the dark blue ellipse. Moreover, for tan β ≈ 10 the points approximately lie at the center of the ellipses, indicating that the observed excesses are well described. We emphasize that this is a non-trivial result. Given the fixed mass spectrum in this scenario, and the fact that the presence of h 96 gives rise to the additional decay channel A → Zh 96 , the location of the band on which the points are found in Fig. 8 is fixed approximately (the band corresponds to the left branch of points in Fig. 4). Hence, a type II N2HDM CP-odd Higgs boson A with a mass of m A = 400 GeV, in combination with a CP-even Higgs boson h 1 at around 96 GeV, yields a predicted pattern that precisely matches the one of the excesses observed by ATLAS in the Higgs searches in the τ + τ − final state.
In the right plot Fig. 8 we show the same parameter points, however with the color of the points indicating the value of χ 2 96 . One can see that a large fraction of points that lie within the 1σ ellipse of χ 2 τ + τ − also have very small values of χ 2 96 . Consequently, the parameter space covered in our scan contains parameter points in which the excesses at 96 GeV are accommodated simultaneously with the τ + τ − excess at 400 GeV. One can see a slight correlation between the cross section σ(gg → A → τ + τ − ) and χ 2 96 . Larger values of the cross section correspond to, on average, lower values of χ 2 96 . In Fig. 9 we show the results of our scan in the plane µ CMS -µ LEP . In the left plot one can see that ∆χ 2 125 , indicated by the colors of the points, ranges from values of ≈ 3.5 to 23. In the scan in the low tan β regime of type II also smaller values closer to zero could be found (see right plot of Fig. 7). Thus, the scenario in the high tan β regime might be associated with larger deviations of the properties of the SM Higgs boson at 125 GeV compared to the SM prediction. However, for the points with the smallest values of ∆χ 2 125 the deviations of the signal rates of h 125 compared to the SM prediction are much below the current experimental uncertainties, and also more precise future measurements at the HL-LHC might not be sufficient to entirely probe this scenario. One can also see in the right plot of Fig. 9 that many points with low values of χ 2 τ + τ − lie within the 1σ ellipse regarding χ 2 96 . The best-fit point, indicated by the magenta star, is close to the central point of the ellipse. We conclude that the type II N2HDM is perfectly capable of accommodating the excesses at 96 GeV in combination with the τ + τ − excess for values of tan β ≈ 10.
Since we saw in Sect. 3.3 that for the high tan β region the CP-odd Higgs boson can also provide an explanation for the excess in the Zh final state, we show the corresponding signal cross sections for the points of this scan in Fig. 10. As before, the subset of points that are in the wrong-sign Yukawa coupling regime provide relatively large contributions to the excess, as can be seen in the right plot of Fig. 10, while the other parameter points do not give rise to a signal in the A → Zh channel. Overall, the values for the cross sections are slightly smaller here than what was found in Sect. 3.3. The reason is that here A has an additional decay mode to a Z boson and 200 300 400 500 600 700 800 900 1000 Type II: Expected ±1σ observed 9 10 11 tan β 200 300 400 500 600 700 800 900 1000 Type II:

NMSSM interpretation
Our analysis in the previous section has shown that the N2HDM of type II provides an attractive framework for accommodating either the tt or the τ + τ − excess at 400 GeV simultaneously with the observed excesses at 96 GeV. In SUSY models, due to the holomorphy of the superpotential, the presence of two Higgs doublet fields with opposite hypercharge, usually denoted as H u and H d , is required for the mass generation of both up-and down-type fermions. Identifying H d = Φ 1 and H u = Φ * 2 , there is a correspondence between the Yukawa structure of the (N)2HDM type II and SUSY extensions of the SM. Thus, it is an interesting question whether the excesses can also be realized in SUSY models. It should be noted in this context that SUSY yields further restrictions on the Higgs sector compared to non-supersymmetric models with a similar structure of the extended Higgs sector.
The simplest SUSY extension of the SM is the Minimal Supersymmetric Standard Model (MSSM) [67,68]. The Higgs sector of the MSSM comprises only the two Higgs doublet fields H u,d introduced above. In the MSSM the quartic scalar couplings of the Higgs sector are at lowest order functions of the gauge couplings due to SUSY relations. The rigid structure of the scalar potential leads to the fact that deviations w.r.t. the SM of the couplings of the lightest Higgs boson h 1 = h 125 scale with the inverse of the mass parameter M A [69] of the CP-odd Higgs boson. Consequently, it is possible to suppress such deviations in the decoupling limit M A M Z , in which h 125 is aligned with the SM vev and effectively indistinguishable from a SM Higgs boson, whereas for M A 1 TeV deviations at the level of 1-10% for the coupling coefficients c h 1 ff are present [70]. 9 Hence, the signal-rate measurements of h 125 at the LHC, which at the present level of accuracy show no significant deviations from the SM prediction, can be cast into a lower limit on M A . Taking into account the most recent LHC results [74,75], the current limit was found to be M A 600 GeV, (largely) independently of the value of tan β [73,76]. Consequently, the MSSM cannot account for the presence of a CP-odd Higgs boson state around 400 GeV without being in strong tension with the signal rates of h 125 , and a SUSY interpretation of the observed excesses therefore needs to be based on non-minimal SUSY models.
The Next-to Minimal Supersymmetric Standard Model (NMSSM) [77,78] extends the Higgs sector of the MSSM by a complex singlet scalar field. The presence of this singlet field is motivated in particular by the so-called µ-problem of the MSSM [79], in which a dimensionful parameter µ is present in the superpotential. While this parameter would naturally be expected to be of the order of the ultraviolet cutoff of the model, it must be of the order of the EW scale to allow for a phenomenologically viable Higgs sector. In the Z 3 symmetric NMSSM this problem is absent, because the µ term is forbidden by the global symmetry. Instead, it is generated dynamically when the scalar component of the singlet superfield acquires a vev, and the Z 3 symmetry is broken spontaneously (see Ref. [78] for details). Besides the fact that the singlet field is complex in the NMSSM, leading to the presence of a second CP-odd Higgs boson in addition to the CPeven scalar singlet state, and that the singlet field is charged under a Z 3 and not a Z 2 symmetry, the Higgs sector and the Yukawa sector resemble the one of the N2HDM type II introduced in Sect. 3.1. Thus, from the point of view of the Higgs phenomenology, the NMSSM is a promising model for investigating the possible realization of the excesses around 400 GeV in the context of SUSY.
We briefly summarize some of the key features of the NMSSM that are relevant for our analysis. We will focus in particular on the Higgs sector of the model. The superpotential of the Z 3 symmetric NMSSM is given by where W MSSM, ¡ µ denotes the superpotential of the MSSM except the aforementioned µ-term. The second term is the portal coupling between the gauge singlet superfieldŝ and the Higgs doublet superfieldsĤ u,d (superfields are denoted here by a hat). The (effective) µ-term arises from this term when the scalar component ofŝ acquires the vev s = v S , µ := λv S . The third term proportional to κ is the singlet self-coupling and, again once v S is non-zero, gives rise to bilinear mass terms both for the scalar and the fermionic component ofŝ. The complete Higgs potential is then derived from the usual F -and D-terms in combination with the terms arising from the soft SUSY-breaking Lagrangian.
The fermionic superpartners of the Higgs fields (called Higgsinos and singlino), together with the superpartners of the gauge bosons (called gauginos) give rise to five massive neutral fermion states (called neutralinos) and two massive charged fermions (called charginos). The gauginos obtain their masses via the soft SUSY-breaking gaugino mass parameters M 1 , M 2 and M 3 , while the higgsinos obtain their masses from the µ parameter. Compared to the MSSM, the only additional fermion is the fifth neutralino (singlino). Moreover, extending the matter sector of the SM, the NMSSM contains the scalar partners of the leptons (called sleptons and sneutrinos) and quarks (called squarks).
The superpotential as defined in Eq. (9) conserves R parity [80]. As a result, the lightest SUSY particle (LSP) is stable. In case the LSP forms a sizable part of the dark matter abundance, bounds from direct detection experiments exclude a substantial part of the NMSSM parameter space. Since we are primarily interested here in the phenomenology of the Higgs sector, we will not include these limits in our analysis. Accordingly, the NMSSM can be considered to be a low-energy effective model of a more complete theory in which R parity violating effects prohibit the existence of a stable SUSY particle, but in which the Higgs phenomenology remains effectively unchanged. This can be realized if, for instance, the R parity violating operators are non-renormalizable and suppressed by the inverse of a high energy scale such as the Planck scale, or if another symmetry suppresses the amount of R parity breaking. An example for the latter is the so-called µ-fromν supersymmetric standard model (µνSSM) [81,82] (see Ref. [83] for a recent review). It was shown that besides the fact that there are more than one gauge singlet scalar fields present in the µνSSM, the phenomenology of the SM-like Higgs boson can be practically unchanged w.r.t. the NMSSM [26,27,84].

The Higgs sector: alignment without decoupling
In the NMSSM an alignment limit exists in which no decoupling of the BSM Higgs bosons and no large cancellations between lowest-order contributions and radiative corrections to scalar masses and their couplings are required [30]. We briefly describe here the necessary conditions that have to be fulfilled in order to account for the presence of a Higgs boson resembling the discovered particle state at 125 GeV even for M A ≈ 400 GeV. We will also demonstrate that only for low values of tan β the alignment conditions can be satisfied. Thus, we expect that a description of the tt excess can be realized in the exact alignment limit of the NMSSM, while a description of the τ + τ − excess will require departures from the alignment limit.
In the following, we will denote the CP-even Higgs bosons either with h 1,2,3 , which is the massordered notation, such that m h 1 < m h 2 < m h 3 , or we will use h 125 for the discovered Higgs boson, h S for the Higgs boson with dominant singlet component and H for the "heavy" CP-even Higgs boson. The latter notation is useful because in the alignment limit the singlet field does not mix with h 125 . For the two CP-odd Higgs bosons we either use the mass-ordered notation A 1,2 or, in order to make connection to the notation of the MSSM, we write A S for the singlet-like state and A for the MSSM-like CP-odd Higgs boson. As for the CP-even scalars, the latter notation does not imply a mass ordering of the particles.
In the alignment limit one of the doublet fields is aligned in field space with the Higgs doublet vev v (with v 2 = v 2 1 + v 2 2 ). Then the tree-level couplings of the corresponding Higgs state, h 125 , acquire their SM values. This can be achieved if the mass matrix of the Higgs bosons is diagonal in the Higgs basis, i.e., after a rotation of the 2 × 2 submatrix of the doublet components by the angle β. In order to ensure that the non-diagonal entry of the mass matrix between the two doublet fields vanishes, as a first condition the following relation has to be fulfilled [30], This condition was derived taking into account the dominant one-loop corrections to the diagonal mass term of h 125 , stemming from the (s)top-sector. As long as tan β is not much larger than 10, the other one-loop corrections (entering also the non-diagonal mass terms) are suppressed by a factor µ/M S . Since we will choose M S µ in our numerical discussion, they can safely be neglected. The same is true also for corrections beyond the one-loop level in the considered parameter region, where the prediction for the mass of the SM-like Higgs boson is largely dominated by the tree-level contribution (see the discussion below). Given the values for the SM vev v ≈ 246 GeV, the mass of the Z boson M Z ≈ 91 GeV and m h 125 ≈ 125 GeV, the first alignment condition determines λ as a function of tan β. A second alignment condition arises from the fact that also the nondiagonal mass matrix entry relating the state h 125 and the gauge singlet field s has to vanish in the alignment limit. This translates into the requirement [30] M 2 A sin 2 2β 4µ 2 + κ sin 2β 2λ = 1 .
Here, M 2 A is the squared CP-odd mass parameter defined by where A λ is the soft SUSY-breaking trilinear coupling corresponding to the λ term in the superpotential of Eq. (9). At tree level M A is equal to the mass of the MSSM-like CP-odd Higgs boson A, so that it is a useful input parameter for our analysis. In order to make a distinction between M A and the loop-corrected physical masses of the CP-odd Higgs bosons, we will denote the latter with the lower-case letter m in the following, i.e. as m A 1 , m A 2 or m A S , m A . Since M A will be fixed to M A ≈ 400 GeV in order to account for the tt or the τ + τ − excess, and taking into account that the value of λ is given by Eq. (10) for a fixed value of tan β, the second alignment condition can be regarded as a relation for κ in terms of µ and tan β.
In the left plot of Fig. 11 we show the values of λ as a function of tan β obtained from Eq. (10). One can see that relatively large values of λ > 0.6 are required. The two black dashed vertical lines indicate the tan β values that we will choose as benchmark scenarios to realize the tt or the τ + τ − excesses in our numerical discussion in Sect. 4.3 and Sect. 4.4, respectively. Using the corresponding values for λ, we show in the right plot of Fig. 11 κ as a function of µ for the two values of tan β. We also indicate with a gray shaded region the values of µ 104 GeV for which the mass of the lightest chargino is below the limits that were obtained from chargino searches at LEP [85]. Here we have assumed that M 2 is sufficiently larger than µ, such that the light chargino mass is determined by µ. One can see that for tan β = 8 (orange line) the second alignment condition can only be satisfied for rather large values of the singlet self-coupling κ. The required values are roughly an order of magnitude larger than the maximum value κ max GUT (indicated by the orange dashed line) based on the condition that there should be no Landau poles below the GUT scale, approximately given by [86] κ 2 + λ 2 ≤ 0.7 2 .
We thus conclude that for values of tan β required to realize the τ + τ − excess one cannot satisfy the alignment conditions for κ values that are in agreement with Eq. (13). On the other hand, for tan β = 1.7, which is the value used for the benchmark scenario addressing the tt excess (see Sect  11) for tan β = 1.7 (blue) and tan β = 8.0 (orange), and assuming that λ is chosen such that it satisfies the first alignment condition. The gray shaded region is excluded by LEP searches for charginos [85]. The dashed lines indicate the maximum values of κ satisfying the condition that there are no Landau poles below the GUT scale given in Eq. (13).
a region of parameter space in the alignment limit can be explored that is in agreement with the LEP limits and the perturbativity constraint of Eq. (13). For the description of the τ + τ − excess at tan β = 8 in Sect. 4.4 we will allow for a departure from the alignment limit in order to satisfy Eq. (13).
Focusing on the low tan β regime and thus the tt excess, one can easily estimate the structure of the Higgs boson spectrum. Obviously, one has to require M A ≈ 400 GeV such that the CP-odd Higgs boson A can play the role of the particle explaining the tt excess. Moreover, the tree-level mass of the charged Higgs bosons H ± is related to M A via the relation For the considered values of λ, the charged Higgs bosons are a few GeV lighter than the CP-odd state A. The CP-even doublet state H is also close in mass to the A boson as long as its singlet component is small. For the squared mass of the SM-like Higgs boson at tree level an upper bound is obtained that is approximately given by In the alignment limit this inequality is nearly exhausted, and values of m 2 h 125 ≈ 125 GeV or even slighly larger are obtained for small tan β and large λ. Thus, no large radiative corrections to m h 125 are required in order to obtain a physical mass of 125 GeV. Finally, the masses of the singletlike particle states h S and A S are controlled by κ and the corresponding soft trilinear parameter A κ . Due to the small values of κ in order to satisfy Eq. (11) (see also Fig. 11), these two Higgs bosons have masses substantially below M A . The dependence of their masses on A κ is opposite and approximately given by where for simplicity possible mixing contributions have not been spelled out in Eq. (16) (those mixing contributions are included in our numerical analysis below). Thus, depending on the value of A κ both h S and A S can be lighter than 125 GeV/2, giving rise to both a lower and an upper limit on A κ in order to avoid large values of BR(h 125 → h S h S , A S A S ). Moreover, the singlino mass term is given by

Experimental constraints
The experimental constraints that we take into account in our NMSSM analysis are focused on the collider phenomenology. In particular, as discussed above, we do not intend to reproduce the observed dark matter relic abundance, nor do we apply constraints from dark matter direct detection experiments. Furthermore, we do not exclude parameter points based on constrains from flavor observables, as the theoretical predictions of these observables in SUSY models depend on various different sectors of the model, while the focus of our analysis is the Higgs sector phenomenology. The impact of constraints of this kind on the parameter space investigated in our analysis would depend on the parameter settings of other BSM contributions that are not relevant for Higgs physics. While we do not impose those constraints for excluding parameter points, in our numerical discussion we will point out possible tensions with constraints from observables beyond the Higgs sector. Tensions with observables from the flavor sector are expected in particular for relatively small values of m H ± M A and low values of tan β [87], see the discussion above. It should be noted that we treat the mentioned upper limits on possible BSM effects from other sectors in the same way as the 4.2 σ discrepancy between the experimental results for the anomalous magnetic moment of the muon (g − 2) µ [88,89] and the SM prediction [90]. While in principle SUSY models are capable of explaining this discrepancy [91,92] (see also Ref. [93] for an analysis of the EW MSSM sector, Ref. [94] for a specific NMSSM analysis, Ref. [95] for an analysis in the above discussed µνSSM, and Ref. [96] for a recent review), since the Higgs-boson sector is only marginally involved we also do not apply this constraint favoring a non-zero BSM contribution.
As done for the N2HDM analysis, we confronted all parameter points with the current limits from searches for additional Higgs bosons using HiggsBounds. There are several collider searches which are especially relevant for the low-M A region of the NMSSM. We briefly summarize the most important ones here: (i) For low values of tan β 2, the LHC searches for H ± decaying into a tb pair exclude parameter regions with m A ≈ m H ± ≈ 400 GeV in the (N)2HDM type II [97]. In the alignment limit of the NMSSM, in which such values of tan β allow for the realization of the tt excess, one finds m H ± < M A because of the relatively large values of λ (see Eq. (14)). Thus, compatibility with the experimental limits on the charged Higgs rates requires a reduction of BR(H ± → tb) compared to the (N)2HDM prediction which can occur if additional decay channels into SUSY particles are kinematically open. Nevertheless, even for a considerably suppressed BR(H ± → tb) the LHC charged Higgs boson searches [97][98][99] belong to the most constraining BSM Higgs boson searches for the analysis dedicated to the tt excess.  [19].
(iv) The presence of rather light singlet states h S and A S , in particular as predicted in the alignment limit, opens up the possibility of Higgs cascade decays H(A) → A S (h S ) Z and H ± → h S /A S W ± . CMS and ATLAS searched for the former processes in multilepton and llbb final states [100][101][102]. In our analysis, only the constraints from the searches regarding the signature H(A) → A S (h S ) Z were capable of excluding parameter points. On the hand, for our analysis the presence of the decay mode H ± → h S /A S W ± turned out to be relevant, because additional H ± decay channels into BSM scalars have the potential to further reduce BR(H ± → tb) and help to avoid the constraints from charged Higgs boson searches in the tb final state discussed above. The experimental signatures mentioned above were also investigated under the assumption that the Higgs boson in the final state is h 125 , giving rise to exclusions of parts of the here analyzed parameter space [9]. This happens when a second BSM Higgs boson h S or A S is present in the mass window (125 ± 10) GeV, yielding another contribution that has to be added to the one of h 125 , such that there can be a relevant enhancement of the experimental signature. This case will be discussed in more detail below. Finally, also purely scalar cascade decays of the form H(A) → h 125 h S (A S ) are relevant, where CMS searched for such a signature assuming that h 125 decays into a τ + τ − pair and the other (BSM) final state particle decays into a bb pair [103].
In addition to the presence of additional Higgs bosons, also SUSY particles can be in the reach of the LHC or previous colliders. The corresponding search limits yield further constraints on the parameter space. We list below the most relevant searches and their impact on our analysis: (i) In our analysis in the alignment limit, the tree-level mass of m (0) h 125 ≈ 125 GeV is determined by the alignment conditions (see Sect. 4.1) for low tan β values, which places a limit on the size of the radiative corrections to the mass of the SM-like Higgs boson. As a consequence, the SUSY breaking scale M S , and therefore the stop masses m t , cannot be too large (see Ref. [30] for more details). This implies that the LHC stop searches can be relevant, as they provide (under several assumptions) a lower limit on m t . Taking into account the uncertainties on the predictions for the Higgs-boson mass(es) in the NMSSM [104][105][106], we use an interval of m h 125 = (125 ± 4) GeV in our analysis. In order to not exceed this mass interval given the value of tan β = 1.7 used in the analysis regarding the tt excess in Sect. 4.3, we find that values of M S ≈ m t 1.2 TeV are required, which is of the order of the current experimental lower limit on the stop masses in simplified scenarios [107]. Here it should be noted that in our analysis, due to the compressed electroweakino spectrum and the presence of Higgs bosons with masses much below M S , the actual experimental lower limit on the stop masses could be substantially weaker than 1 TeV. For larger values of tan β, as required for the τ + τ − excess discussed in Sect. 4.4, the tree-level enhancement of m h 125 given by the term ∼ λ 2 (see Eq. (15)) is suppressed, such that in this case larger radiative corrections are needed to achieve a value of m h 125 ≈ 125 GeV. Consequently, for the τ + τ − excess larger values of M S are required and LHC searches for stops have no significant impact.
(ii) As already mentioned in Sect. 4.1, LEP searches for charginos yield lower limits on the chargino masses of up to m χ ± 1 > 104 GeV [85]. In order to ensure compatibility with the LEP bounds, we demand m χ ± 1 > 104 GeV in our analysis. (iii) LHC searches for light neutralinos in the context of the NMSSM are very challenging when the chargino and neutralino spectra are compressed. In addition, the searches for the neutralinos suffer from the fact that also the background estimation depends on the precise form of the spectrum of the SUSY particles. Therefore, a model interpretation of a particular experimental search has to be done not only for each SUSY model separately, but also for each parameter point within a certain model. Besides, neutralino searches critically depend on whether R parity is assumed to be conserved, or not. Taking the above mentioned considerations into account, we cannot apply general exclusions on the parameter space of the NMSSM from neutralino searches in our analysis. The only exception are exclusion bounds based on the presence of a neutralino with a mass below m h 125 /2, since this can give rise to large values of BR(h 125 → χ 0 1 χ 0 1 ) (see Sect. 4.3), which are excluded by global constraints on the signal rates of the h 125 .
Finally, as we did in the N2HDM analysis, we perform a χ 2 test regarding the signal rates of h 125 using HiggsSignals (which as discussed above excludes large values of BR(h 125 → χ 0 1 χ 0 1 )). An important difference w.r.t. the N2HDM arises from the fact that the Higgs-boson mass could be chosen as an input parameter in the N2HDM, whereas m h 125 is predicted in the NMSSM as a function of other model parameters. The model predictions are affected by a theoretical uncertainty in particular for larger values of λ due to radiative corrections beyond the one-loop level (see Ref. [108] for a recent account of this subject). This is why, as already mentioned above, we allow for an uncertainty of ∆m theo h 125 = 4 GeV in our analysis. In order to make sure that the mass uncertainty does not give rise to unacceptably large modifications of the predicted branching ratios for the decays of h 125 into SM particles (which could occur if the branching rations are evaluated with mass values that are several GeV away from 125 GeV, but still within the ±4 GeV interval), we recalculated the corresponding decay widths by requiring m h 125 = 125 GeV for all points in which the prediction of NMSSMTools for m h 125 lies within the mentioned uncertainty band. This recalculation is especially relevant for BR(h 125 → W W * , ZZ * ), which have a very sensitive dependence on m h 125 . For the recalculation of the branching ratios we made use of the effective coupling coefficients provided by NMSSMTools, by which the partial decay widths for h 125 as predicted by the SM were rescaled. The values for the SM decay widths were taken from Ref. [109].

Benchmark scenario with tan β = 1.7
As discussed in Sect. 4.1, the alignment-without-decoupling limit of the NMSSM is a theoretically well motivated scenario for realizing the tt excess in the context of SUSY. Thus, we will investigate in this section the possibility of accommodating the tt excess in this limit taking into account the various different collider constraints mentioned in the previous section. As in the N2HDM analysis, we consider only the CP-odd Higgs boson A as the origin of the excess, accounting for the fact that the contribution of CP-even states is much smaller compared to the one of a CP-odd state. Moreover, since we restrict our discussion to the CP-conserving case the contributions of CP-even and CP-odd states do not interfere with each other. In addition, we pointed out that in this limit the presence of a second light CP-even Higgs boson arises naturally. Thus, although not being the main focus of this analysis, for the subset of points that feature a particle candidate in the relevant mass range for the LEP and the CMS excesses (see Sect. 2.4) we will check in a second step whether also a simultaneous realization of these excesses in combination with the tt excess is possible.
As can be seen from the first alignment condition shown in Eq. (10), the choice of the parameter tan β plays a key role. Starting from the value of tan β, and given M A ≈ 400 GeV, the remaining parameter values can either be derived or scanned over. We begin by choosing an initial value of tan β = 1.7 for our scan. This value arises from the following considerations. As described in Sect. 4.1, in the alignment limit the singlet CP-odd state A S is substantially lighter than the doublet state A. Accordingly, also the mixing of these states is small. Then one finds that c Att ≈ 1/ tan β (as in the N2HDM type II), such that regarding the tt excess even smaller values of tan β ≈ 1 would be favored. However, the presence of the charged Higgs bosons H ± , having a mass slightly below M A in the alignment limit, yields a lower bound on the possible values of tan β based on LHC searches for charged scalars (see Sect. 4.2). Our numerical analysis has revealed that, by taking into account the possibility that H ± can partially decay into SUSY particles and lighter neutral Higgs bosons plus a W boson, values of tan β = 1.7 provide points that are compatible with the constraints arising from the LHC searches, while also allowing for sizable values of c Att .
Given the chosen value for tan β, one derives a value of λ = 0.6617 from Eq. (10). Taking into account that loop corrections yield a (physical) mass m A of the CP-odd Higgs boson that is slightly below the (tree-level) input parameter M A , we chose to generate parameter points with 410 GeV ≤ M A ≤ 430 GeV, in steps of 5 GeV. Given a value for M A and the ones for tan β and λ, we used the second alignment condition shown in Eq. (11) to obtain values of κ as a function of µ. The parameter points were then generated by scanning over µ and A κ in steps of 1 GeV. We covered the parameter ranges for which points could be found that are in agreement with the experimental and theoretical constraints. The range of µ has a lower limit (for each value of M A ) arising from the restriction that none of the neutralinos should become substantially lighter than 125/2 GeV, since otherwise the decay h 125 → χ 0 1 χ 0 1 would become kinematically allowed, spoiling agreement of the properties of h 125 with the LHC measurements. An upper limit on µ is given by the aforementioned constraints from the charged Higgs boson searches. These constraints can only be evaded if there are sizable branching ratios for the decays H ± → χ 0 1 χ ± 1 , h S W ± and/or A S W ± . However, the branching ratios of these decays get reduced for increasing values of µ, as the masses of the final state particles, in particular the mass of the Higgsino-like chargino χ ± 1 , increase with µ. Upper and lower limits on A κ , in dependence of the other parameters, are given by the fact that either h S or A S becomes lighter than 125/2 GeV (see Eq. (16)). In this case the large value  of λ gives rise to unacceptably large values of BR(h 125 → h S h S /A S A S ). In addition, we find that all points with A S 100 GeV are excluded due to constraints from the search for the signature A → A S h 125 [103].
We summarize the parameter ranges of the scan in Tab For this choice several decays of the kind H ± → χ 0 i χ ± j can become relevant, which play a role in avoiding the constraints from charged Higgs boson searches. As discussed above, we checked that the lighter chargino mass is above the lower limit from LEP constraints. The value for the gluino mass M 3 is large enough to be compatible with the limits from direct searches given the compressed neutralino-chargino spectrum, independently of the fact whether R parity is assumed to be conserved or not [110].
Taking into account the experimental constraints listed in Sect. 4.2, we found an allowed parameter region that is shown in Fig. 12. The left plot shows the µ-κ plane, where the color coding indicates M A , whereas the right plot shows the κ-A κ plane, where the color coding indicates m h S . Here, we applied the same condition as in the N2HDM analysis, demanding that χ 2 ≤ χ 2 SM , where the total χ 2 is defined in Eq. (8). It should be noted that the contribution of χ 2 96 is included in χ 2 , even though it is not the main focus of this analysis. For parameter points that feature a CP-even scalar in the mass interval 94 GeV ≤ m h 1 ≤ 98 GeV, we calculate χ 2 96 given the predicted signal strengths of h 1 , whereas for the other points we set χ 2 96 = χ 2 SM, 96 . Thus, only a subset of points has a contribution of χ 2 96 below the SM one. However, even though the formal definition of χ 2 is identical, one should keep in mind that in the NMSSM analysis the physical mass m A of the CP-odd Higgs boson is not an input parameter, such that it can differ slightly from the value m A = 400 GeV used in the N2HDM. This leads to the fact that in the NMSSM χ 2 tt is a function of the three predicted quantities c Att , Γ A and m A (see also Sect. 2). As before, the result for the SM is given by χ 2 SM,tt = 13.98. As a consequence of the alignment conditions, we found that the properties of h 125 are very well in agreement with the experimental measurements. Even parameter points in which the decay of h 125 into two neutralinos is kinematically allowed were found to predict χ 2 ≤ χ 2 SM . A more detailed discussion of the properties of h 125 can be found in App. B.
In the left plot of Fig. 12 one can see that the allowed range of κ is below the limit from the GUT condition shown in Eq. (13) in most cases. This indicates that the considered scenarios correspond to model realizations that may be well defined up to the GUT scale. The parameter range of µ is close to the EW scale, which corresponds to the parameter region of the NMSSM that is favored by naturalness arguments [111,112], in which no large fine tuning (corresponding to a "little hierarchy problem") is required in the EW sector of the model. Concerning the colors of the points, indicating M A , one should note that each point in the plot corresponds to a subset of parameter points with different values of A κ . One can see that for increasing values of M A also larger values of µ and κ are required in order to fulfill the second alignment condition.
In the right plot of Fig. 12 one can see that only points with negative values for A κ are allowed. For values of A κ > −100 GeV we find that the condition χ 2 < χ 2 SM is only fulfilled for a small number of points, which then are excluded by LHC searches involving the light CP-odd Higgs boson A S with m A S 148 GeV. In particular, the CMS search for the signature A → A S h 125 excludes all such points with m A S 100 GeV, and points with slightly larger values of m A S are excluded by searches for H → A S Z (see the discussion in Sect. 4.2). On the other hand, m 2 h S receives additional contributions proportional to κv S that can compensate the negative contribution ∼ A κ , so that m 2 h S > 0 even for values of A κ ≈ −500 GeV (see Eq. (16)). Nevertheless, one can see that, given a certain value for κ, a lower bound on A κ is found where m h S (indicated by the color coding) decreases below ≈ 125/2 GeV, giving rise to unacceptably large values of BR(h 125 → h S h S ). In addition, even though A S is overall heavier if A κ has large negative values, for the lower values of κ it can still be relatively light, with masses just above or even below 125 GeV. Hence, for small values of κ the searches for A S already mentioned above become relevant again, and give rise to the fact that none of the points with 0.05 ≤ κ ≤ 0.07 pass the HiggsBounds test. Also the two   remaining isolated lines of points at the lower end of the κ range, which correspond to the two points at the lowest κ values in the left plot of Fig. 12, barely escape the constraints from the searches for H → A S Z. As discussed above, points with m h S ≈ 96 GeV can contribute to the LEP and CMS excesses. These points are found for larger κ and smaller |A κ |. Given the condition χ 2 ≤ χ 2 SM in our scans, for these points somewhat larger values of ∆χ 2 125 are allowed as compared to the points for which no candidate at 96 GeV is present. 11 We now turn to the discussion of the description of the tt excess, which is the main motivation for the scan presented here. On the left-hand side of Fig. 13 we show the value of c Att of the parameter points in dependence of the (physical) mass m A . The color coding shows Γ A /m A in %. Since the singlet admixture of A is very small, c Att is given to a very good approximation by the inverse of tan β, such that c Att ≈ 1/1.7 ≈ 0.59. Consequently, the points are all located within a tight branch at roughly this value of c Att , slightly below the experimental best-fit value of c Att ≈ 0.8 assuming m A = 400 GeV and Γ A /m A = 2.5%. This results in a description of the excess at a confidence level of about 1.5 σ for the smallest values of m A , as further discussed below. 12 The mass of the CP-odd Higgs boson A ranges between 400 GeV and 430 GeV, where the range is defined by the scan range used for the tree-level parameter M A . In order to facilitate the comparison between the N2HDM and the NMSSM we also show the subset of N2HDM parameter points featuring 2.0% ≤ Γ A /m A ≤ 3.0% from the scan performed in Sect. 3.4. The N2HDM 11 It should be remembered that due to the uncertainty in the prediction of the Higgs-boson masses, we calculate χ 2 96 as described in Sect. 2.4 for points with 94 GeV ≤ m h S ≤ 98 GeV, while for the other points we set χ 2 96 = χ 2 SM,96 . 12 We wish to stress that the fact that the NMSSM parameter points happen to lie along the curve for the expected limit is a coincidence from which no direct information about the fit result can be derived. points lie at larger values of c Att ≈ 0.8, reflecting the fact that the N2HDM allows for an even better description of the excess. The reason for this is that one cannot further reduce the value of tan β in the NMSSM without violating constraints from the charged Higgs boson searches. In the N2HDM these can be avoided by increasing the mass m H ± 400 GeV. In the NMSSM this is not possible, because of the relations among the Higgs-boson masses, as explained in Sect. 4.2. Thus, also the maximum value of |c Att | that is possibly realized remains below the measured best-fit value (c Att ≈ 0.8). The collider phenomenology of the H ± state in the NMSSM will be discussed in more detail below.
It is also important to note that the tt excess is most pronounced for the smallest values of m A analyzed in the experimental analysis. In Fig. 13 this is reflected by the fact that the smaller m A , the larger is the difference between the observed (blue) and expected (black dashed) exclusion limit. It should be remembered here that χ 2 tt is defined as the χ 2 difference compared to the best fit value, where the latter was found for m A = 400 GeV, Γ A /m A = 4.5% and c Att = 0.94 (see also Fig. 1). Thus, even though the theoretical predictions for the values of c Att are independent of m A , the improvement of the χ 2 tt values in comparison to χ 2 SM,tt is substantially larger for values of m A ≈ 400 GeV, where the smallest χ 2 tt values found are χ 2 tt ≈ 5, i.e. somewhat worse than in the N2HDM (see above). This is clearly visible in the right plot of Fig. 13, in which we show the values of χ 2 tt in dependence of m A : the smallest values of χ 2 tt ≈ 5 (corresponding to a fit at about 1.4 σ confidence level) are obtained for the points with the lowest values of m A . 13 Thus, even though the NMSSM offers a substantially better fit to the data than the SM, the NMSSM analysis does not reach the same level of fit quality as the N2HDM analysis, where values of χ 2 tt < 1 could be achieved (see, for instance, Fig. 3).
In addition, the total width Γ A of the CP-odd Higgs boson, indicated by the color coding in units of m A in Fig. 13, correlates with m A , where for each branch displayed in the right plot of Fig. 13 (each branch corresponds to a different value of M A ) smaller values of m A give rise to smaller values of Γ/m A . This correlation has its origin in the phase space factor for the tt decay width, which grows with increasing values of m A and overcompensates the factor m A in the denominator. Since the experimental analysis was carried out for different hypotheses on Γ A /m A , χ 2 tt depends on m A both directly and via its impact on the ratio Γ A /m A . On the other hand, one can see that the overall variation of Γ A /m A of the parameter points is small compared to the step sizes of the different width hypothesis used in the experimental analysis, as shown in Fig. 1. Thus, given the experimental uncertainties, it would have been sufficient to only take into account the expected and observed exclusion limits (and the resulting χ 2 tt values) obtained under the assumption Γ A /m A = 2.5% (this is the reason why the left plot of Fig. 13 is displayed for Γ A /m A = 2.5%).
As discussed at the beginning of this section, the alignment conditions lead to the presence of a light CP-even Higgs boson h S (= h 1 ), which is almost entirely singlet-like. As a result, h S can be a candidate for describing the excesses at ≈ 96 GeV. We show in Fig. 14 the parameter points that fulfill the additional constraint 94 GeV ≤ m h S ≤ 98 GeV in the plane of µ CMS and µ LEP , where the color coding indicates m h S . We calculate the signal strengths via hence making use of the fact that the cross section ratios can be approximated via the respective effective coupling coefficients squared. One can see that the CMS excess can easily be accommodated. On the other hand, in our analysis the LEP excess cannot be described. This was expected as a consequence of the dominantly singlet-like character of h S and the resulting suppression of the Higgsstrahlung production cross section. 14 Regarding the CMS excess, the obtained relatively large values of 0.35 < µ CMS < 0.6 have their origin in an enhancement of the diphoton branching ratio of the state h S , which reaches values that are an order of magnitude higher than the SM prediction. This enhancement compensates the suppression of the gluon fusion production cross section that is approximately given by c 2 h S gg ≈ 0.07. The large values of BR(h S → γγ) are caused by two different contributions. Firstly, the light chargino, which is close in mass to h S , provides a (positive) BSM contribution to c h S γγ . Even more important, however, are the patterns of the effective couplings of h S to top and bottom quarks. In our scan we find for the coupling coefficients 0.220 < |c h S tt | < 0.237 and 0.074 < |c h S bb | < 0.105. The dominant component for the decay width of the diphoton decay is the diagram with top quarks in the loop. Thus, this partial decay width scales with c 2 h S tt . The dominant component of the total width of h S is given by Γ(h S → bb), such that the total width is approximately proportional to c 2 h S bb . As a result, we observe an enhancement of BR(h S → γγ) by factors of 4.5 < c 2 h S tt /c 2 h S bb < 7.3, which can be further increased by a possible chargino contribution. Accordingly, the CMS excess can be described with an almost singlet-like Higgs boson.
Besides the candidates for describing the discussed excesses, there is another very prominent feature of the scenario presented here that offers promising possibilities for future collider searches. As was already mentioned in Sect. 4.2, in the alignment limit the CP-odd Higgs boson A at ≈ 400 GeV is accompanied by slightly lighter charged Higgs bosons H ± . Due to the small value of tan β, the coupling of H ± to top quarks is relatively large. Experimental upper limits were obtaind . Also shown are the current observed (black) and expected (black dashed) upper limits at the 95% C.L. and the expected 1σ (green) and 2σ (yellow) exclusion regions obtained by ATLAS [97]. The red line indicates the previous ATLAS result for the observed upper limit using 35 fb −1 [98].
on the product σ(pp → tbH ± ) × BR(H ± → tb) from LHC searches performed both by ATLAS [97] and CMS [113]. The currently most stringent limits for the mass region relevant for our scenario was reported by ATLAS in Ref. [97], using the full 13 TeV Run 2 data set with 139 fb −1 . We show in Fig. 15 the parameter points of our scan with the corresponding cross section times branching ratio on the vertical axis and the charged Higgs boson mass on the horizontal axis. We also include the ATLAS observed and expected upper limits at the 95% C.L. Since these limits are included in the most recent version of HiggsBounds, which was used to check for the constraints from direct searches of BSM Higgs bosons, our set of parameter points lies below the observed limit.
As discussed in Sect. 4.2, those experimental limits could be evaded if decay modes of H ± to SUSY particles are kinematically open. In Fig. 15 the color of the points indicates the values of the lightest neutralino mass, m χ 0 1 . One can observe that the signal rate of H ± decreases with decreasing m χ 0 1 . This indicates that in particular the decay modes H ± → χ 0 1 χ ± 1,2 play an important role in this context, weakening the collider bounds from the searches ultilizing the tb final state. 15 Nevertheless, even with such a reduction of the tb mode our analysis shows that the signal rate is still rather close to the current observed upper limit. As a consequence, there are good prospects that future charged Higgs boson searches at the LHC and the HL-LHC, both in the tb final state and via dedicated searches exploring decays into final states involving BSM particles, will be able to probe this scenario.
Another possiblity for testing the scenario presented here is given by the precise measurements of low energy observables. As was mentioned in Sect. 4.2, these are particularly sensitive to the presence of the relatively light charged Higgs bosons and electroweakinos. We used NMSSMTools to obtain predictions for flavour physics observables, including their theoretical uncertainties, and the anomalous magnetic moment of the muon. We find that most of the flavour observables are  predicted to be in agreement with the experimentally measured values, despite the presence of H ± below 400 GeV. However, there are a few observables that show deviations at the ≈ 2σ level. As an illustrative example, we briefly state the largest discrepancies for the best-fit point of this scan. Similar values are found for the other points of the scan. For the BR(b → sγ) decay, the predicted range for the best-fit point within the estimated theoretical uncertainties is 3.67 · 10 −4 ≤ BR theo (b → sγ) ≤ 4.82 · 10 −4 [114][115][116][117][118][119][120][121], which lies just above the experimental 2σ interval 3.02 · 10 −4 ≤ BR exp (b → sγ) ≤ 3.62 · 10 −4 [122]. A similar discrepancy, however this time in the opposite direction, is found for the decay B d → µ + µ − , where the predicted range 1.34·10 −11 ≤ BR theo (B d → µμ) ≤ 9.15·10 −11 [114,[123][124][125] is ≈ 2σ smaller than the experimental range 11 · 10 −11 ≤ BR exp (B d → µμ) ≤ 71 · 10 −11 [122]. Here the substantially larger experimental uncertainties should be kept in mind.
Summarizing the results of the scan presented here, we find that the alignment-withoutdecoupling limit of the NMSSM is a theoretically well motivated framework that is capable of describing, at least approximately, the observed tt excess in the context of SUSY. In addition, by choosing the singlet scalar self-coupling A κ appropriately, a simultaenous description of the CMS excess at about 96 GeV is possible (but not of the LEP excess). We have found that the latter scenario is compatible with the experimental results for the properties of the h 125 state even for the case where the decay of h 125 into BSM particles, for instance into two neutralinos, is kinematically open.

Benchmark scenario with tan β = 8
In this section we analyze a parameter region in the NMSSM aiming towards a possible description of the τ + τ − excess. As already discussed in Sect. 4.1, this requires larger values of tan β as compared to the scan in Sect. 4.3, which leads to the fact that the alignment-without-decoupling conditions cannot be met anymore for perturbative values of λ and κ. This is why we consider here the usual decoupling limit, i.e. larger values of M A , in order to ensure that the h 125 properties are in agreement with the signal rate measurements at the LHC. Hence, we use a wider range of M A in our scan, extending up to M A = 560 GeV. Physical masses of the neutral BSM Higgs bosons of about 400 GeV are then achieved via mixing of the doublet states with the singlet states. Moreover, in the decoupling limit the deviations from the alignment limit scale with the factor ∆c h 125 V V = M 2 Z sin 4β/(2M 2 A ) [69] at tree level, such that also the larger value of tan β = 8 (compared to tan β = 1.7 in Sect. 4.3) already facilitates SM-like behavior of h 125 even for M A ≈ 400 GeV. 16 However, for the coupling of h 125 to down-type fermions the decoupling behavior is delayed for larger values of tan β, so that sizable deviations of c h 125 bb and c h 125 τ + τ − from the respective SM value can occur (see e.g. Ref. [69]).
We show the set of input parameters for our scan in Tab. 4. By comparing to the left plot of Fig. 11 one can see that the value for λ is far below the value required to fulfill the alignment condition given in Eq. (10). In order to avoid a sizable mixing of h 125 and h S we require that h S is much heavier than 125 GeV by using a large value of κ = 0.58 > λ. In combination with the value A κ = −200 GeV and the parameter range of µ, we find masses of 280 GeV < m h S < 442 GeV, 286 GeV < m A S < 367 GeV, 352 GeV < m H < 526 GeV and 364 GeV < m A < 522 GeV. The masses of the states with dominant doublet component (H and A) are more closely related to the value of M A , and therefore can be even larger than 500 GeV in this scan. In order to calculate χ 2 τ + τ − we sum up the cross sections of all neutral scalars in the mass interval (400±40) GeV, where interference effects can be neglected due to CP conservation and a sizable mass splitting of A S /A and h S /H in each parameter point (see also the discussion below). Thus, for the upper end of these ranges both singlet-like states can contribute to the τ + τ − excess, in addition to the doublet states H and A. Here it should be noted that the pairs A/A S and H/h S are mixed, such that the classification into doublet and singlet states is only approximate. Therefore, in the following we will adopt the mass-ordered notation h 2,3 and A 1,2 , with h 1 = h 125 .
Taking into account that we want to exploit the mixing of the singlet states with the heavy doublet states in order to obtain physical masses of m A,H ≈ 400 GeV even for considerably larger values of M A , and also that the singlet-like states can contribute to the τ + τ − excess for m h S ,A S 400 GeV, it is apparent that in the scenario considered here the excesses at 96 GeV cannot be addressed. We emphasize that attempting to fit the excesses at 96 GeV in combination with the τ + τ − excess is in any case not very promising, since (as discussed above) the description of the τ + τ − excess at about 400 GeV in the NMSSM relies on decoupling effects in order to obtain a phenomenologically viable state h 125 . Demanding a CP-even singlet-like state at 96 GeV that is mixed with h 125 in order to provide a description of the observed excesses would give rise to unacceptably large modifications of the couplings of h 125 compared to the SM. 17 The parameter ranges as defined in Tab. 4 lead to the presence of at least two neutral Higgs bosons with masses of ≈ 400 GeV. As mentioned before, for M A > 500 GeV the masses m h 3 and m A 2 can still be close to 400 GeV as long as large mixings with the lighter states are present. The parameters related to the squark sector are chosen to obtain a SM-like Higgs boson mass of ≈ 125 GeV. For a value of tan β = 8 the additional contribution to the tree-level mass of h 1 proportional to λ 2 is small. Thus, sizable radiative corrections to m h 1 are required, in analogy to the case of the MSSM. Compared to the scan presented in Sect. 4.3, we therefore increase the stop soft SUSY-breaking mass parameters to m 2 t L,R = (2.5 TeV) 2 . Furthermore, the large value of A t = 6.2 TeV yields a large stop mixing, which further increases the radiative corrections to the mass of h 1 . 18 For completeness, we list in Tab. 4 also the values of the remaining soft trilinear parameters A b,τ and the soft gaugino mass parameters M 1,2,3 . However, besides the value of M 3 = 2.7 TeV on which m h 1 depends mildly, these parameters are of no relevance in the following discussion.
We performed a simple grid scan over the parameters µ and M A . While the former is propor- tional to the singlet vev v S and thus closely related to the mass scale of the singlet states, the latter defines the approximate mass scale of the heavy doublet states. Since in this scan no Higgs boson at ≈ 96 GeV is present, we used the definition of the total χ 2 as defined in Eq. (7), i.e. without χ 2 96 , in order to quantify the quality of the fit for each parameter point. We emphasize that for the chosen value of tan β = 8 no sizable contribution to the tt excess is present. However, we keep the quantity χ 2 tt in χ 2 for completeness, as was done in Sect. 4.3 for the χ 2 τ + τ − contribution (χ 2 SM,τ + τ − = 9.99). 19 In comparison to the N2HDM analysis of the τ + τ − excess, in which only one particle gave rise to the excess, in our NMSSM analysis we have always at least two Higgs bosons with masses close to ≈ 400 GeV. For the calculation of χ 2 τ + τ − , as mentioned above, we take into account the signal contributions of all neutral BSM Higgs bosons within a mass interval of (400 ± 40) GeV, corresponding to a mass resultion of 10%. We checked for all parameter points that the mass differences of A 1 /A 2 and h 2 /h 3 are much larger than the sums of their total widths. We thus can neglect interference effects (see the discussion in Refs. [127][128][129]) and simply sum the individual contributions. They are each obtained as the product of the cross sections for bb and gg production multiplied by the branching ratio for the decay into a τ lepton pair.
We start the discussion of the scan by presenting the parameter ranges of µ and M A that pass the experimental constraints. In Fig. 16  3.5, and some points even yield a better description of the properties of h 125 than the SM (corresponding to negative values of ∆χ 2 125 ). Accordingly, the signal rates of the h 125 (= h 1 ) state agree very well with the experimental data. As discussed before, in this scan we rely on the decoupling effects in order to obtain a SM-like Higgs boson at ≈ 125 GeV. This is reflected in the fact that ∆χ 2 125 decreases with increasing values of M A . One can also see that, regarding the values of ∆χ 2 125 , there is a slight preference for lower values of µ. Comparing to the right plot of Fig. 16, in which the colors of the points indicate the value of c h 1 bb , one can see that the overall pattern of deviations in c h 1 bb is very similar to the one of ∆χ 2 125 . For the lowest values of M A for a given value of µ the coupling coefficient c h 1 bb deviates from the SM prediction by up to ≈ 6%. This leads to a small increase in the total width of h 1 , which in turn reduces the branching ratios BR(h 1 → γγ, W W * ), and thus increases ∆χ 2 125 . The deviations of c h 1 tt and c h 1 V V (not shown here) remain below 1% and are therefore much smaller than the ones of c h 1 bb . This is due to the fact that those couplings show a faster decoupling with increasing values of tan β and/or M A , while the decoupling behavior of c h 1 bb is delayed (see the discussion above).
The parameter points shown in Fig. 16 furthermore passed the HiggsBounds test. As expected, for the present scan that is targeted to a description of the τ + τ − excess we find that for all the parameter points the most sensitive constraint is the one from the searches for neutral Higgs bosons in the τ + τ − final state performed by ATLAS [7]. The precise shape of the excluded regions in Fig. 16 arises from a complicated interplay between the masses of the Higgs bosons, the mixing between h 2 /h 3 and A 1 /A 2 , and finally the number of Higgs bosons that HiggsBounds combines into a single signal that is confronted with the search limits. For instance, the lower right triangular shaped region is excluded because in this parameter region HiggsBounds combines the signal contributions of all four neutral BSM Higgs bosons into a single signal, owing to the fact that, roughly speaking, their mass differences decrease for decreasing values of the ratio M A /µ. On the contrary, for larger values of this ratio, the contributions of just the two CP-even Higgs bosons h 2,3 with masses m h 2,3 above the region of ≈ 400 GeV where the excess was observed yield a signal rate that for many points is excluded by the limits from the search for a τ + τ − resonance. This is the reason for the fact that the highest values of M A that are displayed in Fig. 16 are excluded by HiggsBounds.
In order to get an idea of the composition of the Higgs bosons that contribute to the description of the τ + τ − excess in this scenario, we show in Fig. 17 the individual values of the signal cross sections for each neutral BSM Higgs boson (i.e., the cross sections for h 125 are not displayed). The colored points fit the τ + τ − excess within the 1σ level, i.e. χ 2 τ + τ − < 2.3, whereas the gray points lie outside of the 1σ ellipse regarding χ 2 τ + τ − . In each row of Fig. 17, we also show the expected (black dashed line with green and yellow bands) and observed (black dots) exclusion limits from the ATLAS search using the full Run 2 dataset [7]. In addition, the red points indicate the expected HL-LHC exclusion limits taking into account 6 ab −1 [66]. For each parameter point there are always at least two (h 2 , A 2 ), sometimes three particles (h 2,3 , A 2 ) that contribute to the excess within the considered mass interval of (400 ± 40) GeV. For some points, even A 1 has a mass of m A 1 > 360 GeV. However, its signal contribution is small compared to the other Higgs bosons as can be seen in the top row of Fig. 17. For bb production, as shown on the left-hand side of Fig. 17, the most significant contribution is given by h 2 (third row). For the gg production mode there is no parameter point for which a single Higgs boson individually produces a sufficiently large signal, but signal rates of up to about 20 fb −1 are achieved when summing over the contributions of two or more Higgs bosons within the required mass window (see also the discussion below). Due to the value of tan β = 8 chosen for this scan, we obtain cross sections for the bb production mode that are roughly twice as large as for the gg production mode. This is well in line with the signal interpretation of the τ + τ − excess, so that one can expect a good fit to the data. The best fit point is highlighted in Fig. 17 with a magenta star. For this point, the signal contribution describing the observed excess arises from the two Higgs bosons A 2 (second row) and h 3 (fourth row) with 200 300 400 500 600 700  masses of m A 2 ≈ m h 3 ≈ 400 GeV. The remaining two BSM Higgs bosons A 1 (first row) and h 2 (third row) are roughly 100 GeV lighter and have much smaller cross sections for this parameter point and therefore do not contribute to the signal. The best fit point has χ 2 = 97.60, which is substantially below the SM value of χ 2 SM = 108.38. The difference of these values arises almost entirely from the values χ 2 τ + τ − = 0.29 and χ 2 SM,τ + τ − = 9.99, respectively. 20 It is also important to note that with the projected sensitivity the HL-LHC will be very well suited for probing the considered scenario. It will either very significantly confirm or rule out the τ + τ − excess. In the latter case, it would exclude the parameter region corresponding to the signal interpretation in our NMSSM analysis (see in particular the left plot in the second row of Fig. 17). On the other hand, in case the observed excess is indeed due to one ore more BSM Higgs boson(s) at ≈ 400 GeV the HL-LHC has excellent prospects for discovering those new states.
In Fig. 18 we show the parameter points with χ 2 ≤ χ 2 SM in the plane of the cross sections for the two production modes of the τ + τ − excess. One should keep in mind that in order to obtain χ 2 τ + τ − we took into account the contributions of all BSM Higgs bosons within the mass intervall (400 ± 40) GeV, which is of the order of the step sizes of 50 GeV of the different mass hypotheses of the experimental analysis [7]. The colors of the points indicate the values of the NMSSM parameters that were varied in this scan, where M A is displayed in the left and µ in the right plot. The best fit point is indicated by a magenta star, and the contours indicate the 1σ and 2σ ellipses of χ 2 τ + τ − [7]. One can see that the majority of points lies within the 1σ ellipse of χ 2 τ + τ − . Taking into account the scan ranges, the points in the 1σ ellipse cover the whole range of µ. However, for M A we find that only values up to M A 460 GeV yield a description of the excess at the 1σ level. For larger values the particles h 2,3 and A 2 become heavier than 440 GeV, such that only a smaller overall signal is achieved. In the left plot of Fig. 18 one can furthermore observe that there is a slight tendency towards smaller values of M A . Similarly, a tendency towards smaller µ values can be observed in the right plot. These slight preferences are however statistically insignificant in view of the current experimental uncertainties.

Prospects and outlook
As a final step of our analyses of the N2HDM and the NMSSM with respect to possible interpretations of the observed excesses, we briefly discuss the most relevant future experimental searches and measurements which can probe those scenarios. In fact, we would like to point out that there are good prospects for either strengthening the evidence for a possible signal or ruling out the discussed scenarios in the near future.
For the low tan β region, the searches for A → tt will obviously be important. The most promising channels are gluon fusion production, in which the excess was found, so far only utilizing the first year Run 2 data corresponding to 35.9 fb −1 , as well as the production of A in association with two top quarks, leading to a final state of four top quarks. The latter is already published taking into account the full Run 2 data by both ATLAS and CMS [62,63], and it is interesting to note that ATLAS measured a cross section for four-top production roughly a factor of two larger than the SM prediction [63]. In addition, in the N2HDM the parameter region in the wrong sign Yukawa coupling regime will further be probed via the process gg → A → Zh 125 , where the ATLAS analysis using the full Run 2 data set puts important constraints on the parameter space already. For the scenarios including the singlet-like Higgs boson at ≈ 96 GeV, also the searches with A decaying into a Z boson and another BSM Higgs boson can be important. Moreover, the searches for the charged Higgs bosons in the tb final state are very promising. In the NMSSM, the current limits exclude parameter points with valuess for c Att of the size of the experimental bestfit values, while in the N2HDM the charged Higgs bosons can be somewhat heavier, so that the parameter region yielding the best description of the tt excess is not affected by the current limits. For the case of the NMSSM also dedicated searches exploring decays of the charged Higgs boson into final states involving BSM particles can be promising (see also the discussion in Sect. 4.3).
For the high tan β region, i.e. tan β ≈ 8, the parameter region corresponding to the signal interpretation will fully be probed by the HL-LHC Higgs boson searches in the τ + τ − final state. If the CMS search including the full Run 2 dataset does not confirm the excess observed by ATLAS, the cross section values that are preferred by the τ + τ − excess in the corresponding ATLAS search could already be excluded on the basis of the Run 2 data from both collaborations. On the other hand, if the observed excess is indeed a first indication of a signal of one or more more BSM Higgs boson(s), the prospects for discovering these new states in future runs at the LHC and the HL-LHC will be excellent. Furthermore, the cascade decays A → Zh 125 , where the A boson is produced in the bb production mode, will also probe this scenario. It is important to note in this context that in contrast to the ATLAS search in the gg production mode, the ATLAS search assuming bb production has not yet been updated to include the full Run 2 data set.
As already mentioned before, independently of the value of tan β the presence of relatively light charged Higgs bosons is a common prediction of the scenarios discussed in this paper, where in the NMSSM the charged Higgs bosons are even lighter than the CP-odd Higgs boson at 400 GeV, while in the N2HDM an upper bound of m H ± 750 GeV applies because of the theoretical constraints. In the NMSSM the constraints from flavor physics do not result in a firm lower limit on m H ± , since in general the theoretical predictions for the flavor observables in SUSY models depend on various different sectors of the model, and may be weakened without changing the Higgs-boson phenomenology discussed here. In the N2HDM, on the other hand, we included a lower limit of m H ± = 550 GeV, based on flavor physics constraints [59]. More recently, a new theoretical calculation suggested a lower limit of m H ± > 800 GeV [130], which however is still under debate in particular in view of the results of Ref. [131] pointing to possible underestimates of theoretical uncertainties that could have the effect of even weakening the bound on m H ± below the value that we have adopted in our scan. In any case it is obvious that new results concerning flavor physics observables can have an important impact on the investigated parameter space.
Regarding the compatibility of the predicted signal rates of h 125 with the experimentally measured values, our results indicate that projected HL-LHC measurements of the properties of h 125 alone might not be sufficient to exclude (or confirm) the scenarios describing the excesses at 400 GeV. The situation is somewhat more promising when additionally the presence of a singletlike Higgs boson at 96 GeV is considered. In this case, the mixing of the singlet-like Higgs boson with h 125 , as required to simultaneously accommodate the LEP and CMS excesses, gives rise to modifications of the signal rates of h 125 compared to the SM predictions that could be observable at the HL-LHC.
The NMSSM is additionally constrained by experimental searches targeted specifically to the SUSY particles. For the scenario in the alignment-without-decoupling limit, it was found that the stops cannot be much heavier than 1 TeV in order to obtain a theoretical prediction for the mass of h 125 in agreement with experiment. The HL-LHC has a high sensitivity for probing this mass region. Another possibility in this context are searches for the light electroweakinos, where in particular searches focusing on compressed electroweakino spectra are promising. More exotic signals that are present in our analysis arise from the decays of the charged Higgs bosons into pairs of a neutralino and a chargino, which as mentioned above could be probed via dedicated searches. Finally, a small fraction of the analyzed parameter points feature a neutralino with masses smaller than 62.5 GeV, such that the decay of h 125 into a pair of neutralinos is kinematically allowed. Such points will be further probed by direct searches for invisible decays of h 125 and by the global constraints from the signal-rate measurements of h 125 at the HL-LHC.

Conclusions
Various searches at the LHC for BSM Higgs bosons with masses above 125 GeV showed excesses over the background expectation in the Run 2 data. As a remarkable feature, several of these excesses were found around a mass scale m φ ≈ 400 GeV of a hypothetical new Higgs boson φ (or more than one). In particular, we focused on possible interpretations of the two most striking excesses: CMS reported a local excess of more than 3 σ in the channel A → tt in their first year Run 2 data [6], while ATLAS reported a local excess of about 3 σ in the channel φ → τ + τ − in their full Run 2 data [7]. In both cases the analysis of the other experiment for the same period of data taking is not yet available. In addition, a local excess of more than 3 σ was found in the ATLAS search for a heavy resonance decaying into a Z boson and h 125 assuming bb-associated production and utilizing 35 fb −1 of data [8]. While the various excesses reach the level of 3 σ local significance, all stay below 3 σ global significance. Also the searches at the LHC for BSM Higgs bosons with masses below 125 GeV show an interesting excess at about 96 GeV in the channel pp → φ → γγ: CMS reported a local excess of about 3 σ in their first year Run 2 data [17] and a similar deviation of 2 σ local significance in their Run 1 data at a comparable mass [20]. The ATLAS results based on the data of the first two years of Run 2 [19] are not sensitive to the excess observed at CMS. Furthermore, LEP reported a local excess of about 2 σ in the channel e + e − → Z φ → Z bb [14].
In this paper we have analyzed whether the observed excesses can be described by models comprising an extended Higgs sector, where we have concentrated on the Next-to-2HDM (N2HDM) and the Next-to-MSSM (NMSSM) as minimal scenarios for accommodating the experimentally observed patterns. Concerning the excesses at about 96 GeV from CMS and LEP taken in isolation, it is known that several models can fit both excesses simultaneously at the level of roughly 1σ or less, including the N2HDM and the NMSSM. In the present paper we have analyzed the question whether the excesses at about 400 GeV can be accommmodated in these models and whether the parameter space that is preferred by those analyses in addition permits a signal interpretation also of the excesses at about 96 GeV. In our analysis we have taken into account the experimental constraints on the properties of the Higgs boson at 125 GeV, the limits from BSM Higgs-boson searches at the LHC and previous colliders, constraints from electroweak precision data and flavor observables, as well as theoretical constraints such as vacuum (meta)stability and perturbativity.
We first investigated whether the N2HDM can fit the two excesses at 400 GeV such that the relevant parameter region is in agreement with the theoretical and experimental constraints. While we find that both excesses cannot be fitted simultaneously, each of them can be described with a very good χ 2 by the type II N2HDM, while complying with the other constraints. The tt excess, independently of the N2HDM type, can be described for 1 tan β 2, whereas a description of the τ + τ − excess in the N2HDM type II requires 6 tan β 11, making them mutually incompatible. In both cases the description of the excesses mainly occurs as a consequence of the presence of a CP-odd Higgs boson with m A ≈ 400 GeV, while the other Higgs bosons are found to be not much heavier in those scenarios owing to theoretical constraints. Interestingly, even though not directly aimed for in our analysis, a sizable fraction of the parameter points that fit the τ + τ − excess additionally yield a contribution that would be compatible with a description of the A → Zh 125 excess that was found at a similar mass.
In a second step, we have demonstrated that for each of the parameter regions that are preferred by the excesses at about 400 GeV the N2HDM can simultaneously also describe both excesses at about 96 GeV. This happens through the presence of a singlet-like scalar at 96 GeV giving rise to signal rates that are compatible with the excesses observed at CMS and at LEP. More specifically, in type II either the tt or the τ + τ − excess at about 400 GeV can be described in their respective tan β range, where in the latter case also a sizable contribution to the Zh excess can be present. In type IV the tt excess is compatible with the 96 GeV excesses, where the CMS excess can be described at the level of roughly 1σ.
In the NMSSM, which has a type II Yukawa sector, as a consequence of the underlying symmetry relations of the model one has much less freedom regarding the choice of parameters in the Higgs sector as compared to the N2HDM. Nevertheless, we have demonstrated that also the NMSSM can fit each of the excesses at 400 GeV individually, while complying with the BSM Higgs-boson searches and Higgs-boson measurements as well as the other constraints. Concretely, we discussed exemplary scenarios in which we scanned over M A , µ, and in case of low tan β also over κ and A κ . The tt excess is described at the level of 1.5 σ for low tan β in the alignmentwithout-decoupling limit, which is a theoretically and, in view of the constraints from the Higgs boson measurements, experimentally well motivated scenario. The τ + τ − excess, on the other hand, can be described for larger tan β values somewhat away from the alignment limit, while the properties of the Higgs boson at about 125 GeV are nevertheless compatible with the experimental data within the current uncertainties.
In the NMSSM parameter region which can give rise to the tt excess naturally a light singletlike scalar is present. Requiring its mass to be ≈ 96 GeV, we demonstrated that the CMS excess at about 96 GeV can be described simultaneously with the tt excess at about 400 GeV, whereas in this scenario hardly any signal contribution would have been generated at LEP.
If the observed excesses indeed turn out to be first indications of one or more BSM Higgs bosons, besides the channels featuring those excesses particularly promising channels regarding a possible discovery of a BSM Higgs boson would be the search for charged Higgs bosons with masses around 400 GeV and searches for the heavier doublet-like Higgs bosons decaying to a Z boson and either the SM-like Higgs bosons or the CP even/odd singlet-like scalar with masses below or slightly above 125 GeV. The scenarios describing the τ + τ − excess will entirely be probed by the φ → τ + τ − searches at the HL-LHC. Thus, there are very good prospects in the near future for clarifying the tantalizing hints that have been observed in the Higgs searches via new results from searches, more precise measurements of the properties of the Higgs boson at 125 GeV and further information, for instance from flavor physics.
Here c Aff are the effective coupling coefficients of A to the SM fermions, which are defined as the coupling of A relative to the one of a hypothetical SM Higgs boson of the same mass, see Sect. 2. In the N2HDM, the coefficients are given depending on the type either by tan β or 1/ tan β, as depicted in Tab. 5. Taking into account that for all four types |c Att | = 1/ tan β holds, one can conclude that regarding the tt excess all four types are equivalent. A sufficiently large coupling |c Att | ≈ 1 is obtained for values of tan β not much larger than one.
The situation is different regarding the τ + τ − excess. Here, also the couplings of A to bottom quarks and, obviously, τ leptons are important. In the following, we will separately discuss for each Yukawa type whether it could realize the τ + τ − excess. We will also distinguish between the low tan β region, tan β ≈ 1, in which a simultaneous explanation of the tt excess would be possible, and the large tan β region, tan β ≈ 10.
In the N2HDM type I, the τ + τ − excess cannot be realized for any value of tan β. The reason is twofold: Firstly, accommodating the observed pattern of the excess requires that the bb associated production cross section σ(bb → A) is roughly of the same size or larger than the gg production cross section σ(gg → A). The latter is mainly mediated via the diagram with the top quark in the loop. Due to the fact that the Yukawa coupling Y t is an order of magnitude larger than Y b , the condition σ(bb → A) σ(gg → A) requires that |c Abb | |c Att |, see above. However, in type I c Abb = c Att holds. In addition, BR(A → τ + τ − ) is always tiny in type I, because c Att = c Aτ + τ − , and for m A = 400 GeV the decay A → tt is kinematically open. Moreover, it was shown in Ref. [21] that the excesses at 96 GeV cannot be realized in type I. Thus, the only excess that can be accommodated in type I is the tt excess.
In type II the situation is different. Here one has |c Abb /c Att | = tan 2 β. Thus, the suppression of Y t can be compensated for tan β ≈ 10, so that a possible excess can occur also for the bb production mode. Moreover, |c Aτ + τ − /c Att | = tan 2 β holds, and therefore BR(A → τ + τ − ) can be sizable due to the enhancement of c Aτ + τ − while BR(A → tt) is suppressed. The question whether the τ + τ − excess can be correctly reproduced depends on whether there is a range of tan β in which both the enhancement of σ(bb → A) and the one of BR(A → τ + τ − ) are of a size such that the rates in both production modes are compatible with the observed patterns. This issue will be addressed in our scan together with the question whether both the tt and the τ + τ − excess can be accommodated simultaneously. Since it is possible in type II to realize the excesses at about 96 GeV [21], we will subsequently perform also a dedicated scan with m A = 400 GeV and m h 1 ≈ 96 GeV for this type.
In type III the τ + τ − excess cannot be accommodated for any value of tan β, because c Att = c Abb , so that the requirement σ(bb → A) σ(gg → A) cannot be realized (as in type I  excesses at 96 GeV cannot be accommodated in type III [21]. Consequently, only the tt excess can be explained in type III. Also in type IV the τ + τ − excess cannot be realized. The coupling coefficient c Aτ + τ − scales with the inverse of tan β, such that tan β < 1 would be required to enhance BR(A → τ + τ − ). On the contrary, in analogy to the argument given for type II, the condition σ(bb → A) σ(gg → A) requires values of tan β ≈ 10. Thus, depending on the value of tan β, either σ(bb → A) or BR(A → τ + τ − ) is too small. Since type IV is capable of explaining the excesses at 96 GeV [21], there is the possibility to accommodate both excesses at 96 GeV together with the tt excess at 400 GeV in type IV.
As a consequence of the above considerations we restrict our parameter scans investigating exclusively the excesses at 400 GeV in the N2HDM to type II. For the other three types the τ + τ − excess cannot be fitted, and the discussion of the tt excess would be qualitatively very similar to the case of type II (the only difference that could arise is that the total width of A can be different in type I, III and IV compared to type II, leading to slightly different preferred tan β values due to different best-fit values of c Att , see Sect. 2). The results of our scan for type II of the N2HDM are described in Sect. 3.3. Concerning the question whether there is the possibility of fitting both the tt excess and the τ + τ − excess at 400 GeV simultaneously, based on our qualitative discussion above we can already anticipate that there will be a tension between the tan β values required to fit either of them. Indeed, we find that the excesses point towards different regions of parameter space with either tan β 3 for the tt excess or tan β 6 for the τ + τ − excess.
In a second step, we will analyze whether the excesses at about 96 GeV can be realized in combination with either of the excesses at about 400 GeV. Thus, we perform two different scans with tan β low = [0.5, 4] (regarding the tt excess at 400 GeV) and tan β high = [6, 12] (regarding the τ + τ − excess at 400 GeV). The results are detailed in Sect. 3.4 and Sect. 3.5. As discussed above, in type IV one can potentially fit the excesses at 96 GeV in combination with the tt excess for tan β = tan β low . A corresponding parameter scan is described in App. A.2.

A.2 Higgs bosons at 96 GeV and 400 GeV for low tan β in type IV
We present here the results for a scan in the low tan β regime in the type IV (flipped) N2HDM. The input parameters were chosen to be identical to the scan in type II as shown in Tab. 2. We show the results of the scan in type IV in the tan β-c Att plane in Fig. 19 and in the µ LEP -µ CMS plane in Fig. 20.
Comparing Fig. 19 with the corresponding type II results shown in Fig. 6, one can see that the values of c Att found in the scan are considerably smaller in type IV. This is related to the smaller predicted widths of A for the parameter points that fulfill χ 2 ≤ χ 2 SM in type IV compared to type II, as indicated in the left plot of Fig. 19. The experimental result shows a less pronounced excess for Γ A /m A ≈ 1.5% . . . 2.0% [6], and the minimum of χ 2 tt is shifted towards smaller couplings (see also Fig. 1), such that points with low values of χ 2 tt require smaller couplings. The smaller values of the width Γ A in type IV compared to type II have their origin in the modified leptonic couplings. For values of tan β > 1 the decay width for the decay A → τ + τ − is suppressed in type IV, while it is enhanced in type II (see also Tab. 5). In the right plot of Fig. 19, in which the color coding indicates the value of m H ± , one can observe that m H ± is substantially smaller than the upper limit of the scan range 1 TeV for all the parameter points. The same result was found already in Sect. 3.4 for type II. The reason for the preference for relatively small values of m H ± lies again in the theoretical constraints on the absolute values of the quartic couplings (where m A is kept fixed at 400 GeV), which are independent of the Yukawa type. Comparing Fig. 20 to the results of type II (see Fig. 7), one can see that the points are substantially shifted towards lower values of µ CMS . As explained in Ref. [21], in type IV BR(h 1 → τ + τ − ) is enhanced in the parameter space in which the CMS excess can be accommodated, giving rise to smaller values of BR(h 1 → γγ). While it is not possible to reach χ 2 96 values close to zero, corresponding to the center of the displayed ellipse, several points lie within the 1σ ellipse of χ 2 96 while in addition accommodating the tt excess. This is indicated by the low values of χ 2 tt shown in the left plot of Fig. 20. The lowest values of χ 2 tt are at the level of ≈ 1.1. This is slightly larger than the lowest values found in type II (≈ 0.25), as shown in Fig. 7. This difference has its origin in the smaller values of the width of A in type IV, which yield a worse fit to the data, as can be seen in Fig. 1. Taking the issues mentioned above into account, we conclude that there is a slight tension between fitting the CMS excess at around 96 GeV and the tt excess at around 400 GeV in type IV, but a realization of all three excess is possible at the 1σ level of χ 2 96 . Comparing the right plot of Fig. 20 to the one in Fig. 7, one can see that the lowest values of ∆χ 2 125 found in the scan for type IV are larger than the ones for type II. This suggests that larger deviations of the properties of the Higgs boson at ≈ 125 GeV compared to the SM prediction can be expected in type IV. The higher values for χ 2 125 in combination with a worse fit to the CMS excess also manifest themselves in the χ 2 value of the best fit point χ 2 = 102.3, which is substantially larger than the corresponding value χ 2 = 97.93 found in the type II analysis. Consequently, since in both the type II and the type IV analysis the condition χ 2 ≤ χ 2 SM was required to be fulfilled for each point, the total number of points within and near the 1σ ellipse for χ 2 96 in Fig. 20 compared to Fig. 7 is substantially smaller. In summary, our analysis in this section has revealed that for the low tan β region of the N2HDM type II provides a better description of the observed data than type IV. 125 rather close to zero, as indicated by the color coding of the points, and the properties of the particle state at 125 GeV are in good agreement with the experimental constraints given the current uncertainties. We find ∆χ 2 125 ≈ 2 as lowest values, see the discussion below. One can furthermore see two small clusters of points at c h 125 V V ≈ 0.994 for which the values of ∆χ 2 125 are considerably larger compared to the surrounding points. These clusters arise from additional decay channels of the state h 125 , as will be discussed below, and they appear distinctively in Fig. 21 because they correspond to the two isolated points (lines of points) at the lower end of κ in the left (right) plot of Fig. 12. As mentioned before, a point is considered to be allowed when the condition χ 2 ≤ χ 2 SM is fulfilled. Thus, for points with ∆χ 2 125 considerably above zero, it is expected that they provide predictions describing at least one of the excesses, such that the penalty of χ 2 125 within the total χ 2 is compensated by either χ 2 tt < χ 2 SM,tt , χ 2 96 < χ 2 SM,96 , or both. In order to further elucidate our results for ∆χ 2 125 and to determine the origin of the larger values of ∆χ 2 125 in the two clusters mentioned above, we show in Fig. 22 a selection of h 125 branching ratios. On the left-hand side, one can see that both BR(h 125 → W W * ) and BR(h 125 → γγ) agree with the SM prediction (indicated by the blue dashed lines) at the level of ≈ 5-10%. For all the points the value of BR(h 125 → W W * ) is slightly below the SM prediction. This has two reasons: firstly, the coupling coefficient c h 125 V V ranges between 0.983 and 0.997, as can be seen in the left plot of Fig. 21, such that a corresponding suppression is also expected for the partial widths of the decays to W W * and ZZ * . Secondly, as shown in the right plot of Fig. 21, the coefficients c h 125 bb are a few percent above the SM prediction for all of the points. This leads to a small enhancement of the partial decay width to b quarks, and therefore also to an enhancement of the total width of h 125 , which in turn yields a further suppression of BR(h 125 → W W * ). Overall, the small deviations of the branching ratios for h 125 → V V * , γγ compared to the SM predictions give rise to a penalty of ∆χ 2 125 ≈ 2-6. The two isolated clusters of points mentioned above are also well visible in the left plot of Fig. 22. They have an even slightly smaller BR for the decay to W W * than the rest of the points. This can be understood from the right plot of Fig. 22, where we show the plane BR(h 125 → χ 0 1 χ 0 1 )-BR(h 125 → τ + τ − ) with the color coding indicating again ∆χ 2 125 . The blue dotted lines indicate the SM prediction. Most of the points are found at BR(h 125 → χ 0 1 χ 0 1 ) = 0. However, the two isolated clusters of points yield a neutralino with a mass m χ 0 1 62 GeV, such that the decay h 125 → χ 0 1 χ 0 1 becomes kinematically allowed. The corresponding partial decay width contributes to the total decay width of h 125 and reduces the branching ratios for decays into SM particles. While for h 125 → W W * this results in a larger deviation from the SM prediction, for the decays h 125 → γγ, τ + τ − this brings the prediction closer to the SM value. Overall, depending on the size of BR(h 125 → χ 0 1 χ 0 1 ), this yields values of ∆χ 2 125 ≈ 5-8. In the right plot of Fig. 22 we have also indicated the current upper limit on BR(h 125 → inv.) as reported by ATLAS [132]. For our analysis this bound does not yield a restriction since the global experimental constraints on the signal rates of h 125 , as tested by HiggsSignals, already exclude all points potentially featuring values of BR(h 125 → χ 0 1 χ 0 1 ) 0.07. [5] ATLAS collaboration, Combined measurements of Higgs boson production and decay using up to 80