Applying exclusion likelihoods from LHC searches to extended Higgs sectors

LHC searches for non-standard Higgs bosons decaying into tau lepton pairs constitute a sensitive experimental probe for physics beyond the Standard Model (BSM), such as supersymmetry (SUSY). Recently, the limits obtained from these searches have been presented by the CMS collaboration in a nearly model-independent fashion – as a narrow resonance model – based on the full 8TeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8\,\, \mathrm {TeV}$$\end{document} dataset. In addition to publishing a 95%C.L.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$95~\%~\mathrm {C.L.}$$\end{document} exclusion limit, the full likelihood information for the narrow resonance model has been released. This provides valuable information that can be incorporated into global BSM fits. We present a simple algorithm that maps an arbitrary model with multiple neutral Higgs bosons onto the narrow resonance model and derives the corresponding value for the exclusion likelihood from the CMS search. This procedure has been implemented into the public computer code HiggsBounds (version 4.2.0 and higher). We validate our implementation by cross-checking against the official CMS exclusion contours in three Higgs benchmark scenarios in the Minimal Supersymmetric Standard Model (MSSM), and find very good agreement. Going beyond validation, we discuss the combined constraints of the ττ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau \tau $$\end{document} search and the rate measurements of the SM-like Higgs at 125GeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$125\,\, \mathrm {GeV}$$\end{document} in a recently proposed MSSM benchmark scenario, where the lightest Higgs boson obtains SM-like couplings independently of the decoupling of the heavier Higgs states. Technical details for how to access the likelihood information within HiggsBounds are given in the appendix. The program is available at http://higgsbounds.hepforge.org.


Introduction
The search for Higgs bosons [1][2][3][4][5][6] continues to be a cornerstone of the physics program at the Large Hadron Collider (LHC). After the discovery of a Higgs boson by ATLAS [7] and CMS [8] it is crucial to find out whether the detected particle is part of a Higgs sector that contains several physical states. Higgs sectors of this kind are predicted in many theories of physics beyond the Standard Model (SM). For the understanding of the mechanism of electroweak symmetry breaking two complementary experimental endeavors are important: On the one hand the precise determination of the properties of the Higgs signal detected at around 125 GeV, and on the other hand the search for additional Higgs bosons. Both are crucial in the quest to identify the underlying physics. The existing limits from the Higgs searches at LEP, the Tevatron and the LHC already put very important constraints on the parameter spaces of different models that provide a Higgs-like state compatible with the detected signal. More data on both the detected signal and on searches for additional Higgs bosons will further enhance the sensitivity for discriminating possible scenarios of new physics from the SM and from each other.
In order to facilitate the available experimental information from the Higgs searches at LEP, the Tevatron and the LHC, expressed in terms of relatively model-independent cross-section limits for testing a wide variety of theoretical models, the program HiggsBounds [9][10][11][12] has been developed. Experimental information on the Higgs signal detected at a mass value of around 125 GeV is utilized in the sister program HiggsSignals [13] for testing the theoret-ical predictions from any kind of Higgs sector. The experimental information on the detected signal incorporated in HiggsSignals is turned into a χ 2 likelihood, which is suitable for the inclusion into global fits (see, e.g., Refs. [14][15][16][17][18][19]), where in addition many other observables are taken into account. In contrast, exclusion limits have traditionally been presented in terms of 95 % C.L. limits, which a priori only provide the information whether a particular parameter point is excluded or not at the 95 % C.L. by the considered search channel. In a global fit, where the predictions of a model are confronted with a large number of observables, it would usually be too restrictive to disregard a certain parameter point just because it falls outside of the 95 % C.L. region of a single search channel. In fact, testing a large variety of observables one would expect that the measured values of some observables lie outside of the respective 95 % C.L. regions for purely statistical reasons. It would therefore be very desirable if also negative experimental outcomes from Higgs searches were provided in terms of the likelihood information in the relevant parameters, instead of a simple binary rejection or acceptance at a certain confidence level (C.L.). Up to now, likelihood information was available in HiggsBounds only for the results from the LEP Higgs searches [12], while for all search channels at the Tevatron and LHC only 95 % C.L. limits have been accessible. We report here on significant progress in this direction for the LHC Higgs boson search in the τ + τ − final state, which plays a central role in the search for additional Higgs bosons. Many models that can accommodate a SM-like Higgs boson at 125 GeV, such as the Minimal Supersymmetric Standard Model (MSSM) or the various types of Two-Higgs-Doublet Models (2HDM), predict additional Higgs bosons that decay predominantly into SM fermions. Therefore, LHC searches for new neutral Higgs bosons decaying to τ + τ − play a crucial role. In particular within the MSSM, these searches lead to large excluded regions in the parameter space. The highest experimental sensitivity occurs for smaller values of the CP-odd Higgs boson mass, M A , and larger values of tan β, the ratio of the two vacuum expectation values [20,21].
One complication that arises for this search channel is the fact that two different production modes, gluon fusion and b quark associated production, can both be important. Their individual contributions to the signal rate can vary strongly over the parameter space. Since the acceptances of these two channels can also be very different, a two-dimensional cross section interpretation for the τ + τ − final state is desirable as a basis for (close to) model-independent exclusion limits. Recently, the CMS collaboration has published the likelihood information for their Higgs boson search in the τ + τ − final state [20]. The likelihood is given as a function of the two relevant Higgs production channels, gluon fusion and b quark associated production, for various mass values of the narrow resonance assumed for the signal model.
In this paper we investigate the application of this new experimental information for testing the theoretical predictions of extended Higgs sectors and its incorporation in global fits. We develop a simple algorithm that maps an arbitrary model with in general several neutral Higgs bosons onto the narrow resonance model. In this way the corresponding value of the exclusion likelihood from the CMS search for the tested model can be determined. We furthermore describe the inclusion of this likelihood information into the publicly available Fortran code HiggsBounds [9][10][11][12]. For nearly any model under consideration, HiggsBounds provides an evaluation of the exclusion likelihood for a model parameter point based on the information from [20]. While the new likelihood information goes well beyond the standard test whether a particular parameter point is excluded at the 95 % C.L., the likelihood information can also be employed to run Higgs-Bounds in this "standard" mode. In this case, Higgs-Bounds determines the parameter region that is excluded at the 95 % C.L. based on all available searches, including the new τ + τ − result from CMS. The new version of HiggsBounds can be used together with its sister program HiggsSignals [13] in order to take into account both the information from search limits and from the detected signal for a comprehensive test of Higgs phenomenology. Both codes are available at: http://higgsbounds.hepforge.org.
The paper is organized as follows. In Sect. 2 we briefly summarize the experimental results that are used as input for our investigations. Details of the employed algorithm and the implementation of the exclusion likelihood of Ref. [20] into HiggsBounds are given in Sect. 3. The validation in various MSSM Higgs benchmark scenarios is discussed in Sect. 4. As an application, in Sect. 5 we investigate the constraints on a certain benchmark scenario in the MSSM that are obtained from using the new exclusion likelihood in combination with the information on the detected signal incorporated in HiggsSignals. We conclude in Sect. 6. Finally, all relevant information needed to run HiggsBounds to obtain the likelihood information for the τ + τ − Higgs search channel for any parameter point under investigation are contained in an Appendix.

Experimental results
This section briefly summarizes the experimental results from the CMS non-standard Higgs search in the τ τ final state [20], that we have used as starting point for our investigation and that we have implemented in HiggsBounds. The search analysis of CMS is carried out in two separate selection categories: One requiring the presence of at least one b-tagged jet, and one without the presence of a b-tag. The former category is enriched by the production of a Higgs boson, denoted generically by φ, in association with two b quarks, gg → bbφ, while the latter is dominated by the gluon fusion process, gg → φ. Hence, the search features sensitivity to the two different production modes separately, which enables the presentation of the search results in terms of individual signal strengths in both production modes for all tested Higgs boson masses. Separate information on the two production modes is an indispensable ingredient for enabling the presentation of (close to) model-independent exclusion limits or measurements, in case of a discovery.
The data is further classified in categories defined by the two τ lepton decay modes: eτ h , μτ h , eμ, μμ and τ h τ h , where τ h denotes a hadronically decaying τ lepton. Using a maximum likelihood technique, an estimator for the true τ τ invariant mass, m τ τ , is reconstructed from the momenta of the visible τ decay products and the missing transverse energy in the event. The uncertainty of the m τ τ reconstruction is estimated to be around 20 % when averaged over all decay modes [20].
The resulting m τ τ spectrum in all categories (b-tag and τ decay) separately is then subject to a profile likelihood analysis [22], where the background parametrization, obtained from control region data and Monte Carlo simulation, and the signal shape parametrization are fitted simultaneously to the reconstructed mass spectrum. The fit is performed individually for test masses m φ between 90 GeV and 1 TeV, and the results are interpolated between the test masses.
CMS interprets the results both in a nearly modelindependent way 1 for a single narrow resonance φ, and, in a model-specific context, in the MSSM, where three neutral Higgs bosons h, H and A potentially comprise the signal. The latter interpretation performs a likelihood ratio hypothesis test for the two hypotheses of a single SM-like Higgs boson at m h = 125 GeV with exact SM properties and, alternatively, for the signal consisting of all three neutral Higgs bosons of the MSSM. In the latter case, the m τ τ distributions of the h/H/A → τ τ decays are combined before the calculation of the likelihood. Note that, whereas the modelspecific limits for the MSSM are based on the full integrated luminosity of the combined 7 + 8 TeV dataset, the results for the single narrow resonance model are obtained from only the 8 TeV dataset.
Since HiggsBounds is designed to test any extended Higgs sector, with any coupling properties of the 125 GeV Higgs candidate (if the model under consideration provides such a candidate; obviously the phenomenological interest in other models is rather limited) and any masses and properties of the remaining Higgs spectrum, the nearly modelindependent single resonance results are chosen for the implementation in HiggsBounds. On the one hand, this is the only possibility unless one is willing to adopt further model-dependent assumptions. On the other hand, it will in general yield weaker, i.e. more conservative, limits than a dedicated analysis taking into account the full structure of the considered model. For example, in the MSSM this will yield a conservative limit whenever either m H/A ≈ m h , or, more generally, whenever the model predicts that more than one Higgs boson have a non-negligible signal yield and contribute in different regions in m τ τ . Since the likelihood is constructed for single resonances, such a case cannot be properly reconstructed from the likelihood. However, if e.g. the heavy MSSM Higgs bosons H and A contribute at the same point in m τ τ , their signal rates can be added before interpreting the likelihood. In this case the implementation is not necessarily conservative. A detailed study on the applicability of these limits to the MSSM benchmarks is presented in Sect. 4.
The profile likelihood analysis follows the standard implementation: The test statistics is given by Here N is the number of measured events, b and s(m) the number of expected background and signal events for a given resonance mass hypothesis m, μ the signal strength modifier, and θ are the nuisance parameters describing the systematic uncertainties.θ μ maximizes the likelihood in the numerator given a certain value of μ, whereas the likelihood reaches its global maximum atμ andθ , which is given in the denominator. The constraint 0 ≤μ ≤ μ is employed to not penalize the model for a possible excess of the data over the signal plus background prediction. It should be noted that the signal yield contains two independent components, corresponding to the two production modes gg → φ and gg → bbφ. Thus, μ and s(m) are two-component vectors.
Using toy Monte Carlo techniques or asymptotic expressions for large statistics, the expected probability distributions P(q μ |hypothesis) can be constructed for the test statistics given above. In the model-independent analysis used here, the hypothesis either consists of H 1 = μ · s(m) + b for the case of the presence of a single narrow resonance with a given mass m and signal yield μ · s(m), and of H 2 = b for the background. Using these hypothesis definitions and the observed value of the likelihood ratio, q obs μ , given by Eq. (1) with N given by the actual observed number of events, N = N obs , the likelihood ratio technique can be used to define the CL s as In a stand-alone search, the criterion CL s ≤ α is then used to exclude the presence of a signal at 1 − α confi-dence level (C.L.). For the expected limit, the observed data in the calculation of q obs μ is replaced by the median of the background-only expectation for q μ . A model-independent limit, e.g. at the 95 % C.L., can then be derived by varying μ until CL s = 0.05. The value μ at which this happens then represents the signal strength modifier which is just allowed at the 95 % C.L. For any given model, HiggsBounds reconstructs the predicted signal yield s(m) from the theoretical input provided by the user, and obtains the corresponding value of the test statistics q μ=1 (or simply denoted q model ) from the CMS likelihood data. The details of this procedure will be described in the following section. Note that HiggsBounds directly employs the expectation and observation of the test statistics, q exp μ and q obs μ , respectively, as the provided CMS data does not allow for a full reconstruction of the CL s value. Nevertheless, in the limit of large numbers, the test statistics can be approximated by a chi-squared differences function above minimum, q μ ≈ χ 2 = χ 2 − χ 2 min , as Thus, for example, we approximately obtain the twodimensional limit at the 68 and 95 % C.L. when the test statistics q μ takes the values 2.28 and 5.99, respectively.
As an example, we show the observed likelihood distribution, q obs μ , for test masses of 125 and 300 GeV in Fig. 1. We also indicate the approximated 68 and 95 % C.L. limits by contour lines. It should be kept in mind that these contours are only for illustrational purposes. In HiggsBounds the full likelihood information q μ is used, and the specific limit at a certain C.L. can be easily obtained from this information.
The implementation of the likelihood for the CMS φ → τ τ search differs in two significant ways from the implementation of the LEP Higgs search χ 2 in HiggsBounds, which is already available since version 2.0.0 [10]. In the LEP implementation, each Higgs search channel, comprised of one production mode and one Higgs boson decay mode, is treated separately, thus no combination of production modes is applied or possible for the user. In addition, in the LEP implementation the χ 2 is estimated from the CL s+b value in each channel at the given signal strength prediction using Gaussian approximations. In contrast, for the CMS φ → τ τ search the exact values of the test statistics q μ as presented by CMS are used and properly combined for both production modes.

Likelihood reconstruction for extended Higgs sectors
For the construction of the exclusion likelihood from the H → τ τ search, we make use of the following quantities: For each neutral Higgs boson, h i (i = 1, . . . , N ), in a model with N neutral Higgs bosons, we have a prediction of the mass, m i (where the relevant range is currently m i ∈ [90, 1000] GeV), the gluon fusion production cross section, σ (gg → h i ), the cross section for production in association with b quarks, σ (gg → bbh i ), and the branch- The main algorithm for the likelihood reconstruction proceeds as follows: 1. Signal rates of multiple Higgs bosons that cannot be resolved by the experimental analysis are added. We thus combine the signal predictions for two Higgs boson h i and h j ( j = i), if Each Higgs boson can appear in different such combinations. For each combination, also called Higgs cluster and labeled with capital characters in the following, we evaluate the physical quantities as follows: We assume that the total rates are given by the incoherent sum of the signal rates of the individual Higgs bosons in the cluster, The cluster mass, m I , is determined by a signal strengths weighted mass average The sums in Eqs. (6) [23,24]. An extension of HiggsBounds that enables the implementation of interference effects of nearby resonances in a generalized narrow-width approximation is currently under development, see also Ref. [25]. 2. In the second step, the expected and observed likelihood values, q exp model and q obs model , respectively, for each Higgs cluster h I are evaluated from the experimental likelihood data grid. The likelihood is first evaluated for the rate values σ (gg → h I → τ τ ) and σ (gg → bbh I → τ τ ), obtained through Eqs. (6), (7), in the mass-neighboring data slices, i.e. at the nearest grid mass values below and above m I , denoted by m − and m + , respectively. The likelihood value at the predicted cluster mass m I is then obtained through linear interpolation: Here q +/− denote the values of the test statistics obtained at the neighboring grid above or below the predicted mass Following this algorithm, the full likelihood from the CMS φ → τ τ analysis for both the expected and observed exclusion can be directly obtained within HiggsBounds for any tested model. This is carried out via Fortran subroutines. For a technical documentation see Appendix.
The use of this likelihood information is complementary to the other type of information contained in a full Higgs-Bounds application, which considers exclusion limits from many other Higgs searches from the LEP, Tevatron and LHC experiments. As an alternative to using the full likelihood, we therefore also provide the option to reconstruct a limit at 95 % C.L. and use this in the "standard" HiggsBounds operation. For clarity, we now repeat some elements of how this works [12]: In the statistical procedure, HiggsBounds first determines the most sensitive analysis to the model by picking the analysis application, for which the ratio between the model-predicted signal rate, S predicted , over the expected upper limit on the signal rate, S 95 % CL expected , is maximized. After the most sensitive analysis has been determined, the model prediction is confronted with the observed exclusion limit of this particular analysis, S 95 % CL observed . The model is considered to be excluded at the 95 % C.L., if For the CMS φ → τ τ analysis described above, S 95 % CL expected and S 95 % CL observed are a priori not known and need to be determined from the implemented likelihood distribution. In a numerical procedure, we therefore scale the model-predicted gg → φ → τ τ and gg → bbφ → τ τ rates with a universal factor μ until the obtained expected/observed likelihood q exp/obs μ values are equal to 5.99, corresponding to the two-dimensional 95% C.L. interval. 2 The so-obtained scale factors, μ exp/obs 95 % CL , are then identified with the expected/observed 95 % C.L. upper limits on the signal rate, respectively, which enter Eqs. (10) and (11). In this way, the likelihood-based results from the CMS φ → τ τ analysis can be incorporated in the standard HiggsBounds run.

Validation
Besides the nearly model-independent limits, CMS has also presented model-specific interpretations of their search results. This has been done for the MSSM, employing the benchmark scenarios proposed in Ref. [26] (see also Ref. [27]). Here, we validate our likelihood implementation in HiggsBounds against the CMS results for three of these scenarios: The m max h , the light stop and the low-M H scenarios (see Ref. [26] for details). The comparison of the reconstructed 95 % C.L. exclusion line with the official CMS result provides a non-trivial test of our implementation: Firstly, it checks whether the exclusion likelihood agrees over a wide range of different compositions of the gluon fusion and b quark associated Higgs production rates obtained in the MSSM parameter space, which are mapped onto the twodimensional likelihood grids (for fixed Higgs mass) in our reconstruction. Secondly, it tests whether our simple criterion of combining signal rates of Higgs bosons which have similar masses (overlapping within 20 %) is a reasonable approximation. Thirdly, the validation also tests whether the results obtained from the statistical hypothesis test of a single narrow resonance model can be mapped reasonably well onto the full neutral Higgs spectrum of the MSSM (and beyond).
Some deviations can be expected at the transition between regimes with different contributing Higgs combinations. As explained above, the implementation in HiggsBounds is based on the CMS results for the single narrow resonance interpretation, and the contributions of different Higgs bosons of a considered model can only be combined if their mass differences are such that they would appear as a single resonance in the CMS search. In contrast, in the dedicated CMS analyses carried out in specific MSSM benchmark scenarios it was possible to properly combine the contributions from different Higgs bosons at any given mass constellation since these have been simulated and tested with their particular masses at every parameter point. Therefore, the dedicated CMS analysis is expected to have a higher sensitivity than the HiggsBounds implementation if multiple Higgs bosons with different masses each give a non-negligible contribution to the signal yield. Furthermore, due to the simple criterion used in HiggsBounds for including/excluding the contributions of additional Higgs bosons, the considered rates in HiggsBounds may change quite abruptly in a transition region, where the selection of the tested Higgs boson combination changes. The single resonance approximation is expected to work best when the signal can be described as a single resonance formed by one or several Higgs bosons, while contributions of other Higgs bosons besides those associated with the resonance are negligible.
For predictions in the MSSM benchmark scenarios we employ the (M A , tan β) grids of Higgs production cross sections and branching fractions for the MSSM benchmark scenarios provided by the LHC Higgs Cross Section Working Group (LHCHXSWG) [28]. 3 For the gg → bb(h/H/A) production process we employ Santander-matching of the 4and 5-flavor scheme (FS) cross sections [57].
The results for the m max h scenario in the (M A , tan β) plane are shown in Fig. 2. 4 In Fig. 2a we show the distribution of the observed exclusion likelihood, q obs MSSM (in color), as obtained from HiggsBounds. The corresponding 95 % C.L. exclusion limit (orange, solid contour), which fulfills q obs MSSM = 5.99, is shown together with the CMS result obtained from a dedicated analysis in this benchmark scenario [20] (green, dashed contour). As mentioned in Sect. 2, the latter is based on the full combined 7+8 TeV dataset, whereas the exclusion information implemented in HiggsBounds is only based on the 8 TeV dataset. However, this fact is expected to lead to only minor differences in the excluded parameter regions. As can be seen, there is very good agreement between the exclusion limit reconstructed with HiggsBounds and the CMS result. Small deviations can be observed in the low M A region, M A 150 GeV, where all three neutral MSSM Higgs bosons contribute substantially to the signal yield. Here, the result reconstructed with HiggsBounds excludes a slightly smaller area of parameter space. The Higgs-Bounds result can thus be regarded as a conservative estimate of the actual exclusion limit.
In Fig. 2b we display, for every parameter point in the (M A , tan β) plane, the Higgs boson or combination of Higgs bosons (cluster) that has been selected to obtain the observed exclusion likelihood by the algorithm described in Sect. 3. It can be seen that all three neutral Higgs bosons are combined in most of the parameter region with M A 160 − 170 GeV and tan β 5, whereas at larger M A values only the two heavier Higgs bosons, H and A, which are nearly mass degenerate, are combined to yield the most sensitive constraint. At low tan β and large M A values, however, the combined signal rate of the heavier Higgs bosons becomes so small that it is instead the light Higgs boson, with mass around 120 − 125 GeV, that is selected to give the observed exclusion likelihood. This is because its expected exclusion likelihood is larger than that obtained for H/A. The observed exclusion likelihood obtained for the light Higgs boson with mass around 125 GeV is non-zero because the best-fit point in the two-dimensional cross section grid at m φ = 125 GeV is not identical with the SM prediction, cf. Fig. 1a. This leads to the small, but non-zero q obs MSSM values that are visible in Fig. 2a at large M A and small tan β.
Next we look at the light stop benchmark scenario, for which the cross section predictions and their associated theoretical uncertainties have been discussed in detail in Ref. [58]. This scenario features a relatively low SUSY particle mass scale, M SUSY = 500 GeV, and large stop mixing, X t = 2 M SUSY , leading to a lightest stop with a mass of ∼ 325 GeV. This leads to a reduction of the gluon fusion cross section of the light Higgs by around 10 − 15 % with respect to the SM prediction [26]. The results of applying our exclusion likelihood implementation in this scenario are shown in Fig. 3 (with colors similar to Fig. 2). The agree-  [20], as displayed in Fig. 3a, is very good for pseudoscalar Higgs masses M A 250 GeV. Similarly as in the m max h scenario the reconstructed exclusion limit obtained from HiggsBounds is slightly weaker for lower M A values than the CMS result from the analysis of the benchmark scenarios. As one can see in Fig. 3b, the reconstructed likelihood in the low-M A region obtained from HiggsBounds is mainly based on a combination of the H and A signals, while in most part of this parameter region the light Higgs boson at 125 GeV is not covered by the 20 % mass overlap criterion used in Higgs-Bounds. In contrast, in the dedicated CMS analysis in this scenario also the contribution from the light Higgs boson h is properly combined with the contributions from the other neutral Higgs bosons. The latter analysis therefore has a slightly higher sensitivity in this region, which means that the exclusion bound that we find here is slightly conservative compared to the dedicated CMS result. In addition to the excluded parameter region at values of tan β 5, the light stop scenario features an additional small excluded area at lower tan β values, namely tan β 2, and M A ∼ 145 − 190 GeV. The exclusion contour evaluated with HiggsBounds matches very well with the CMS result also in this region, where gluon fusion is the dominant production mode.
Finally, we test our implementation against the results obtained in the low-M H scenario, where the heavy CP-even Higgs boson is interpreted as the discovered SM-like Higgs boson at around ∼ 125 GeV and the light CP-even Higgs is largely decoupled from the SM gauge bosons [59]. Unlike the other benchmark scenarios, which use M A as a free parameter, this scenario is defined as a two-dimensional parameter plane in tan β and the Higgsino mixing parameter μ. The mass of the pseudoscalar Higgs, M A , is fixed to 110 GeV, which leads to a lightest Higgs mass that varies mostly between ∼80 and ∼105 GeV, but reaching even lower values at very low tan β and very high μ. Since in the MSSM a low value of M A implies also a light charged Higgs boson, this scenario served in particular as a benchmark for the LHC searches for charged Higgs bosons in top quark decays. In fact, the parameter space for this scenario in the MSSM is meanwhile essentially excluded [60-62] (see also Ref. [63]). Nevertheless this benchmark scenario is still very useful for our validation since all three neutral Higgs bosons have similar masses and thus contribute non-trivially to the analysis.
The comparison of the exclusion limits that have been reconstructed with HiggsBounds with the CMS results is shown in Fig. 4. It can be observed in Fig. 4a that there is rather good agreement between the exclusion limit obtained with HiggsBounds and the CMS result for μ values up to μ ∼ 2600 GeV. At this value of μ (depending on tan β) the reconstructed exclusion limit develops an "edge", and for higher μ values the reconstructed exclusion limit is significantly weaker than the one of the CMS result. The reason for this behavior is that at large μ the light Higgs mass becomes smaller than 88 GeV and is hence not combined with the heavier Higgs bosons A and H in HiggsBounds. This can be seen in Fig. 4b. In this parameter region, the tested signal rate is therefore significantly smaller than in the case of a full combination of h, H and A, and the resulting exclusion limit is accordingly weaker. In contrast, in the CMS analysis the signal yield of the light Higgs h has been properly taken into account even for very low mass values, and the possibly decreasing signal efficiency is partially compensated by the increasing production cross section, leading to the sig-nificantly stronger exclusion at high μ values obtained by CMS. 5 In almost all the remaining parameter space, all three Higgs bosons are combined by HiggsBounds, as can be seen in Fig. 4b, and the reconstructed exclusion contour resembles the CMS result for μ values below μ ∼ 2600 GeV. The slight deviations observed could result from mass-dependent selection efficiencies for the h and H signal yields, which cannot be accounted for in the Higgs-Bounds implementation since this information is not publicly available. Overall, even for this rather extreme scenario in the MSSM we find that the exclusion likelihood reconstructed with HiggsBounds approximates the results of a dedicated analysis reasonably well for large parts of the parameter space.

Example application: "Alignment without decoupling"
We now go beyond the validation with official CMS results and illustrate the usefulness of our exclusion likelihood implementation for another MSSM scenario. We consider here a scenario where the couplings of the light CP-even Higgs boson become SM-like for a certain range of tan β values, independently of the masses of the remaining Higgs spectrum. The existence of this so-called alignment limit was first pointed out in Ref. [64] for the 2HDM. After the Higgs discovery this possibility has gained attention through a series of papers [65][66][67][68][69], see also the "τ -phobic" benchmark scenario in Ref. [26]. In the MSSM the alignment limit can be realized independently of the decoupling of the heavier Higgs states through a cancellation between tree-level and higher-order contributions in the Higgs sector. This cancellation can occur at relatively large values of tan β and μ M S , with M S = √ mt 1 mt 2 being the stop mass scale. In the approximation tan β 1, and taking into account for simplicity only the dominant corrections at one loop, the alignment condition reads [69] tan β = (12) 5 A better agreement in the large μ parameter region could be obtained by increasing the mass overlap value of 20 % in the criterion for forming Higgs boson combinations to a sufficiently high value. However, firstly, there is no strong physics motivation to choose a value well beyond the quoted mass resolution of ∼ 20 % of the experimental τ τ analysis. Secondly, values larger than 20 % might lead to too aggressive exclusions in some scenarios. We therefore stick to the value of 20 % as the default setting.
Here, M Z and m t are the Z boson and top quark mass, respectively. M h denotes the light CP-even Higgs boson mass in the above approximation, and v ≈ 246 GeV. A t is the trilinear soft-breaking term in the stop sector.
Solutions of Eq. (12) with tan β > 0 exist if μA t (A 2 t − 6M 2 S ) > 0. Typically, in order to achieve the correct Higgs mass M h ∼ 125 GeV for not too large values of the stop masses, the stop mixing is chosen in the region where the prediction for M h is maximized, i.e. |X t | ∼ |A t | ∼ ± √ 6M S (at one-loop). Therefore, for μA t > 0 (μA t < 0), the alignment condition has viable solutions for values of |A t | that are slightly above (below) the value where the prediction for M h is maximized. By increasing |μA t /M 2 S | it is possible to lower the tan β value at which alignment occurs. An MSSM benchmark scenario of this kind for BSM Higgs searches at the LHC has recently been proposed in Ref. [69].
Here, we investigate the benchmark scenario proposed in Ref. [69], which is essentially a modification of the m mod+ h scenario [26] to allow for alignment independent of decoupling. This so-called m alt h scenario is defined by the parameter values (in the on-shell scheme) In contrast to the benchmark scenarios of Ref. [26], the parameters μ and M Q are adjustable parameters in the m alt h scenario. For convenience, the slepton, sbottom and first and second generation squark soft-breaking mass parameters are set to M Q , however, these can easily be adjusted to higher values in order to avoid constraints from SUSY searches at the LHC as their influence on the Higgs phenomenology is negligible here. We follow the suggestion to set M Q to 1 TeV per default and, if necessary, increase this value until a light Higgs mass of m h ≥ 123 GeV is obtained. In practice, this is only relevant at very low values of tan β in the benchmark scenario defined by Eq. (13). The parameter μ is then adjusted according to a chosen ratio μ/M Q . We focus here on the choice μ/M Q = 3, implying rather large values of μ, where alignment independent of decoupling occurs at tan β ∼ 10.
The MSSM predictions are obtained using the public computer codes FeynHiggs-2.10.2 [31][32][33][34]36] for the Higgs masses, couplings and branching fractions, and SusHi-1.4.1 [30] for the gluon fusion and b quark associated production cross sections of the three neutral Higgs bosons.
The numerical results for this benchmark scenario are displayed in Fig. 5. The observed exclusion likelihood from the CMS φ → τ τ search as obtained from HiggsBounds is shown in color in Fig. 5a, and  parison, the green contour shows the exclusion line obtained in Ref. [69] using results from the same CMS analysis, however, following a more simplistic approach. 6 As can be seen from the figure, the more advanced implementation of the observed exclusion likelihood in HiggsBounds leads to a somewhat stronger 95 % C.L. exclusion limit over most of the parameter space. The relative behavior seems to be different in the tt threshold region, M A ≈ 2m t ≈ 345 GeV, where in particular the gg → A cross section is enhanced. However, the approximations made in Ref. [69] appear to be least reliable in this region. The gray dotted line shows the exclusion limit obtained by CMS for the m mod+ h scenario. As discussed in Ref. [26], the excluded regions in the benchmark scenarios are significantly affected if decay modes of the heavy Higgs bosons H and A into SUSY particles are kinematically open and unsuppressed. The presence of such decay modes leads to a sizable reduction of the H/A → τ τ branching fractions and therefore to a smaller excluded region. In the alignment scenario μ is very large, leading to a negligible Higgsino component in the light neutralinos and chargino. The branching fractions for the Higgs decays to neutralinos and charginos are therefore essentially absent. In addition, the heavy Higgs decays to gauge bosons, H → W + W − and H → Z Z, are also suppressed, as the responsible coupling ∝ cos(β − α) vanishes in the alignment limit. As a result, the H/A → τ τ branching fraction is significantly higher in the alignment scenario than in the m mod+ h scenario, which leads to a much larger excluded region in the alignment scenario, see also the discussion in Ref. [69].
In order to illustrate the complementarity between the constraints from the CMS φ → τ τ search and the constraints obtained from the signal rate measurements of the discovered Higgs boson, we show in Fig. 5b the likelihood distribution, χ 2 HS , obtained from a χ 2 test of the light Higgs boson signal rates against a combination of the latest rate measurements from the LHC [70-78] and the Tevatron [79,80], using the public computer code HiggsSignals-1.3.0 [13] (see also Refs. [19,81]). The 95 % C.L. preferred region lies within the orange contours in Fig. 5b. It is given by the χ 2 difference with respect to the minimal χ 2 value (located in the alignment region and indicated as gray asterisk in Fig. 5b), χ 2 HS ≡ χ 2 HS −χ 2 HS,min ≤ 5.99. It can be seen that the χ 2 distribution becomes independent of M A at around tan β ≈ 10, indicating that the couplings of the light Higgs become SMlike independently of the decoupling of the heavier Higgs states.
Since we now have the exclusion likelihood q obs MSSM from the CMS φ → τ τ search available, we can perform a statistical combination with the constraints from the Higgs rate measurements by constructing the global χ 2 function χ 2 tot = q obs MSSM + χ 2 HS . The resulting χ 2 tot distribution 7 is shown in Fig. 6. The constraints from the φ → τ τ searches at the LHC are highly complementary to the rate measure-7 Again, χ 2 tot is the χ 2 difference with respect to the minimal χ 2 value (obtained at M A = 500 GeV, tan β = 4, i.e. in the lower right corner of Fig. 6), now based on the global likelihood χ 2 tot . The contours indicate the 1σ , 2σ and 3σ allowed regions ments, since they are particularly sensitive at higher values of tan β where the production process gg → bbφ is enhanced. In the m alt h scenario with μ = 3M Q , the combination of both constraints yields a lower limit of M A 350 GeV at the 95 % C.L. Thus, alignment of the light Higgs boson occurring without the simultaneous decoupling of the heavier Higgs states is ruled out for this scenario. The alignment without decoupling limit can be pushed to lower values of tan β in this scenario, where the constraints from the φ → τ τ searches are less significant, only by choosing even more extreme values of μA t /M 2 Q , which potentially leads to problems with vacuum stability [82,83].

Conclusions
LHC searches for non-standard Higgs bosons decaying into tau lepton pairs constitute a sensitive experimental probe for BSM physics. Recently, the CMS collaboration published the likelihood information for their Higgs boson searches in the τ + τ − final state [20]. The likelihood is given as a function of the two relevant Higgs production channels, gluon fusion and b quark associated production, for various mass values of the narrow resonance assumed for the signal model. In this paper we have shown how this experimental information can be utilized to test large classes of theoretical models. In particular, we have developed a simple algorithm that maps an arbitrary model with multiple neutral Higgs bosons onto a model with a single narrow resonance, for which the corresponding exclusion likelihood from the CMS search can be determined. We have described the inclusion of this method into the new version of the publicly available For-tran code HiggsBounds (version 4.2.0 and higher). For nearly any model under consideration, HiggsBounds provides an evaluation of the exclusion likelihood for a model parameter point based on the information from Ref. [20]. Similarly, if requested, HiggsBounds can also perform a test of whether or not a given parameter point is excluded at the 95 % C.L. based on all available searches, including the new τ + τ − result. The approach to test BSM models with exclusion limits is complementary to testing the compatibility of a given model with the observed Higgs signal (and possible future signals of additional Higgs bosons). The latter kind of information is contained in the sister program HiggsSignals, and both programs can be used together in order to obtain the combined likelihood information from the search limits and the Higgs rate measurements. Both codes are available at http://higgsbounds.hepforge.org.
We have validated our implementation of the τ + τ − search results into HiggsBounds by comparing the 95 % C.L. exclusion contours obtained with HiggsBounds with the ones obtained by CMS from dedicated analyses in three Higgs benchmark scenarios [26] in the MSSM. We found very good agreement in the parameter regions where the sensitivity of the search is dominated by a single combination of Higgs bosons that can be identified with a single narrow resonance assuming an experimental mass resolution of 20 %. As expected, the largest but still relatively small deviations occur in parameter regions where all neutral MSSM Higgs bosons are relatively close in mass and contribute comparably to the signal yield.
As an application, we have discussed the combined constraints of the τ τ search and the rate measurements of the SM-like Higgs at 125 GeV in a recently proposed MSSM benchmark scenario, where the lightest Higgs boson obtains SM-like couplings independently of the decoupling of the heavier Higgs states. Here we combined the χ 2 analysis of the rate measurements for the Higgs signal, evaluated with HiggsSignals, with the exclusion likelihood from the non-observation in the τ + τ − search channel, evaluated with HiggsBounds. We have shown that the combined information yields very significant constraints on the available parameter space in this scenario and in fact disfavors the "alignment without decoupling region" in the studied benchmark model.
We encourage ATLAS and CMS to continue providing their search results including the relevant likelihood information. This will greatly facilitate the application of the search results for testing BSM models.

Appendix: User guide: how to obtain the exclusion likelihood with HiggsBounds
The exclusion likelihood information for the CMS φ → τ τ analysis [20] is implemented in HiggsBounds from version 4.2 on. As described in Sect. 3 this information is used in a standard HiggsBounds run to reconstruct the expected and observed 95 % C.L. exclusion limit, which is then considered alongside all other available Higgs search limits in the full HiggsBounds application. This leads to the global information whether the tested parameter point is allowed or excluded at the 95 % C.L. Beyond this information the value for the exclusion likelihood for the model parameter point under investigation, q model , can also be obtained directly via HiggsBounds Fortran subroutines, enabling the user to incorporate this information e.g. in a global parameter fit. In the following we document the relevant subroutines that make this information accessible.
The main routine that runs the algorithm presented in Sect. 3 to obtain the exclusion likelihood is: HiggsBounds get likelihood(int analysisID, int Hindex, int nc, int cbin, dble M, dble llh, char( * ) obspred) The (mandatory) input argument 8 analysisID specifies the analysis for which the likelihood should be obtained. At the moment, the CMS φ → τ τ analysis based on the full 8 TeV dataset (analysisID = 3316) is the only available likelihood, but the framework is easily extendable for future experimental results. The output values provide information about the selected Higgs boson combination (or Higgs cluster): -Hindex gives the index i of the Higgs boson h i , which provided the initial seed to form the dominant Higgs cluster (cf. Sect. 3, item 1), 8 Here, and in the following, input arguments and optional input arguments are highlighted in dark blue and green, respectively. The remaining arguments are output values.
nc gives the number of Higgs bosons contained in the combination, cbin is a binary code (bitmask) for the identifiers of the participating Higgs bosons. The binary code is given by summing over 2 (i−1) for all involved Higgs bosons. 9 For example, in the MSSM (with common indexing h 1 = h, h 2 = H , h 3 = A), the combination H +A would give cbin = 6, whereas a cluster formed only by the light Higgs h gives cbin = 1.
The output value M gives the averaged mass value, calculated according to Eq. (8), at which the likelihood value has been evaluated. The computed value of the likelihood, q model , is returned as llh. The final argument is an optional input, obspred, which takes a string value that can be either 'obs' or 'pred', specifying whether the observed or expected (or predicted) exclusion likelihood should be evaluated, respectively. The default behavior if this argument is not provided is that the routine returns the observed likelihood, following the algorithm described in Sect. 3.
In addition to the main subroutine we provide the following two auxiliary routines, which may be helpful to understand the obtained results. Table 1 Implemented 95 % C.L. exclusion limits from LHC searches for BSM Higgs bosons with τ τ final states in HiggsBounds-4.2. The analysisID is used as a unique identifier for an individual analysis and can be used to deactivate/activate them in HiggsBounds (see text) analysisID Experiment Luminosity and CM-Energy Additional notes References 3316 CMS 19.7 fb −1 at 8 TeV Using −2 ln L reconstruction [20] 2014049 ATLAS 19.5 − 20.3 fb −1 at 8 TeV Profiled limit on gg → bbφ process [21] 20140492 ATLAS 19.5 − 20.3 fb −1 at 8 TeV Profiled limit on gg → φ process [21] run, since these are better described by the likelihood information. In particular, the 95 % C.L. limits from previous BSM φ → τ τ searches can be deactivated if instead the CMS φ → τ τ exclusion likelihood is used. In order to do so, we provide two new subroutines, contained in the Fortran module 'channels'.

HiggsBounds deactivate analyses(int(:) analysisID list)
This routine should be called before the subroutine run_higgsbounds in order to deactivate the analyses specified by the integer array analysisid_list. for convenience, the analysis identifiers of the currently implemented lhc φ → τ τ searches are given in Table 1.

HiggsBounds activate all analyses()
This subroutine can be used at any time to re-activate all previously deactivated analyses for the succeeding Higgs-Bounds run.
In order to demonstrate the use of these subroutines, we provide the example program HBwithLHClikelihood, included in the /example_programs/ directory of the HiggsBounds distribution. This program shows how to obtain the observed exclusion likelihood from the CMS φ → τ τ analysis in the MSSM m max h scenario, such that the user should be able to directly reproduce Fig. 2. The example can be compiled by calling 'make HBwithLHClikelihood' in the HiggsBounds main folder, and run from the example_programs folder by calling './HBwith-LHClikelihood'. Following a successful run, the gnuplot scripts 'plot_mhmax_llh.gnu' and 'plot_ mhmax_llh_comb.gnu' in the same folder then reproduce Fig. 2.