HiggsBounds-4: Improved Tests of Extended Higgs Sectors against Exclusion Bounds from LEP, the Tevatron and the LHC

We describe the new developments in version 4 of the public computer code HiggsBounds. HiggsBounds is a tool to test models with arbitrary Higgs sectors, containing both neutral and charged Higgs bosons, against the published exclusion bounds from Higgs searches at the LEP, Tevatron and LHC experiments. From the model predictions for the Higgs masses, branching ratios, production cross sections and total decay widths - which are specified by the user in the input for the program - the code calculates the predicted signal rates for the search channels considered in the experimental data. The signal rates are compared to the expected and observed cross section limits from the Higgs searches to determine whether a point in the model parameter space is excluded at 95% confidence level. In this document we present a modification of the HiggsBounds main algorithm that extends the exclusion test in order to ensure that it provides useful results in the presence of one or more significant excesses in the data, corresponding to potential Higgs signals. We also describe a new method to test whether the limits from an experimental search performed under certain model assumptions can be applied to a different theoretical model. Further developments discussed here include a framework to take into account theoretical uncertainties on the Higgs mass predictions, and the possibility to obtain the $\chi^2$ likelihood of Higgs exclusion limits from LEP. Extensions to the user subroutines from earlier versions of HiggsBounds are described. The new features are demonstrated by additional example programs.


Introduction
The search for Higgs bosons [1][2][3][4][5][6] is, and has been, a major cornerstone of the physics programmes of past, present and future high energy colliders. This has become even more important in view of the recent discovery of a Higgs signal by ATLAS [7,8] and CMS [9][10][11]. Determining the properties of this newly observed state and comparing the measurements to explicit theories beyond the Standard Model (SM) is one of the present challenges. These theories often contain enlarged Higgs sectors with multiple Higgs bosons. Even in the presence of a signal, it is therefore important that the LEP, Tevatron and LHC experiments present exclusion limits from the non-observation of Higgs bosons in various channels. These are very useful for constraining the available parameter space of those models which are able to fit correctly the observed Higgs signal. Such constraints will need to be taken into account also in the future interpretation of the Higgs results in the context of models of new physics.
In this paper we describe new developments in version 4 of the publicly available Fortran code HiggsBounds [12][13][14], which has been designed for exactly this purpose. For the complementary approach, to test whether a model is compatible with the observed LHC Higgs signal (and possible future signals of additional Higgs bosons), we have also developed the sister program HiggsSignals, which has been described in Ref. [15]. It is highly recommended to use these two programs in parallel to obtain the most complete test for extensions of the SM Higgs sector.
The experimental analyses implemented in Higgs-Bounds usually take one of two forms. Dedicated analyses have been carried out in order to constrain some of the most popular models, such as the SM [7][8][9][10][11]16] and various benchmark scenarios in the Minimal Supersymmetric Standard Model (MSSM) [17][18][19][20]. In addition, model-independent limits on the cross sections of individual signal topologies (such as e + e − → h i Z → bbZ) have been published. In the former type of analyses several search channels (or signal topologies) are typically combined in order to maximize the discovery/exclusion reach. However, the re-interpretation of these results in the context of different models than those already investigated by the search analysis requires detailed knowledge of the individual efficiencies (or signal contaminations) of the investigated search channels. In contrast, the latter type of analysis can be used easily to test a wider class of models.
HiggsBounds-4 has been designed to facilitate the task of comparing Higgs sector predictions with existing exclusion limits, thus allowing the user to quickly and conveniently check a wide variety of models against the state-of-the-art results from Higgs searches. Version 4 differs significantly from previous versions of the code (described in [12,13,21]) in several respects. The code now fully supports testing models against exclusion limits from the LHC, which are implemented for analyses performed at center-of-mass energies of both √ s = 7 and 8 TeV. The main algorithm of HiggsBounds has been extended to ensure a reliable application of exclusion limits in the presence of a signal (as is now observed in the LHC data). The model-likeness test, which tests whether a given model fulfills the assumptions of a particular Higgs search to a sufficient degree, has been fully rewritten to enable in particular the limits from SM Higgs searches at the LHC to be applied in a wider context. We introduce an option to take into account theory uncertainties on the Higgs mass predictions, which are relevant, for instance, for the lightest Higgs boson mass in the MSSM. An alternative statistical treatment for the LEP constraints (in the form of a χ 2 output) is provided. Finally, we describe an improved input/output for supersymmetric (SUSY) models that can now be given in the SLHA format [22,23]. The main focus of this updated documentation is to provide a detailed description of these new developments, to show relevant physics examples of where improvements can be expected, and to introduce the user to how the improved HiggsBounds code can be used in practice.
This paper is structured as follows. In Sect. 2 we give a general introduction to the statistical approach employed in the HiggsBounds code, and describe, in particular, the way in which this approach has been extended in Higgs-Bounds-4. Section 3 gives a thorough description of the different methods of providing theory input for Higgs-Bounds, and their extension to LHC7/8 predictions. This is followed by Sect. 4, which contains a discussion of the major new developments, including numerical examples. Finally, Sect. 5 contains the technical details on how these new features can be used in practice, extending the original Higgs-Bounds manual [12,13] with a description of the new subroutines, data files and example programs that are provided. In an Appendix we list and provide references to all the experimental analyses that provide results implemented in the code, including the analyses added in the latest public version (HiggsBounds-4.1).

General approach of HiggsBounds-4
The general concept of HiggsBounds, including details on the treatment of limits from LEP and the Tevatron, has already been published in [12][13][14] (see also Ref. [21]). From a conceptual point of view, the extension of HiggsBounds to include LHC limits is straightforward. The technicalities of this implementation, and how it modifies the user input, is discussed in Sect. 3. Our aim here is to give a brief introduction to the purpose of the code and the methods it uses. We also introduce one conceptual change with respect to previous versions, which has been prompted by the application of HiggsBounds to models which feature a Higgs boson with a mass close to the observed LHC Higgs signal.
The basic input for HiggsBounds (which the user has to provide) are the relevant physical quantities predicted for the Higgs sector of the model under consideration. The necessary predictions for each Higgs boson H i (i = 1, . . . , n H 0 +n H ± ) in the model are, schematically, i.e. the Higgs boson mass, its total decay width (it is assumed that the narrow width approximation holds), its decay branching ratios, and the production cross sections, normalized to a particular reference value. Here, P denotes a specific Higgs production process. If P exists in the SM, its cross section, σ SM (P(H )), evaluated at the same mass value, M H = M H i , is typically used as the reference cross section, σ ref . In some cases it can also be necessary to supply additional predictions, such as the BR(t → bH + ), or the CP properties of the neutral Higgs bosons. Variations on the input format are offered, which allow the user to specify a simpler set of input quantities when certain basic approximations are valid. A complete list of the options for giving model input is given in Sect. 3.
In addition to the model predictions, the other important ingredient of HiggsBounds is the experimental data. Exclusion limits from (negative) Higgs searches are collected from the experimental publications, with the aim of keeping the code as up-to-date as possible with the latest developments. Currently the code includes results from LEP, the Tevatron and the LHC experiments. More information on which analyses are available in HiggsBounds is provided in Appendix 7.1. The data for these analyses is contained in tables of expected exclusion limits at 95 % C.L. in the absence of a signal (based on Monte Carlo simulations), and the corresponding observed limits, as a function of the Higgs boson mass. The list consists both of analyses for which modelindependent limits were published, and of dedicated analyses carried out specifically under the assumption of the SM (like most LHC searches to date), or for Higgs bosons with certain CP properties. These limits can be applied to models with Higgs bosons which show these characteristics to a sufficient degree 1 at the parameter point being considered.
We call the application of the limit from a particular Higgs search to one of the Higgs bosons of the model under study (or to two of the Higgs bosons, for searches involving two Higgs bosons) an "analysis application", which we denote by X in the following. 2 Each analysis application has a correspond-ing signal cross section prediction σ (X ), which Higgs-Bounds uses to calculate the relevant quantity Q model (X ) for which the experimental limit is given; typically this is a conveniently normalized cross section times a branching ratio. The corresponding experimental quantities are denoted Q expec (X ) and Q obs (X ) for the expected and observed limits, respectively. If two Higgs bosons have a narrow mass separation, δ M = M h i − M h j , then their predicted cross sections are added for certain analyses where the mass resolution is limited and interference effects are expected to be negligible. The settings for the maximal δ M h can be varied by the user separately for LEP, Tevatron, and LHC analyses (the default values are 0 GeV for LEP and 10 GeV for Tevatron/LHC).
HiggsBounds operates by considering, for each analysis application, the ratio of the model predictions, Q model (X ), to the experimental limits. To ensure that the result can be interpreted as an exclusion at 95 % C.L. (which is the same confidence level as adopted by the individual analyses), it is crucial that the model prediction is only compared to the experimentally observed limit for one particular analysis application. In a first step, HiggsBounds therefore uses the expected experimental limits to determine the analysis application X 0 with the highest statistical sensitivity to exclude the model point under consideration, In the second step, HiggsBounds then performs the exclusion test for the Higgs boson and analysis combination represented by X 0 , by computing the ratio to the observed limit If k 0 > 1, HiggsBounds concludes that this parameter point of the tested model is excluded at 95 % C.L. 3 The statistical method as described here (in the following referred to as the classic method) has been the only mode of operation available in previous HiggsBounds versions. For HiggsBounds-4, we have extended this method to Footnote 2 continued Higgs searches, A 1 and A 2 . Then there are six possible analysis applications, 3 If we had instead compared the predicted cross sections directly to the experimentally observed limits for all available search channels and considered the model excluded if at least one of them gave exclusion at 95 % C.L., the result would in general not correspond to an exclusion at 95 % C.L. The combined probability of yielding a false exclusion from any of the individual comparisons of Q model to Q obs would also yield an overall probability for false exclusion higher than that from applying a single limit. perform better in situations where a Higgs boson signal is present (as is now the case in the LHC data). The problem of the classic method arises for models with multiple Higgs bosons. If one of these has a mass close to that of the observed signal (which is likely, since any reasonable model should also explain this signal), its analysis applications will test the model predictions against limits (for various channels) in the signal region. In this region, the expected limits (based on the background-only hypothesis) will continue to improve with more experimental data and optimized analysis methods, whereas the observed limits can never be expected to reach exclusion at the SM level (provided a true signal of near-SM strength is what is observed). For model points where the most sensitive analysis application X 0 is a test of the signallike Higgs boson, the classic HiggsBounds method would therefore never yield exclusion. Moreover, constraints on the remaining Higgs spectrum (with less expected sensitivity) are not applied. Even if the exclusion remains formally valid at 95 % C.L., it could be anticipated that this problem would eventually become serious enough to limit the usability of the code.
Among the several possible ways that the HiggsBounds algorithm could be extended to address this problem, all involving different compromises, we have opted for a solution which involves a slight violation of the strict testing of only one experimental limit. We call this the full HiggsBounds method. In summary, this method performs the original HiggsBounds test separately for each Higgs boson in the model. In the full HiggsBounds method, the first step is to evaluate the most sensitive analysis application X i for each Higgs boson H i according to This is followed by a straightforward exclusion test on the individually most sensitive analysis applications The result of these tests contains more information than the single test of HiggsBounds classic (such as exclusion/nonexclusion by individual Higgs bosons), which is now made available to the user (see Sect. 5 for details). A combined HiggsBounds exclusion is also calculated, where the result is interpreted as model exclusion if k i > 1 for any of the k i . The combined (single-number) output is then calculated as By the construction of the full method, it follows directly that the two methods are equivalent for models with a single Higgs boson. It is also clear that the full method can only give stronger exclusion than the classic method. This is consistent with the fact that the exclusion of the full method will correspond to a limit at somewhat lower statistical confidence level than 95 %. Still, the deviation from the strict 95 % C.L. should be minor in this approach compared to the alternative (naive) testing of all Higgs bosons versus all observed limits, since the number of Higgs bosons in a model in general is much smaller than the number of implemented experimental analyses. Furthermore, a non-negligible dilution of the 95 % C.L. interpretation of the combined result is only expected in the case where more than one test X i leads to a ratio k i ≈ 1.
To illustrate the difference between the classic and full methods of HiggsBounds, we show in Fig. 1 figure). The exclusion bounds, as evaluated by Higgs-Bounds, are shown separately for LEP exclusion (blue) and the LHC (red). When evaluating the limits in this figure, a theory uncertainty of 3 GeV is taken into account in the evaluation of the lightest Higgs mass, see Sect. 4.2 for details , and H ± (right). The colour coding is the same as in Fig. 1 on how this is done. As can be seen from this figure, the full method gives the strongest exclusion, corresponding to the most accurate application of the existing limits in this scenario (as also used in [24]). The difference to the classic method can be seen in particular for high M A and high tan β (the decoupling regime). Here the applicability of the classic method is limited, since the globally most sensitive channel is a search for the lightest (SM-like) Higgs boson, which cannot be excluded when is mass its in the signal region, M h 125 GeV. This is in contrast to the results in the full method, which can be further illustrated by looking at the contribution of individual Higgs bosons as shown in Fig. 2 for the same MSSM example. The first panel shows the exclusion contributed by h. The narrow unexcluded region around M A = 135 GeV results from a particular channel ( pp → V H, H → bb) being the most sensitive. For this channel, the observed limit is not strong enough here to exclude the lightest Higgs. The second panel shows the exclusion for H/A. They are treated together, since their masses are close to degenerate over most of the parameter space. The dominant exclusion therefore comes from the same search channels and their signal rates are added. Finally, the last plot shows the exclusion from H ± . The exclusion region presented for the full method in Fig. 1 consists of the union of the three different exclusion regions shown here. In the HiggsBounds distribution we provide an updated example program, HBwithFH, which can be used to test MSSM parameter points for exclusion using either the full or the classic HiggsBounds methods.

Theoretical predictions
The theoretical model predictions, which are compared to the experimental data in the analysis applications, are computed from the user input. A detailed explanation of this input has been given for previous version of HiggsBounds [12,13]. For completeness, we include the full input specification here, describing both the original input and the updates in HiggsBounds-4. For the input of the required theory predictions to HiggsBounds, the user can choose between four different formats. In the command line version of the program, these are labelled by the variable whichinput, which can take the values hadr (hadronic), part (partonic), effC (effective couplings), or SLHA (SUSY Les Houches Accord). The type of input required for each of these formats is summarized in Table 1, and the following subsections contain more detailed information about each type of input. In particular, we describe the internal treatment of LHC cross sections (at 7 and 8 TeV). Technical details on the subroutines which should be used for the different input options can be found in Sect. 5.

Hadronic input
The hadronic input option whichinput=hadr requires model input in the most general form and therefore contains the lowest degree of approximation. The complete input in this format involves the specification of: (1) The masses for the neutral Higgs bosons, h k (k = 1, . . . , n H 0 ), and the (singly) charged Higgs bosons, H ± k (k = 1, . . . , n H ± ), in the model.
their total decay widths, the neutral Higgs branching ratios without SM equivalents the charged Higgs branching ratios to SM particles the top quark branching ratios normalized ratios R σ (P) (see below) of hadronic Higgs production cross sections at the Tevatron for the processes and, finally, normalized ratios R σ (P) of hadronic Higgs production cross sections at the LHC (both for √ s = 7 and 8 TeV, given as separate inputs) for the processes It is important to note that this corresponds to an exhaustive list of the possible input. In certain (most) cases only a subset of these inputs may be required. For example, if the user is only interested in limits from the LHC, no LEP or Tevatron cross sections need to be given as input (the corresponding values can be set to zero). Likewise, for models with only neutral Higgs bosons, no input involving either the charged Higgs sector or top decays will be used. The meaning of most of these quantities should be pretty clear already from the notation; for those that require further clarifications we provide some more details below.
For input items (8), (9) and (10), the normalized cross section of a Higgs production process P is simply defined by Where a SM equivalent exists, the reference cross section , where e denotes the electromagnetic coupling constant, and s w (c w ) the sine (cosine) of the weak mixing angle, respectively. 4 Similarly, for the process P = e + e − → H + j H − j , the reference cross section is the cross section of the process e + e − → H + H − in a two-Higgs-doublet model (e.g. the MSSM) at tree-level (i.e. with s-channel γ and Z exchange). This reference cross section depends solely on the mass of the charged Higgs boson and SM quantities. As a consequence, the leading-order prediction in the MSSM is R σ (e + e − → H + j H − j ) = 1.
The branching ratio to "invisible", BR model (h k → invisible), is defined as the branching ratio of a neutral Higgs boson into particles which are only infered in the detector by their contribution to the missing transverse energy. Examples of this includes the lightest neutralino in the MSSM [25,26], inert scalars [27][28][29], or majorons in supersymmetric models with spontaneous breaking of R-parity [30][31][32].
The hadronic input is the most versatile, since it allows the user to provide the predictions in the form of the most precise calculations available. It is also the only input format which allows for studying e.g. effects of parton distribution functions or hadronic uncertainties on the Higgs exclusion bounds. On the other hand, this input format is also the most demanding, and in order to make it more convenient for the user to correctly normalize the rate predictions, the Higgs-Bounds library provides a series of Fortran functions which allow the user to access the predictions of certain SM quantities, including the hadronic SM Higgs production cross sections, total decay width, and branching ratios as a function of the Higgs mass.

Parton-level input
The second possibility for specifying the HiggsBounds input (with whichinput=part) is using ratios of partonic cross sections as far as is possible. This input format is in many cases more convenient for the user to calculate than the hadronic option. It requires (at most) the following model predictions (1)-(8) the same as for whichinput=hadr, (8) It should be calculated for a parton-system center-of-mass energy squaredŝ =ŝ 0 , whereŝ 0 denotes the partonic production thresholdŝ 0 = (M h j + M y ) 2 . For this approximation to be valid the dependence of the partonic cross section on s is required to be mild. For the case of single Higgs boson production, M y = 0.
The partonic cross section ratios R h j +y nm can be a lot easier to calculate than hadronic cross section ratios for a wide range of models. In addition, it is often (at least to a good approximation) the case that for all nm. This reduces the number of partonic cross section ratios which need to be provided by the user from 13 to 5. An example of a model of this type is given by the MSSM with real parameters, where the common ratio can be calculated approximately from the normalized (squared) effective coupling of the Higgs boson to a pair of Z bosons, For a more complete discussion of the use of the efffective coupling approximation as input to Higgs-Bounds, see Sect. 3.3.
Internally, HiggsBounds uses the approximate relations to calculate the missing hadronic cross section ratios from the partonic cross section ratios. Here, σ SM ( pp → nm → H + y) denotes the contribution from the partonic initial state nm to the total hadronic cross section for the process pp → H + y in the SM. The hadronic LHC ratios are evaluated separately for 7 and 8 TeV but, as already mentioned, using the same values for R h j +y nm . The hadronic cross section ratios for the SM appearing in Eqs. (9), (10) are provided within HiggsBounds. Further discussion of the applicability of this approximation, and details of how the ratios σ SM ( pp → nm → H + y)/σ SM ( pp → H + y) are calculated for the Tevatron have been presented in [12]. In HiggsBounds-4, these cross section ratios are provided for the LHC with √ s = 7 TeV and √ s = 8 TeV based on the prediction of the LHC Higgs Cross Section Working Group [33][34][35] for the gluon fusion cross section and the H Z and H W cross sections, and the program bbh@nnlo 1.3 [36]. Results on

The effective coupling approximation
With the effective coupling option (whichinput=effC), the user input is greatly simplified and reduced to a smaller number of quantities. From this input, physical predictions corresponding to input with the partonic option whichinput=part are calculated. The effective couplings involve specifying (at most) (1)-(2) the same as for whichinput=hadr, (3) normalized scalar and pseudoscalar (squared) effective Higgs couplings to fermions

normalized (squared) effective Higgs couplings to bosons
neutral Higgs branching ratios without SM equivalents, charged Higgs branching ratios to SM particles, and top quark branching ratios according to (5)-(7) of whichinput=hadr.
The scalar and pseudoscalar components of a Higgs coupling to a fermion pair are defined in the usual way, via the Feynman rule for the coupling of a generic neutral Higgs boson h to fermions: where g s and g p are real-valued scalar and pseudoscalar coupling constants respectively, and 1 and γ 5 are the usual matrices in Dirac space. A CP-even scalar, like the SM Higgs boson, has g p = 0 and a CP-odd scalar has g s = 0. Where these exist, the reference couplings are taken as the SM tree-level equivalents: g SM where M Z is the Z boson mass, M W the W boson mass, and m f the mass of fermion f . The reference coupling (g ref H H Z ) 2 , that does not have a SM equivalent, is defined as The effective couplings (g model h k γ γ /g SM H γ γ ) 2 (and similarly for γ Z ) are loop-induced. They can most conveniently be defined via For the Higgs-gluon-gluon effective coupling, (g model there is a choice of definition. It can be defined either via decay widths as or via partonic cross sections: It has to be understood that both these definitions will involve approximations. It is therefore only recommended to use the input of effective couplings when both definitions result in similar values for (g model h k gg /g SM Hgg ) 2 . However, under certain circumstances, this condition can be relaxed. For example, in a model in which the LEP searches for Higgs bosons decaying into hadrons are not relevant, the effective h k gg coupling can be defined solely by Eq. (18). Conversely, if for some reason the gluon fusion Higgs production mechanism is not relevant, the effective coupling can be defined solely by Eq. (17).
The calculation of LEP and Tevatron cross section ratios from the effective couplings has been described in [22,23].
Here we shall focus on the extension of this procedure to LHC cross sections. Partonic cross section ratios for the LHC are calculated as where (q, q ) ∈ {(u, d), (c, s)} and q ∈ {u, d, c, s, b}. The normalized hadronic cross section ratio for tt production together with a CP-even Higgs boson for the LHC is obtained using For h k production via VBF, the normalized hadronic cross section ratio is calculated using the approximate relation The M H dependence of the relative fractions of VBF events induced from W W and Z Z fusion, denoted as R W W VBF,LHC and R Z Z VBF,LHC , respectively, is mild and, for the Tevatron case, can be approximated by constant values. For the LHC case, we obtain these functions by fitting to SM results produced with HAWK 1.1 [37,38] for pp collisions at √ s = 7 TeV. Note that, for models where (g model independent of R W W VBF,LHC and R Z Z VBF,LHC . The input scheme for decay widths and branching ratios is not affected by the extension to include LHC results in HiggsBounds-4, and the calculation from the effective couplings follows what has been published earlier [12,13]. 5 The only difference compared to previous versions is that we have updated the SM reference values to agree with those recommended by the LHC Higgs cross section working group [33][34][35] for M H between 80 GeV and 1 TeV.
In order to constrain Higgs bosons with masses below ∼ 10 GeV, the effective coupling input option is usually not appropriate because the implemented SM reference branching ratios are rather inaccurate for such low masses. It is therefore strongly recommended to use one of the other input formats and enter the branching ratios for such a light Higgs boson directly. It can also be relevant in this mass region to consider constraints from other colliders than LEP, which are not included in HiggsBounds.

Input using the SUSY LesHouches Accord
For the convenience of users interested in supersymmetric models, an input option using the SUSY Les Houches Accord (SLHA) [22,23] is now offered. This option (available by setting whichinput=SLHA) uses the calculated decay information from the SLHA file and cross sections approximated through the effective couplings approach. It requires the following input to be read from an SLHA file: (i) The masses for the neutral Higgs bosons h k (k = 1 . . . n H 0 ) and singly charged Higgs bosons (ii) the Higgs total decay widths, (iii) neutral Higgs branching ratios with SM equivalents Note that the decay modes have to be specified as twobody decays, irrespectively of whether they are on-or off-shell. (iv) the neutral Higgs branching ratios without SM equivalents (v) the charged Higgs branching ratios to SM particles (vi) the top quark branching ratios (vii) normalized scalar and pseudoscalar (squared) effective Higgs couplings to fermions

(viii) normalized (squared) effective Higgs couplings to bosons
In the SLHA input the effective couplings are only used to calculate the Higgs production cross section ratios (unlike the effective coupling option, where they are also used to calculate the branching ratios). The Higgs decay branching ratios are taken directly from the corresponding decay blocks in the SLHA file. In the case of incomplete input, Higgs masses which are not specified in the SLHA file will be set equal to minus one (such that these Higgs bosons are not tested against any limits), whereas any other input that is not specified will be set equal to zero. Table 2 lists the PDG codes of the particles that can be considered as neutral Higgs bosons by HiggsBounds. The setting of nHzero determines how many of these are used, starting from the top of Table 2. For example, with nHzero=3, the properties of particles with the PDG numbers 25, 35, and 36 are read from the SLHA file. Note that no CP properties are inferred from the PDG numbers of the neutral Higgs bosons. The invisible Higgs branching ratios are obtained from the branching ratios of Higgs bosons into a weakly-interacting lightest supersymmetric particle (LSP).
HiggsBounds finds the weakly-interacting candidate with the lowest mass (considering neutralinos and sneutrinos as candidates) and confirms that this particle is indeed the LSP by comparing its mass against the masses of the charged leptons, the lightest chargino, and the gluino. If the LSP is not neutral, the invisible Higgs branching ratio is set to zero.
To specify the required effective couplings, as described by points (vii) and (viii) above, HiggsBounds requires two extra input blocks which are not part of the normal SLHA. An example of these blocks is shown in Table 3. There are some cases when HiggsBounds is unable to use an SLHA file, including any of the following: • The Block MODSEL indicates that there is R-parity violation, • either Block SPINFO or DCINFO contains an entry with the label '4' (which is used to indicate an unphysical parameter point), • the number of neutral Higgs bosons, nHzero, is greater than 5, • the number of charged Higgs bosons, nHplus, is greater than 1.
The settings for nHzero and nHplus are given as input, either as arguments to the subroutine initialize_ HiggsBounds or on the command line, and are not read from the SLHA file. If HiggsBounds is unable to use an SLHA file (i.e. if one of the situations listed above applies), it might still be possible to run HiggsBounds with one of the other input options. It is nevertheless recommended that the user investigates and understands the reason behind the SLHA failure before proceeding. The supersymmetric spectrum calculator SPheno [39,40] can directly write the HiggsBounds SLHA input blocks to its output SLHA file. 6 For FeynHiggs [41][42][43][44][45]  we provide a stand-alone program, HBSLHAinputblocks fromFH, which creates the necessary SLHA blocks from the FeynHiggs output.

New developments in
HiggsBounds-4

Applying exclusion limits to arbitrary Higgs models
The aim of HiggsBounds is to apply limits derived in Higgs collider searches to models which have not been directly investigated by the experimental analyses. These models can be arbitrary in the sense that they may contain any number of neutral or (singly) charged Higgs bosons, 7 or particles which behave like (elementary) Higgs bosons in Higgs collider searches. Examples of the latter include theories with composite Higgs bosons [46] or dilatons [47]. More specifically, the requirements on the theory in order for the results of HiggsBounds to be reliable are the following: (i) The narrow width approximation must be applicable, such that the predictions for each process can be factorized into Higgs production and decay. (ii) The investigated model should not change the signature of the background processes considerably. Usually, new physics models which show strong deviations from the SM in the background processes of Higgs searches are not considered in the literature, since this would often put them in conflict with SM electroweak precision data [48,49]. Hence, they would most likely not be interesting for HiggsBounds anymore. The presence of such backgrounds would rather correspond to an opportunity for the discovery of physics beyond the SM in other areas. (iii) The investigated model should not significantly change the kinematical distributions of the signal topology X (e.g. the η and p T distributions of the final state particles) from that assumed in the corresponding analysis. For a more detailed discussion of this requirement, see Ref. [12,13].
The above requirements are typically sufficient to ensure the applicability of model-independent exclusion limits, i.e. limits on a cross section of a certain Higgs signal topology, composed of one production and one decay process. If further model assumptions have been made in the experimental analysis, for instance on the CP-properties of the Higgs boson or on the top quark branching ratios, HiggsBounds checks whether the investigated model fulfills them before applying the analysis.
The application of exclusion limits to arbitrary Higgs models becomes less trivial if the experimental analysis combines several Higgs signal topologies under the assumption of a specific model. This is the case for most of the Tevatron and LHC Higgs searches, where a SM Higgs boson is assumed. The exclusion limit is then set on a common signal scale factor for all considered SM Higgs topologies (also called signal strength modifier), usually denoted by μ (sometimes also σ/σ SM is used). For an analysis considering i = 1, . . . , N signal topologies (each consisting of a production mode P and a final state F), the prediction for this quantity can be computed for Higgs boson h k of the investigated model as The channel efficiencies, i , are assumed to be the same for the model and the SM (see requirements (ii) and (iii) above). If these channel efficiencies were published together with the exclusion limit posed by an experimental analysis, the signal strength modifier μ could be computed for a given model without further assumptions. However, these efficiencies have so far been made publicly available only in a very few cases. 8 In HiggsBounds we therefore neglect the channel efficiencies in Eq. (27), leading to an unavoidable model-dependence of the resulting limit, since the calculation of μ via Eq. (27) with all i ≡ 1, is strictly speaking only valid if the model predictions for all signal topologies of the analysis contribute to the total signal rate in (approximately) the same proportions as in the SM.
In order to ensure that an analysis is applied only when this last requirement is fulfilled by the model, Higgs-Bounds performs a SM-likeness test for every Higgs analysis performed under SM assumptions. A test of this kind has been present in all versions of HiggsBounds [12,13]. However, this test was significantly improved in Higgs-Bounds-3.8.0 onwards, as described in Ref. [21], and it is this improved version which we describe here. Neglecting the channel efficiencies, the predicted signal strength modifier μ given in Eq. (27) can be obtained as and are the (SM normalized) channel signal strengths and the SM channel weights, respectively. The requirement that the signal topologies contribute in similar proportions to the total signal rate as in the SM is fulfilled if all channel signal strengths c i are similar to the total signal strength modifier μ (and thus similar to each other). A possible SM-likeness criterion would therefore be to require with δc i = c i − μ and ζ ∼ O (few %), i.e. that the maximal relative deviation of the channel signal strength modifiers from the total signal strength modifier is less than a few percent. In fact, this criterion is very similar to what was used in earlier versions of HiggsBounds. However, this choice was found to be too restrictive in some cases, since it may reject an analysis application which is actually justifiable, leading to overly conservative results. In particular, it is reasonable that channels contributing only a few percent to the total signal rate should be allowed to deviate more from their 8 These efficiencies usually depend on the tested Higgs boson mass. Using a single number for i therefore might appear to be a crude approximation. Nevertheless, for many searches, having access to this information even for one or a few values of the Higgs mass would already provide a better approximation of the full result than in the current situation.
The default setting in HiggsBounds is ζ = 2 %. This is a conservative choice, considering that the uncertainties on the rate predictions for individual channels (even in the SM) are generally larger. With the improved SM-likeness test, the maximal weighted deviation of an individual signal strength modifier from the total signal strength modifier is required to be less than 2 %. Models fulfilling this SM-likeness test for a SM analysis can be safely tested against its exclusion limit.
To illustrate the inclusion of the SM weights ω i in the SMlikeness test criterion, we consider as an example the ATLAS H → γ γ search [50] and test a model with a single Higgs boson with mass m = 125 GeV. We depart from the SM by modifying either the squared effective Higgs coupling ratio to gluons (normalized to the SM), g 2 Hgg , or the coupling to vector bosons, for the production processes (gg → H, VBF, H W, H Z, Htt). In Fig. 4 we show how the total signal strength modifier μ and the c i for the signal topologies are influenced by the modified effective Higgs couplings. Varying g 2 Hgg influences only the gg → H (ggf) cross section. However, due to its large SM weight, ω ggf ≈ 87.7 %, the total signal strength modifier μ follows c(ggf) closely. The failure of the SMlikeness test at g 2 Hgg = 0.835 and 1.225 is therefore eventually caused by the ggf signal topology, although the deviation δc i for the remaining signal topologies VBF, H W , H Z and Htt is much larger here. However, the SM weights of these channels are much smaller. The same effects can be seen when varying g 2 H V V (V = W, Z ). Now, the c i of the VBF, H W , H Z signal topologies are affected by the modified effective coupling, but the total signal strength modifier μ is only slightly influenced due the small weight of these channels. Again, the deviation between μ and c(ggf) eventually causes the SM-likeness test to fail. Due to the inclusion of the SM weights in Eq. (31), subdominant signal topologies are allowed to deviate further from μ.
In comparison with the old SM-likeness test (which was used in HiggsBounds up to version 3.7.0), the new criterion leads to a wider applicability of SM Higgs search results to other Higgs sectors, and thus to a significant improvement of the performance of HiggsBounds. This is shown in  Fig. 2, which has also been produced using the weighted SM-likeness criterion.

Including Higgs mass uncertainties
In several theories with extended Higgs sectors, the Higgs boson masses are not free parameters but can be predicted as a function of the other model parameters up to a certain theoretical accuracy. This is the case, for example, in the MSSM where out of the five physical Higgs states typically only one mass, M A or M H ± , is used as an (on-shell) input parameter. The remaining Higgs masses then become predictions of the  Fig. 1 model, with a theoretical uncertainty that varies within the parameter space and with the sophistication of the theoretical prediction.
We have extended HiggsBounds to be able to take this type of theoretical uncertainty into account when evaluating the Higgs exclusion limits. For theories that have no Higgs mass uncertainties, or where they are negligibly small, this new feature can be left deactivated. In HiggsBounds-4, the Higgs mass uncertainties are taken into account approximately by varying each mass within a user-defined interval. 9 If the tested Higgs boson is unexcluded by the probed limit (in the normal HiggsBounds sense) for any mass in this interval, the tested parameter point of the model is regarded as being unexcluded. This leads to an overall conservative (weaker) result for the exclusion limit when uncertainties are included.
Technically, the number of mass points N considered in the variation can be specified by a variable. The default setting is N = 3 (this corresponds to testing the central mass value, M H , and the two endpoints, M H ± M H , of the specified uncertainty interval). When a sensitive limit varies rapidly with M H , it is advisable to increase N for best results. The mass variation is performed for each Higgs boson independently. In the classic HiggsBounds method this variation is also simultaneous, which leads to a multidimensional computation grid with a complexity growing as O (N n H ), where n H is the number of Higgs bosons with a non-zero mass uncertainty. 10 For the full method, since the limit from each Higgs boson is already considered independently of the others, the complexity remains managable, i.e. O(n H N ). Nevertheless, the user is strongly encouraged to 9 Changes to the relative rates induced by the Higgs mass variation is considered to be "second-order", and are therefore neglected in this procedure. This approximation is valid when the rate predictions vary slowly within the mass uncertainty interval, which sets an upper limit on the reasonable size of the mass uncertainties. 10 To avoid unnecessary calculations, uncertainties smaller than the minimal mass spacing at which the experimental results are available (currently 0.1 GeV for some channels) are not considered.  Fig. 6, which shows the combined exclusion for a SM-like Higgs boson with M H = 0 GeV (solid lines), and similarly for a Higgs boson with SM-like couplings but a theory mass uncertainty of M H = 2 GeV (dashed lines). In this figure, the mass range excluded at 95 % C.L. corresponds to where the limit on σ model /σ SM < 1. Including the mass uncertainties can be seen here as a broadening of the allowed range for the Higgs mass prediction in the model by ±2 GeV around the signal region. It can also be seen that for a given mass point the resulting upper limit on the signal rate is always weaker or equal to the upper limit obtained without theoretical mass uncertainty. Including a theory mass uncertainty in Higgs-Bounds therefore produces overall more conservative limits, which is as expected. 128 GeV is excluded when no theory uncertainty is applied (cf. Fig. 6). The three panels in Fig. 7 show, from left to right, the results when using a mass uncertainty (resulting from the calculation of M h in the model [45]) of M h = 0 GeV, M h = 1 GeV, and M h = 2 GeV. It can be seen that the exclusion at high M A from the limit on the lightest Higgs boson goes down to lower tan β values when M h is small. This illustrates the importance of taking Higgs mass uncertainties into account when interpreting exclusion limits (and compatibility with observed signals, see [15]) in the MSSM and other scenarios for physics beyond the SM.

LEP χ 2 extension
An unfortunate limitation of both the model-independent limits implemented in HiggsBounds, as well as the modeldependent search limits presented by the experiments, is that they are available only at one fixed confidence interval, which is 95 % C.L. for all searches implemented here. The result of an exclusion provided by HiggsBounds based on these searches therefore has a confidence limit of at least 95 % C.L., and in many cases higher. However, the exact level of confidence at which a signal with the properties given to HiggsBounds is either excluded or allowed, is generally unknown.
This has unfortunate consequences for the use of these limits in applications like global fits (see e.g. Refs. [52][53][54] for examples of such fits in the MSSM). There, a model point where the predicted Higgs signal is excluded for exam-ple at 96 % C.L., i.e. with a significance of slightly more than 2 σ , might still be a very good fit if the other properties of the model point in the global fit match the data well. However, the conventional HiggsBounds output only contains information about whether the parameter point is experimentally excluded at at least 95 % C.L. and thus can only be treated as a "hard cut" on the validity of a parameter point.
In order to circumvent this problem, at least for the LEP Higgs searches, the full information on C L s+b and C L s for all Higgs mass combinations in the model-independent LEP searches from [17] have been re-calculated for varying cross sections [55]. These can be written as where σ i,ref is the reference cross section times branching fraction for search i, motivated by the SM Higgs boson or the corresponding cross section for non-SM-Higgs bosons (see [56] for details), and f i is an arbitrary scaling parameter. A logarithmic grid in the scaling parameters f i with 100 points between 10 −3 and 1 is used. Using an interpolation, the actual C L can be calculated for every Higgs production mode at LEP for every physically allowed cross section.
This C L can then be transferred into a quantity whose properties closely follow that of a χ 2 function. This is achieved by assuming that the distribution of −2 ln Q [56] is Gaussian in the asymptotic limit. Transferring the one-sided C L into the two-sided calculation of a χ 2 , the following formula can then be used The resulting χ 2 H can be used as a continuous expression of the agreement between the result of the LEP Higgs boson searches and the model predictions. Note that, in the case of a strong excess in one of the searches, χ 2 H is not only large for models whose predicted cross section times branching fraction is above the observed limit, but also for predictions much smaller than the observed rate in data.

(d)
In cases when the predicted cross section is lower than the minimal (rescaled) value available in the table, the corresponding χ 2 value is set to zero. When the predicted cross section exceeds the tabulated values, no reliable χ 2 value can be calculated, and the value χ 2 = −999 is returned to indicate that a problem has occured. This default behavior can be changed (by setting a flag in usefulbits.f90) to use instead the χ 2 value for the maximal (rescaled) cross section available for that combination of Higgs masses.
An example of the relation between the LEP C L s+b and χ 2 H , also for different values of f , is given in Fig. 8. It can be seen that for C L s+b ≈ 0.5, indicating very good agreement of the signal plus background prediction with the data, fluctuations of χ 2 H around 0 are unavoidable, but numerically irrelevant. In addition, the possibility exists to follow a prescription from [57] to include a mass uncertainty into χ 2 H by folding the full χ 2 distribution with a gaussian G M H with a mass uncertainty M H given by the user, 11 instead of evaluating χ 2 H just at the given M H : 11 Note that this mass uncertainty can be specified completely independently from the uncertainties discussed in Sect. 4.2. This implementation has already been used in global fits of constrained SUSY models [52][53][54]. A non-trivial example of how the LEP χ 2 information can be applied is given in Fig. 9. This figure shows the MSSM low-M H benchmark scenario [24], where the heavier CP-even Higgs boson is interpreted as the LHC signal around M H ∼ 126 GeV. In that case, the lightest Higgs boson, h, is usually below the SM LEP limit and has suppressed couplings to gauge bosons. This is reflected in the figure, where a sizeable χ 2 penalty can be Fig. 9 HiggsBounds results for the LEP χ 2 (colours) in the low-M H scenario of the MSSM [24]. The LEP χ 2 information is shown both on its own (left), and with with the combined LHC exclusion bounds in gray (right). The latest limit from ATLAS charged Higgs searches [58] is not applied here. These results, which are included from Higgs-Bounds-4.1, lead to exclusion over the whole parameter space of this benchmark scenario seen to result in parts of the parameter space, corresponding to regions of low M h (an uncertainty of 2 GeV was used here) where the couplings to gauge bosons is such that the LEP Higgs searches are sensitive to the production of such a state. The sharp edge in the χ 2 distribution in Fig. 9 is obtained at the boundary between two regions of parameter space where the χ 2 contribution comes from the channels e + e − → h Z, h → bb and e + e − → h A → 4b, respectively. Using the LEP χ 2 information together with the Higgs-Bounds exclusion at 95 % C.L. from Tevatron/LHC, Fig. 9 (right) gives the most complete information available from direct Higgs search limits.
Even after the discovery of a SM-like Higgs boson, Higgs boson exclusions still plays, and will continue to play, a vital role in fitting models of physics beyond the SM with an extended Higgs sector. It would therefore be to great advantage if the Tevatron and LHC collaborations could follow the example of the LEP Higgs WG and provide exclusion limits for varying values of f σ ref in addition to the results that are presented at 95 % C.L..

User operating instructions
In this section we describe in detail the two main methods to use HiggsBounds-4: The command line version and the library of subroutines. There is also an online version that provides quick access to all the functionality of Higgs-Bounds, without the need to install the code.

Installation
The HiggsBounds source code, the online version, and the documentation can all be obtained at the URL http:// higgsbounds.hepforge.org The HiggsBounds-4 code is mostly written in Fortran 90 but includes also a few Fortran 2003 features. It has been tested with a variety of Fortran compilers, including the free GNU compiler 12 (gfortran) which accompanies most Linux distributions.
Before compiling the HiggsBounds code, the user should first make changes to the configure script to appropriately reflect the compiler and path settings on the user's system. The code can then be compiled by running ./configure make which creates the main HiggsBounds executable and the library of subroutines, libHB.a. Any program for which the HiggsBounds subroutines should be used can be compiled and linked to the library by adding −L < HBpath > −lHB to the command line, for example, where <HBpath> is the location of the HiggsBounds library. The HiggsBounds subroutines make use of the Fortran file handles 10, 11, 44, 45 and 87, which means that users should avoid these file handles in programs calling the subroutines.
The default behavior of HiggsBounds-4 is to use the full (new) method to generate combined exclusion. To set the classic method as the default, the user can modify the flag run_HB_classic in the file usefulbits.f90 before compiling HiggsBounds. When running the subroutine version of HiggsBounds, it is also possible to access the results from both methods without changing the default behavior, see below.
The library of subroutines and the command-line version share a common set of features, which we will describe first. We will then give the proper operating instructions to use each of these HiggsBounds formats individually.

Common features: Input
Regardless of the operation mode, HiggsBounds requires five basic types of user input: (i) the number of neutral Higgs bosons in the model under study (nHzero) (ii) the number of (singly, positively) charged Higgs bosons in the model (nHplus) (iii) the set of experimental analyses which should be considered (whichanalyses) (iv) the desired input format for the theoretical predictions (whichinput) (v) the theoretical predictions for the Higgs sector of the model (given as arrays) The variables nHzero and nHplus are currently both limited to the range 0-9, but if necessary this could easily be extended in the future. The possible values for the choice of experimental analyses (whichanalyses) are described in Table 4.
HiggsBounds expects the theoretical input to be provided in one of three formats labelled by the variable whichinput. These formats are described in detail in Sect. 3, and their required input is briefly summarized in Table 1. In Appendix 8 ( Tables 11,12,13,14) we assign names and list the full contents of all the possible input arrays for the theory predictions. These names will be used below to describe the input requirements of each version of HiggsBounds individually.

Common features: Output
HiggsBounds provides the user with four types of basic output: with the highest statistical sensitivity (chan) (iii) the number of Higgs bosons that contributed to the theoretical rate for the corresponding process (ncombined) (iv) the ratio k 0 = Q model /Q obs for the process with highest statistical sensitivity (obsratio).
As discussed in Sect. 2, the extended HiggsBounds algorithm now offers similar quantities to be calculated individually for each Higgs boson in the model. When making use of the full method, the corresponding output quantities are promoted to arrays of length n + 1, where n is the total number of (neutral and charged) Higgs bosons in the model. The combined result (contained in element 0 of these arrays) from this extended test can also be used in analogy to the result of HiggsBounds classic. When several Higgs bosons exclude the same point through different searches, the values for chan, obsratio, and ncombined in the combined result refers to the channel giving the strongest exclusion. Table 5 shows how to interpret the possible values of HBresult and obsratio (one entry in the case of arrays), which are complementary. When using either the library of subroutines or the command-line version, the keys associating the reference numbers (as given by chan) with the analysis applications is written in human-readable format in the file Key.dat. In the online version, this information appears directly on the screen. When the SLHA option is used for input, the HiggsBounds results can be added to SLHA files in the form of a new block, called HiggsBoundsResults. An example of this block is shown in Table 6.

Library of subroutines
In this section we list all the user subroutines available through the HiggsBounds library.

initialize_HiggsBounds_int(int nHzero, int nHplus, int whichanalyses)
This is an alternative version of initialize_ HiggsBounds, which takes an integer for the third argument instead of a string constant. This code specifies the set of experimental data that is used by HiggsBounds according to the first column of Table 4.
This subroutine sets the model input for the neutral Higgs sector using the effective couplings (whichinput=effC), as defined in Sect. 3.3. Using this method also excludes the use of either parton-level or hadron-level input. The meaning of the input arrays (all of length n = nHzero > 0) is summarized in Appendix 8, Table 11. If any of the effective couplings are deemed to be irrelevant, the corresponding array may be filled with zeros. However, this needs to be exercised with some caution for quantities relevant in a SM-like Higgs search (as most of the limits reported from the LHC are). It is possible that setting certain couplings artificially to zero could lead to the model failing the SM-likeness test, cf. Sect. 4.1.
This routine is used to set the input for the neutral Higgs sector using parton-level cross sections (whichinput=part), as defined in Sect. 3.2. Using this method excludes the simultaneous use of effective couplings or hadron-level input. The meaning of the input arrays (of length n = nHzero > 0) are summarized in more detail in Appendix 8, Table 12 (cross sections) and Table 13 (branching ratios). As for the effective coupling case, quantities which are not required by any channel that has a competitive sensitivity can be set to zero to simplify the input (and the same caveats about searches for SM-like Higgs bosons apply).
This subroutine sets the input for the neutral Higgs sector using hadron-level cross sections (whichinput=hadr), as defined in Sect. 3.1. Using this method excludes the use of effective couplings or parton-level input. The names for the input arrays (of length n = nHzero > 0) are described in Appendix 8, Table 12 (cross sections) and Table 13 (branching ratios). Similarly to the effective coupling case, quantities which are not required by any channel that has a competitive sensitivity can be set to zero to simplify the input (and the same caveats about searches for SM-like Higgs bosons apply).
The subroutine HiggsBounds_charged_input gives the charged Higgs sector input to HiggsBounds. The use of this subroutine is only required if k = nHplus is nonzero (recall that nHplus is set in subroutine initialize_ HiggsBounds). Currently, only results from searches for light charged Higgs bosons (M H ± < m t ) are available. Once results from heavy charged Higgs searches are presented, this interface will be extended with input of the necessary cross sections. The names used for the input arrays are described in Appendix 8, Table 13.

HiggsBounds_input_SLHA(char(:) SLHAfilename)
This subroutine can be used for supersymmetric theories as an alternative to the other routines which provide model input to HiggsBounds. It is called with a string-type variable, SLHAfilename, which gives the name of an SLHA file (full path should be included if not in the current working directory). The model predictions are then read in from this file, which should contain the two HiggsBounds-specific blocks as described in Sect. 3.4. Furthermore, it will set the mass uncertainties of the neutral and charged Higgs bosons according to the values given in the SLHA block DMASS (if available).

HiggsBounds_set_mass_uncertainties (double(n) dMhneut, double(k) dMhch)
This subroutine allows the user to specify theory mass uncertainties for the neutral and charged Higgs bosons of the model. The implementation and use of these uncertainties when the exclusion limits are evaluated is discussed in detail in Sect. 4.2. The default is for all the uncertainties to be zero. The treatment of mass uncertainties in the limit setting is invoked automatically by setting any of them to a non-zero value. The routine takes two arrays as arguments: dMhneut(n) (of length n = 1 . . .nHzero), which specifies the (absolute) uncertainties for the neutral Higgs boson masses in GeV, and dMhch(k) (k = 1 . . .nHplus) which does the same for the charged Higgs bosons. If either nHzero = 0 or nHplus = 0, the corresponding uncertainty array will not be used (and can therefore be set to arbitrary values).
run_HiggsBounds (double HBresult, int chan, double obsratio, int ncombined) After initializating and setting the model input using one of the methods discussed above, this subroutine is called to perform the main part of the HiggsBounds calculations.
The results from the run is given as output. The combined result, HBresult, is reported according to the description in Table 5. The channel with the highest exclusion power is identified by its code, chan (the channel codes are written to the file Key.dat), and the corresponding ratio of the model prediction to the observed limit in this channel is given by obsratio. Finally, the number of Higgs bosons combined in this prediction is ncombined. The default behavior of this subroutine (which can be controlled by setting a flag in usefulbits.f90) is to use the full exclusion method of HiggsBounds, rather than the classic method employed in previous versions. This subroutine produces the results for a single Higgs boson, which should be indexed by h. The indexing is such that the neutral Higgs bosons correspond to h = 1 . . . nHzero, followed by the charged Higgs bosons of the model, h = nHzero + 1 . . . nHzero + nHplus. To get the results for more than one individual Higgs boson, it is recommended to instead use the subroutine run_HiggsBounds_full for better performance.

HiggsBounds_SLHA_output()
When using the SLHA input, the subroutine HiggsBounds _SLHA_output can be called after using (any of the different) run_HiggsBounds routines in order to write the block HiggsBoundsResults to the SLHA file. The results are written in terms of the combined exclusion, see Table 6 for an example.

finish_HiggsBounds()
The subroutine finish_HiggsBounds should be called once at the end of the program, after all other Higgs-Bounds subroutines. This deallocates the allocatable arrays used within HiggsBounds.

Command-line version
When using HiggsBounds from the command line, the run options are specified in the program call, which should be of the form The variable <prefix> is a string which is added to the front of input and output file names. It may include directory names or other information identifying the run files. If whichinput = SLHA, < prefix > should contain the full name of the SLHA file to use, including the path if it is not in the current working directory. When running HiggsBounds from the command line, the program behaviour (full/classic) is determined by a flag specified in the file usefulbits.f90 (the same as for the subroutine run_HiggsBounds). The default setting is that the full method is used.

Input file format
The arrays containing the theoretical model predictions are read in from text files, with each value given in a separate column (separated by whitespace). The contents of each input file is described in Tables 7 and 8. Note that all these files will not be necessary at the same time. This will be specified below. Each row in the input files starts with a line number, k, which identifies predictions belonging to the same parameter point in different files. The input files must not contain any comments or blank lines. Care should be taken with the order of the array elements in the files. The elements of a onedimensional array, e.g. Mh for nHzero = 3, is given in the order The correspondence between the array elements and the physical input quantities is clarified in Appendix B. Not all of the elements of the two dimensional arrays are required. From the array g2hjhiZ only the lower left triangle (including the diagonal) is required (and similarly for lepCS_hjhi_ratio below), since they are symmetric matrices. From the general matrix ⎛ ⎝ g2hjhiZ(1, 1) g2hjhiZ(1, 2) g2hjhiZ(1, 3) g2hjhiZ(2, 1) g2hjhiZ(2, 2) g2hjhiZ(2, 3) g2hjhiZ(3, 1) g2hjhiZ (3,2) the required elements should be written in the input file using the order g2hjhiZ(1,1), g2hjhiZ(2,1), g2hjhiZ(2,2), g2hjhiZ(3,1), g2hjhiZ(3,2), g2hjhiZ (3,3). The file MHall_uncertainties.dat is optional. If it is provided, HiggsBounds will automatically include the theoretical Higgs mass uncertainties given in the file. This has been described in more detail in Sect. 4.2. The file additional.dat is also listed as optional. If this file is included, it can have any number of columns greater than 1 (as for the previous files, the first entry on each line should still be the line number). This input is particularly useful to keep track of variables which are not required by Higgs-Bounds, but which are helpful when plotting the results from a parameter scan. For example, in the case of the MSSM, additional.dat could be used to store the values of tan β.
As in the subroutine version, the command line version of HiggsBounds expects a subset of the total list of possible input arrays, which depends on the chosen setting of whichinput. The maximal list of files used for each value of whichinput is given in Table 9. Furthermore, some of the arrays will not be relevant for some of the choices for whichanalyses. The command line version of Higgs-Bounds will consider the list of input files appropriate to Table 8 File names and data format for the contents of HiggsBounds input files (part II). The right column shows the order of the input data arrays within one row of the input file (k is the line number). For the order of elements within the arrays, see the text for details. Note that several arrays appear in two different input files. These files are never required simultaneously in one run of HiggsBounds the settings of whichinput and whichanalyses, and then attempt to read only those input files where at least one of the arrays contained in the file will be used. For supersymmetric models, one possible way of generating HiggsBounds input files is to use the model building tool SARAH [59][60][61] in conjunction with the spectrum generator SPheno [39,40]. These codes can directly write out the HiggsBounds input files required for the effective coupling approximation.

Output file format
When the command line version of HiggsBounds is used with whichinput=hadr, part or effC, the output is written to the file < prefix > HiggsBounds_ results.dat. A sample of this output is shown in Fig.  10. The key to the process numbering is written to the file < prefix >Key.dat.
When the command-line version of HiggsBounds is used with whichinput=SLHA, the HiggsBounds results are added to the SLHA file in the form of the block Table 9 List of possible input files for each setting of whichinput. Optional files are marked with (*). The files required can also depend on the setting of whichanalyses, see HiggsBoundsResults, see Table 6. It should be noted that it is not efficient to use the command-line version of HiggsBounds with SLHA input for large parameter scans, since the experimental data tables must be read in again for each SLHA file. If this is a concern, a better option is to use the HiggsBounds subroutines to create a program which can be called from the command line. An example program, HBwithSLHA, demonstrating this is included.

Example input files
The HiggsBounds package includes a full set of sample input files for the case n H = 3, n H + = 1, contained in the folder example_data. Each filename is prefixed with HB_randomtest50points_. To run the command-line version of HiggsBounds with these files as input, use, for example, ./HiggsBounds LandH effC 3 1 example_data/

HB_randomtest50points_
where the values of whichanalyses and whichinput can be varied as desired. The setting nHplus = 0 can also be used if the user does not wish to test the charged Higgs sector, ./HiggsBounds LandH effC 3 0 example_data/ HB_randomtest50points_

Example programs
We provide a number of example programs which demonstrate the different features of the HiggsBounds subroutines. These are available in the subfolder /example_ programs/ of the main installation directory. After the HiggsBounds library (libHB.a) has been compiled (using ./configure; make as described previously), each of the examples can be compiled and with the command make <program name> More generally, any program linking with the Higgs-Bounds libraries can be compiled (assuming gfortran is used) as follows Table 10 List of input files relevant to each setting of whichanalyses (marked by 'y'). The required files also depend on the settings of whichinput, nHzero and nHplus; see Table 9 [41][42][43][44][45]. The model parameters should be specified in the source file (see the code for details). The results are written directly to the screen.

• HBwithCPsuperH
This example is similar to the HBwithFH example above, but uses the spectrum generator CPSuperH [62][63][64] for the theory predictions instead of FeynHiggs. As above, the model parameters should be specified directly in the source file and the results are written directly to the screen.

• HBwithFH_dm
This is an updated version of the HBwithFH example, which demonstrates the use of several new features in HiggsBounds-4. It makes it use of FeynHiggs- When using this example, the user can provide input in the SLHA format with one or more input files. The program run settings are fixed as | <whichanalyses>= LandH|, | <nHzero>= 3| and | <nHplus>= 1|, but this can be changed in the code. The set of SLHA files to be used as input should be named | <stem> .i| where |i| = 1 . . . |npoints|. The program can be called from the command line as: ./example_programs/HBwithSLHA <npoints><stem> The block |HiggsBoundsResults| will be added to each SLHA file. In addition, the HiggsBounds results (|i|, |HBresult|, |chan|, |obsratio|, |ncombined|) are collected in a summary output file called | < stem > −fromHB|. As in the command-line version, the HiggsBounds results are obtained using either the full (default) or classic method, following the setting of the corresponding flag in usefulbits.f90.
In addition to the programs listed here, there are two more example codes specifically for the LEP χ 2 extension. These are discussed in the next section.

Installing and using the LEP χ 2 extension
The usage of the LEP χ 2 information is restricted to the subroutine version of HiggsBounds. In order to enable this feature, the user first needs to download a separate package containing the binary files with the relevant experimental information from the URL http://higgsbounds.hepforge.org.
These files are contained in the tarball csboutput_ trans_binary.tar.gz, which should be extracted to a user-defined directory < clsbtablesdir > (not exceeding 80 characters), such that the following file structure is obtained:

binary
A convenient choice for <clsbtablesdir> might be the HiggsBounds main directory. In the next step, <clsbtablesdir> has to be specified in the script configure-with-chisq, in addition to the usual com-piler settings, cf. Sect. 5.1. Then, the HiggsBounds library can be built with: ./configure − with − chisq make After a successful compilation, new subroutines for the LEP χ 2 extension are available. These are described in the following.

initialize_HiggsBounds_chisqtables()
This subroutine initializes the new arrays and tables needed for the LEP χ 2 extension. It reads in all the relevant experimental information from the binary files installed in <clsbtablesdir>.
HB_calc_stats (double theory_uncertainty, double chisq_withouttheory, double chisq_withtheory, int channel) This routine is run to calculate the LEP χ 2 value. The user can specify a theoretical mass uncertainty (in GeV), theory_uncertainty. Note that this value is only used here, and not in "standard" mass uncertainties for the limits (which can be different). The resulting χ 2 value is reported both including (chisq_withtheory) and without (chisq_withouttheory) this Higgs mass uncertainty. The channel code for the experimental analysis from which the χ 2 value is derived is also given. This subroutine requires a preceeding call to the subroutine run_HiggsBounds_classic, in order to determine the most sensitive analyses for the model. Therefore, the usage of the LEP χ 2 extension always requires a simultaneous run of the standard HiggsBounds program.

finish_HiggsBounds_chisqtables( )
This deallocates the new arrays and tables, and should be called at the end of a run.

Usage
The typical sequence of subroutine calls when using the LEP χ 2 extension is the following: Note that the LEP χ 2 functionality requires a classic HiggsBounds run to determine the most sensitive channel. For a consistent combination of LEP χ 2 extension with Tevatron and LHC limits we recommend to initialize the LEP χ 2 functionality with the option whichanalyses="onlyL" and performing a separate HiggsBounds run to consider the hadronic collider limits, i.e. whichanalyses= "onlyH". Two example programs are provided for the LEP χ 2 extension. They are called HBchisq and HBchisqwithSLHA, respectively, and are both contained in the directory example_programs of the main Higgs-Bounds directory. After setting up HiggsBounds to use the LEP χ 2 extensions, these examples can be compiled with make HBchisq The first example, HBchisq, simply scans over the SM Higgs boson mass within the range M H ∈ [100, 120] GeV and evaluates the LEP exclusion χ 2 value. The second program, HBchisqwithSLHA, runs HiggsBounds on a set of n SLHA files, named <SLHA-filename>. i with i = 1 . . . n. It is called as /HBchisqwithSLHA < number of files > < SLHA − filename > The results, including the LEP χ 2 values, are printed for all parameter points to the new file < SLHA − filename > − fromHB.

Conclusions
We have presented HiggsBounds-4, an extension of the HiggsBounds program which can be used to study exclusion bounds on arbitrary Higgs sectors using experimental results from LEP, the Tevatron and the LHC. It includes the latest LHC results presented in 2013, many of which are based on the full 8 TeV dataset.
We briefly reviewed the options for user input, including a new option (in the case of a supersymmetric Higgs sector) which allows an SLHA file to be used as input. Several improvements and updates of the code have been presented. This includes in particular an improved SM-likeness test that now takes into account the relative weight of the contributing channels, and in this way substantially enlarges the parameter space in which SM analyses can be applied. We have included the option of a theoretical Higgs mass uncertainty, which can be relevant, e.g., in the MSSM. Taking the theory uncertainties into account conservatively broadens the range of non-excluded Higgs mass values. Concerning the LEP limits, we include an option to obtain the full χ 2 informa-tion, i.e. not "only" a hard 95 % C.L. cut. This is particularly useful for fits in the Higgs sector.
In view of the discovery of a Higgs signal at the LHC at ∼ 125.5 GeV we have included the option to test every Higgs boson in the model under consideration individually. In this way we slightly deviate from the pure 95 % C.L. exclusion limit, but we ensure that models do not falsely pass the HiggsBounds test because the spectrum contains one (SM-like) Higgs boson at a mass of ∼ 125.5 GeV.
HiggsBounds can now readily be used together with its new sister code, HiggsSignals [15]. HiggsSignals performs a χ 2 evaluation of the compatibility between the predictions of arbitrary Higgs sectors to measured signal rates. This includes in particular the possibility to test the model predictions against the observed signal at ∼ 125.5 GeV, but also future, hypothetical, signals of extended Higgs sectors. A combined analysis using both codes exploits all the public information on the Higgs signal and the Higgs exclusion bounds obtained at LEP, the Tevatron and the LHC. gram under grant MultiDark CSD2009-00064. The work of T.S. was partially funded by the Bonn-Cologne-Graduate-School. O.S. is supported by the Swedish Research Council (VR) through the Oskar Klein centre.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited. Funded by SCOAP 3 / License Version CC BY 4.0.

Analyses included in HiggsBounds-4
The intention is to keep HiggsBounds continuously updated with the latest experimental results as they become available. Older results, which have been surpassed in sensitivity by newer analyses, are removed. After compiling HiggsBounds, the user can run the command ./AllAnalyses to print an up-to-date list of the implemented experimental results to the screen. Data from the following experimental analyses is included in HiggsBounds-4:  Fig. 2  In addition to the analyses listed above, HiggsBounds-4.1 contains the results from the following experimental analyses: ATLAS Collaboration [58, 106] CMS Collaboration [107,108] In particular the updated ATLAS results from light charged Higgs searches [58] are interesting for constraining the MSSM (and other models with multiple Higgs doublets) in the region M H ± < 160 GeV. In Fig. 11 we show an updated version of the results from charged Higgs exclusion in the M mod+ h scenario presented in Fig. 2(c). The new limit excludes small values of M A for all tan β. The same limit also excludes the whole parameter space of the MSSM low-M H scenario displayed in Fig. 9.

Appendix: Input arrays
This appendix contains tables which define the names of the input arrays the user must give as arguments to the various subroutines and/or data files used to specify the model predictions. The physics definitions for the different input quantities (given in the right-hand columns of these tables) follow Sect. 3 (Tables 11,12,13,14).