HiggsBounds-5: Testing Higgs Sectors in the LHC 13 TeV Era

We describe recent developments of the public computer code HiggsBounds. In particular, these include the incorporation of LHC Higgs search results from Run 2 at a center-of-mass energy of 13 TeV, and an updated and extended framework for the theoretical input that accounts for improved Higgs cross section and branching ratio predictions and new search channels. We furthermore discuss an improved method used in HiggsBounds to approximately reconstruct the exclusion likelihood for LHC searches for non-standard Higgs bosons decaying to $\tau\tau$ final states. We describe in detail the new and updated functionalities of the new version HiggsBounds-5.

We describe recent developments of the public computer code HiggsBounds. In particular, these include the incorporation of LHC Higgs search results from Run 2 at a center-of-mass energy of 13 TeV, and an updated and extended framework for the theoretical input that accounts for improved Higgs cross section and branching ratio predictions and new search channels. We furthermore discuss an improved method used in HiggsBounds to approximately reconstruct the exclusion likelihood for LHC searches for non-standard Higgs bosons decaying to τ τ final states. We describe in detail the new and updated functionalities of the new version HiggsBounds-5.

Introduction
After the discovery of a Higgs boson with a mass around 125 GeV [1,2] the searches for new scalars have intensified and expanded into more and more search signatures. Evidence for the existence of additional neutral and/or electrically charged Higgs bosons would be an unambiguous sign of physics beyond the Standard Model (BSM), where the SM scalar sector is extended by new scalar field(s). These fields could be singlets, doublets or even higher representations of the electroweak gauge group SU(2) L . Well-known examples of such extensions are the real or complex scalar singlet extension [3][4][5][6], featuring one or two additional neutral scalar bosons, respectively, and the Two-Higgs-Doublet Model (2HDM) [7][8][9], which contains three neutral Higgs bosons -typically denoted h, H and A in the CP-conserving case -and a pair of charged Higgs bosons, H ± . The Higgs sector of the Minimal Supersymmetric Standard Model (MSSM), for instance, is at the tree-level a specific version of a 2HDM [10][11][12]. Examples of BSM Models containing both additional scalar doublets and singlets are the Next-to-2HDM [13,14] and the Next-to-MSSM [15,16], while higher representations of scalars are considered e.g. in the Georgi-Machacek model [17].
Since -so far -no additional Higgs bosons have been discovered at the LHC all of these searches have resulted in exclusion limits that constrain the possible parameter space of BSM theories with extended Higgs sectors. Due to the large number of results in many different search channels, the task of testing BSM model predictions against the assembled results from Higgs searches warrants dedicated tools. The tool HiggsBounds [18][19][20][21] has been developed to perform such a check against all available Higgs searches from LEP, the Tevatron, and the LHC. This paper presents the upgrades to the code in HiggsBounds-5 compared to the previous version HiggsBounds-4 described in Ref. [20] and discusses the most important new features.
The general approach of HiggsBounds remains unchanged compared to HiggsBounds-4 and we refer to Ref. [20] for the details. For each Higgs boson of the investigated model, based on the model predictions input by the user, HiggsBounds selects the most sensitive limit by comparing the model predictions to the expected limits of all analyses. It then checks the observed limit of this selected analysis against the model predictions to obtain a bound for each Higgs boson. The model parameter point is considered allowed if none of the Higgs bosons are excluded by the corresponding selected analysis. All of the individual limits implemented in HiggsBounds are exclusion limits at 95% confidence level (C.L.). This procedure ensures that the overall combined limit is still approximately at the 95 % C.L. 1 Likelihood information has been made available by the experimental collaborations for several analyses at LEP and at the LHC. HiggsBounds utilizes this detailed input to reconstruct the corresponding 95 % C.L.
limit in a nearly model-independent fashion or, optionally, return the corresponding χ 2 that can be used in a model fit.
In light of the increasing number of experimental search channels with improved sensitivity one of the most important aspects of HiggsBounds is to provide an input framework that works for a large class of BSM models and can incorporate all of the required model predictions. Most of the major changes in HiggsBounds-5 relate to improvements in the input frameworksuch as allowing additional cross sections and branching ratios to be set as input or providing precise, model-independent parametrizations of required input quantities. The extended input framework is also used by the code HiggsSignals-2 [22][23][24] which tests BSM models against the Higgs rate measurements at the LHC (and the Tevatron).
In Section 2 we describe the extended HiggsBounds input framework in detail, discussing the available input schemes and the production and decay channels supported by HiggsBounds-5. This includes the possibility of providing input beyond the narrow width approximation as detailed in Section 2.3. In Section 3 we review the experimental input used by HiggsBounds and consider possible improvements to the input presentation that could make the experimental results more readily useable. Section 4 discusses functional changes in HiggsBounds-5. This includes parametrizations of V H and charged Higgs production cross sections in the effective coupling approximation, as well as a description of improvements in the derivation of exclusion limits from the available likelihood information compared to the method first discussed in Ref. [21]. We give an overview of the technical changes relevant to users of the code in Section 5 and conclude in Section 6.

Theoretical Input
The input for HiggsBounds consists of the phenomenologically relevant physical quantities of the Higgs sector, i.e. the number of neutral and charged Higgs bosons that should be considered, their masses, total decay widths, production and decay rates. By relying only on these physical quantities (in contrast to model-specific parameters), the code maintains a flexible input framework with rather minimal model assumptions.
HiggsBounds-5 supports three types of input specified in the variable whichinput when initialising the code. These are the hadronic cross section input (whichinput = 'hadr'), the effective coupling input (whichinput = 'effC'), and SLHA (SUSY Les Houches Accord [25,26]) input (whichinput = 'SLHA'). The partonic input mode (whichinput = 'part') present in previous versions of HiggsBounds has been removed and is no longer supported. We describe the available methods of providing input to HiggsBounds in more detail in Section 5.
For whichinput = 'hadr' the inclusive hadronic 2 production cross sections have to be provided to the code. Cross sections for various colliders and center-of-mass (CM) energies are required -namely LEP, Tevatron, and the LHC at 7, 8 and 13 TeV. If the considered production mode also exists in the SM, the input cross section is normalized to the corresponding SM prediction for the same Higgs mass, and otherwise specified in picobarn (pb). In effC and SLHA input the hadronic cross sections are calculated internally from the provided effective couplings whenever possible. The branching ratios (BRs) for all Higgs boson decays are also required as input. In effC mode the BRs for decay modes to SM particles are per default approximated from the provided effective couplings, and only the BRs for Higgs decay modes that are not present for a SM Higgs boson have to be specified explicitly. In contrast, in SLHA input, the BRs for all decay modes are directly taken from the SLHA DECAY blocks.

Neutral Higgs Bosons
The quantities needed to describe the production and decay rates of neutral Higgs bosons are listed in Tables 1 to 3. The hadronic cross sections in Table 1 have been extended by separate input for gluon fusion and bb associated Higgs production, whereas HiggsBounds-4 only required the sum of the two processes (denoted as single Higgs production). This is particularly relevant for exclusion limits from analyses that have specific requirements on the b-jet multiplicity in the event, as is e.g. the case in searches for heavy BSM Higgs bosons decaying to τ + τ − (see also Section 4.3). Furthermore, the cross sections for processes of Higgs production in association with a single top quark have been added. We distinguish between the t-channel process qb → tqh j (specified in the 5-flavor scheme) and the s-channel process qq → tbh j (see Ref. [27] for a comprehensive discussion and for the NLO QCD predictions in the SM). Similarly, separate cross sections for gluon-and quark-initiated Z-boson associated Higgs production have been added to the HiggsBounds input. These involve different Higgs couplings (see Section 4.1) and are partly separated in the simplified template cross section (STXS) measurements that are newly included in HiggsSignals-2 (see Ref. [24]). Since HiggsBounds handles the input for HiggsSignals, this also means that this subchannel information is available, such that differential information can be incorporated if it becomes available in a framework similar to the STXS. Finally, the non-resonant double Higgs production cross section has been added as input. Note that it is not normalized to the SM prediction but should instead be given in pb.
The cross section input for LEP is unchanged with respect to HiggsBounds-4. For completeness, these quantities are listed in Table 2.  hj → ss hj → gg hj → e ± µ ∓ (NEW) hj → e ± τ ∓ (NEW) tion for all these quantities. In case this approximation is employed, the effective couplings listed in Table 4 have to be provided. With respect to HiggsBounds-4 we have changed the entire input from squared effective couplings (or scale factors) to the sign-sensitive single effective couplings (or scale factors). This allows us to take into account interference effects e.g. in the prediction for the h j Z production cross section. Furthermore, we removed the effective (squared) h j ggZ coupling present in earlier versions. Instead, the gg → h j Z contribution is derived from the h j tt and h j bb effective couplings. Note that the loop-induced h j gg, h j γγ and h j γZ couplings are still free input quantities not derived from the other coupling parameters.
The scalar and pseudoscalar components of the Higgs couplings to a generic fermion pair ff are defined through  where g s and g p are the real-valued scalar and pseudoscalar coupling constants. As effective couplings, they are both normalized to the SM value of g s for the corresponding fermion given by with electric charge e, the weak mixing angle θ w , fermion mass m f and the W -boson mass M W . The couplings to W and Z bosons are normalized to the corresponding SM tree-level couplings where M Z is the mass of the Z-boson. The loop-induced effective couplings to γγ and Zγ are best defined through the partial decay widths normalized to the SM-value for the same Higgs mass. This can also be used for the gg effective coupling, however it is a better approximation in this case to use the normalized gluon fusion production cross section. Either way, these yield the squared effective coupling whose square-root is the input expected by HiggsBounds-5. The sign of these loop-induced couplings does not enter any observables, so the positive square root can be used without loss of generality.
Finally, the h j h i Z coupling does not have a SM equivalent that could be used for normalization. It is instead normalized to

Charged Higgs Bosons
The HiggsBounds input framework has been broadly extended in the charged Higgs sector. We list all relevant charged Higgs sector quantities in Table 5. HiggsBounds-5 supports direct charged Higgs boson production at hadron colliders, including H ± j production in association with a top or charm quark and a bottom quark as well as flavor-suppressed production in association with lighter quark jets. We also include charged Higgs production in association with a vector boson or a neutral Higgs boson, as well as charged Higgs boson production in vector boson fusion and charged Higgs pair production. Note that all hadronic cross sections are directly given in pb, and not specified as normalized quantities. All input cross sections are required to be summed over the two possible charges. Note that at present there is no effective coupling input for the charged Higgs bosons.
For light charged Higgs bosons with mass below the top quark mass, the most important search channel is top quark pair production with successive decay of one top quark to a charged Higgs boson and a bottom quark. HiggsBounds thus also requires the branching fractions for t → H + j b and t → W + b -where the latter is needed to check for model assumptions -as input. The charged Higgs branching fractions have been extended by the decays to top and bottom quarks, W and Z bosons, as well as neutral Higgs and W bosons.
Note that, thus far 3 , LHC searches have only considered pp → H ± tb and H ± production in vector boson fusion as direct production channels. The remaining cross sections listed in the upper section in Table 5 are therefore only placeholders at the moment, and setting them to zero in a HiggsBounds run will not affect the results until relevant experimental results become available.

Input Beyond the Narrow Width Approximation
All of the input schemes described above rely on the narrow width approximation (NWA) to construct the signal rates in specific collider channels from the provided cross sections and branching ratios. As such, in the NWA the channel rate r p,d is given by where p denotes the production and d the decay mode of the channel. In cases where the NWA is not applicable, the hadron collider channel rates r p,d for neutral Higgs bosons can  be specified directly by the user, which then replace the corresponding values obtained from the NWA. In this way, individual channel rates can be set while keeping the remaining HiggsBounds input unchanged. For instance, this is relevant if one of the neutral Higgs bosons of the model has a very large width, while the narrow width approximation is applicable for the remaining particles. Non-trivial modifications through signal-signal or signal-background interference can also be accounted for by explicitly setting the channel rate. For instance, destructive signal-signal interference of two heavy Higgs bosons can appear in the MSSM with CP-violation [28], leading to sizable differences in the exclusion obtained from BSM Higgsto-τ + τ − searches, as compared to the naive, incoherent combination of the individual Higgs boson signal rates (see e.g. the discussion in [29]). 4 Note that there is currently no way to specify channel rates for charged Higgs processes.
Channel rates can be input using the Fortran subroutine interface (using the subroutine HiggsBounds_neutral_input_hadr_channelrates or, for a single process, the subroutine HiggsBounds_neutral_input_hadr_channelrates_single), see the online documentation for details. These subroutines expect the channel rates to be normalized as is the corresponding production cross section for a SM-like Higgs boson of the same mass. The HBwithchannelrates example programs illustrates the use of explicitly set channel rates.
Besides accounting for width effects in the theoretical input for the channel rates, the experimental limits in various search channels are often provided as a function of the total decay width. In that case, this width-dependence of the limit is fully implemented in HiggsBounds-5 and accounted for in the model testing, irrespectively of whether the NWA was employed or the channel rates have been set directly by the user.

Experimental Input
In this section we describe the experimental results that are used by HiggsBounds, and we address possible limitations of the application of Higgs search limits to a model. In this context we discuss how search limits (at fixed confidence level (C.L.) or as a likelihood) should be presented, in particular how they should be parametrized and what information is required in order to apply an experimental limit to (nearly) arbitrary Higgs models. Finally, we make suggestions for possible future refinements in the presentation of experimental limits.

Experimental data in HiggsBounds and Limitations of Applicability
HiggsBounds currently incorporates results from LEP [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44], the Tevatron , and the ATLAS [2, and CMS  experiments at the LHC. A detailed list of the implemented analyses is returned by the AllAnalyses executable (see Section 5.2). An up-to-date version of this list, together with a bibliography of all implemented results is also available on the webpage, and the InspireHEP cite keys of the analyses are included in the HiggsBounds output (in the Keys.dat file). We expect all users of HiggsBounds to cite the relevant experimental analyses.
The application of the experimental exclusion limits to a model parameter point is described in detail in Ref. [20]. The basic procedure is as follows: In the first step, based on the expected exclusion limit (at 95 % C.L.), HiggsBounds selects the most sensitive analysis for each Higgs boson of the model. In the second step, the model predictions for each Higgs boson are compared with the observed limit from the particular experimental search that is most sensitive to it. If the predicted signal rate exceeds the observed limit for any of the Higgs bosons, the model parameter point is regarded as excluded (at 95 % C.L.). The validity of this test depends, in short, on the following basic assumptions (see Ref. [20] for details): • the narrow width approximation is valid, i.e. the signal rate can be approximated by the product of the Higgs boson production cross section and branching ratio 5 ; • background processes in the experimental analyses are not altered significantly by the signal (new physics) model; • the kinematics of the signal processes are not altered significantly with respect to the signal hypothesis employed in the experimental analysis (typically, a scalar or pseudoscalar boson with renormalizable couplings).
Furthermore, there are experimental analyses that combine different Higgs boson search channels. Such combined limits may require an applicability test, i.e. a check whether the parameter point fulfills the assumptions on which the combination is based. If this applicability test fails, the corresponding search limit is not considered in the HiggsBounds test of the parameter point. Examples of such analyses are searches for a SM-like Higgs boson or an invisibly decaying Higgs boson, where various production modes are combined under the assumption of a signal composition as predicted in the SM. If no further information on the signal efficiencies is given (see below), the application of these search limits to a model requires a SM-likeness test of the parameter point (see Ref. [20] for details). In short, this test checks whether the model-predicted signal composition is similar to the SM prediction, where all relevant search channels quoted in the analysis are taken into account, and their inclusive cross sections are used as weights in the determination of the maximally allowed deviation of the individual channel signal strength from the total signal strength. The details of the SM-likeness test procedure are described in Ref. [20] and have not been changed in HiggsBounds-5. After the LHC discovery of a SM-like Higgs boson, the focus of Higgs searches has somewhat shifted towards more model-independent, less combined search channels. Nevertheless, there are still many searches that combine different production or decay modes, assuming the relative contributions to be equal to those predicted in the SM, as e.g. motivated by predictions in pure scalar singlet extensions of the SM with a non-zero singlet-doublet mixing (see e.g. Refs. [215][216][217][218][219][220][221][222] for phenomenological studies in the LHC Run-2 era).
The SM-likeness test would not be needed if more information was provided publicly for the considered experimental analysis. The limit is typically reported on an inclusive cross section, σ tot , often also as signal strength, µ, i.e. normalized to the corresponding SM prediction. 6 For 5 As described in Section 2.3, exceptions to this assumption are possible in specific cases. 6 In special cases where no signal rate limit can be constructed HiggsBounds can also implement limits on other quantities that rely on additional model assumptions. This is currently the case for the CMS gg → φ → tt search [205] that constrains the effective φtt coupling. a combination of search channels i = 1, . . . , N , this signal strength can be calculated as where σ i (σ SM i ) denotes the inclusive signal rate for search channel i -comprised of one production and one decay mode -in the model (SM), respectively, and i ( SM i ) is the signal efficiency of channel i in the analysis, i.e. the fraction of signal events that pass the event selection, as predicted in the model (SM). If the three basic assumptions listed in the bullet points above are fulfilled, we have -to a good approximationi = SM i . A complication arises, however, if the experimental analysis does not provide information of the signal efficiencies of the involved signal channels as predicted in the SM, SM i . Unfortunately, it has been common practice to not release this information until now. 7 If the SM Besides the usual 95 % C.L. limits, HiggsBounds contains additional exclusion likelihood approximations for several cases. These are available for the main Higgs boson search channels as well as the SM and MSSM Higgs boson search combinations at LEP [20,41]. Moreover, HiggsBounds reconstructs the exclusion likelihood from BSM Higgs boson searches in the τ + τ − final state by ATLAS [127,223] and CMS [164,196], based on the numerical results presented in a single narrow resonance parametrization with two production modes. More details will be given in Section 4.3.

Recommendations for the Presentation of Future Search Results
Recently, a joint effort within the experimental and theoretical community led to the release of recommendations for the presentation of experimental results [224]. Following and partly expanding upon these recommendations, we propose the following guidelines for the publication of limits for LHC searches for new scalar bosons: 1. upper limits on the cross sections of the signal processes should be presented as a function of all relevant kinematical parameters, e.g. the masses and total widths of the involved scalar boson(s); 2. the search results should always contain the expected and the observed limit; 3. if the signal is comprised of several signal channels (i.e. different production and/or decay modes), the limit is set on a common scale factor -the signal strength µ -or a total signal rate. In this case, the signal efficiency of each signal channel should be provided as a function of all relevant kinematical parameters (see point 1); 4. if the limit is presented as a normalized signal rate (e.g. to the SM prediction), the reference signal rate should be quoted by the experimental analysis along with the result, thus enabling the recalculation of the limit on the signal rate's absolute value; 5. the search limit should always be presented at 95 % C.L.; 6. in addition, it would be beneficial to present results as exclusion likelihoods, using the same parametrization as the one used for the 95 % C.L. upper limit (see point 1).
These guidelines should result in a format of the search limit that is to a large extent modelindependent, in the sense that all dependences on kinematic parameters are fully described and can thus be incorporated in a recast of the limit onto specific models. Presenting both the expected and observed result enables a well-defined selection of the most sensitive analysis out of many search results -as already done in HiggsBounds-such that the derived global exclusion result can be interpreted at the 95 % C.L.. As already mentioned in the previous subsection, quoting the signal efficiencies is necessary for the proper determination of the signal strength in the model if several signal channels are involved in the signal process, and would be a better alternative to the SM-likeness test that needs to be applied otherwise.
As already mentioned, the BSM Higgs searches in the τ + τ − final state by CMS [164,196] and ATLAS [127,223] pioneered the publication of (multi-dimensional) exclusion likelihoods in addition to the usual 95 % C.L. upper limits at the LHC. These likelihoods were presented for a simplified model of a single narrow scalar resonance φ, parametrized in terms of its mass, m φ , the signal rate for single scalar production, σ(pp → φ → τ + τ − ), and the signal rate for scalar production in association with bottom-quarks, σ(pp → bbφ → τ + τ − ) (see Section 4.3 for details). This likelihood information has already turned out to be very useful in various phenomenological analyses (see e.g. Refs. [225][226][227][228]). Therefore, we strongly encourage the publication of exclusion likelihoods (in a similar form) also for other BSM Higgs search channels. We list in Table 6 a set of search channels, along with the relevant kinematic and signal rate parameters, which we deem suitable for providing this public information.
The  We also appreciate the efforts of using simplified workspaces provided by the experimental collaborations. In this approach it might be possible to retain the dominant theory nuisance parameters in the likelihood, while all experimental nuisance parameters that are not correlated with the dominant theory nuisance parameters are marginalized. Extensive tests would be necessary to explore the feasibility of this approach and the potential gain in information and precision over the current approach. A simplified approach based on the JSON format [239] has been proposed where only the most relevant theoretical (and, if necessary, experimental) systematic uncertainties are retained separately, and all other uncertainties are combined.
BSM Higgs boson searches at the LHC have so far considered inclusive signal processes, or, at least, have presented the result as a limit on the inclusive cross section. In the future a possible new path in the presentation of search results could be the presentation of limits on signal rates in specific phase space regions -so-called fiducial signal rates -instead of unfolding the result onto the inclusive signal rate. This is analogous to measurements of the discovered Higgs boson's signal rates, where strong efforts have recently been made to define specific phase-space regions for measurements in order to reduce the theory-dependence introduced in the unfolding process. In a similar way, one could think of a generalization of the Simplified Template Cross Section (STXS) framework [240] to the case of search limits. As HiggsBounds is used as a framework for HiggsSignals, which already incorporates STXS measurements of the Higgs boson signal, it would be straight-forward to implement corresponding phase-spacedependent limits also in HiggsBounds.
Lastly, it would be desirable in the future to improve and automize the implementation of experimental search results in HiggsBounds. The first step is that the experimental collaborations provide search limits in a machine-readable format, e.g. via their TWiki pages or via HEPData [241]. This is already done in many cases. As a next step, it would be useful to define a common data format that contains all necessary information about the search limit. Such data files can then be read in automatically by HiggsBounds. Such a data interface would also allow the HiggsBounds user to select specific search limits for their study in a very versatile way.

New Features in HiggsBounds-5
The most important improvement in HiggsBounds-5 over its predecessor HiggsBounds-4 is the inclusion of experimental results from the 13 TeV LHC. However, we will not discuss the newly implemented search limits in detail. Instead, we refer to the output of the AllAnalyses executable (see Section 5.2) to get a complete list of the experimental results included in a specific version of HiggsBounds. Instead, this section will describe the most relevant functional changes within HiggsBounds-5.

Effective Coupling Approximation for φ i V Production
In previous versions of HiggsBounds the cross section ratios for neutral Higgs boson (φ) production in association with a massive gauge boson, pp/pp → V φ (with V = W ± , Z), were obtained in the effective coupling approximation purely from the effective φV V coupling. In HiggsBounds-5 we have extended this approximation beyond the leading order by including contributions proportional to the scalar and pseudo-scalar Higgs couplings to top and bottom quarks, as well as interference effects.
The production of a neutral scalar boson φ in association with W -bosons always (up to nextto-leading order in QCD) takes place via Higgs-strahlung and is necessarily dependent on the coupling g φW W of the respective particle to the W bosons. At next-to-next-to-leading order (NNLO) in QCD, corrections from virtual top-quark loops arise which depend on the scalar coupling to top-quarks, g s, φtt . In case of Zφ production, additional important box-diagrams from the partonic process gg → Zφ need to be accounted for.
In terms of the effective couplings for a scalar particle φ, Eq. (1) and Eq. (3), we define The cross sections can been expanded as follows: Note that W φ production is largely dominated by the leading-order Higgs-strahlung process and only gets minor corrections from virtual top-quark loops. This is why the effect of bottom quarks has been neglected and only the κ W κ t interference term is considered in Eq. (8).
In contrast, for Zφ production all possible combinations κ a κ b have been included. For this expansion we neglect the effects from other possible scalar bosons in the model that may contribute due to non-vanishing φφ Z couplings. If these contributions turn out to be relevant in the investigated model, we advise the user to directly provide the hadronic cross sections instead of using the effective coupling approximation.
We calculate the inclusive W φ and Zφ production cross sections with VH@NNLO-2.0 [242,243] at NNLO in QCD. The mass-dependent expansion coefficientsσ V φ ab [m φ ] are determined by using the CP-violating 2HDM implementation of VH@NNLO to calculate various cross sections for different combinations of the effective couplings and solving the resulting system of linear equations. They are symmetric under a ↔ b.
VH@NNLO does not evaluate contributions from bb → Zφ for CP-mixed scalars which could lead to sizable differences in scenarios with large κb. However, within our approximation g φφ Z = 0 at tree-level holds for a pure CP-even scalar H and CP-odd scalar A with respective effective couplings κ b and κb (and equal masses). We therefore use the VH@NNLO SM Higgs boson implementation to determine the bb → HZ contribution, and consider it for both theσ Zφ bb andσ Zφ bb term in Eq. (9). The cross section calculation for a pure CP-even scalar particle H in VH@NNLO is more precise than the calculation for a CP-mixed scalar boson φ. For a SM-like CP-even Higgs boson Eq. (9) therefore produces a less accurate result than the dedicated calculation for a pure CP-even scalar boson in VH@NNLO. To circumvent this problem, we apply a K-factor approach and rescale the cross section for the CP-mixed scalar φ by the ratio of the more accurate CP-even calculation σ ZH and the scalar terms of the CP-mixed calculation σ φZ , i.e. with The definition in Eq. (11) ensures that for a pure scalar, i.e. for κt = κb = 0, our approximation coincides with the more accurate SM-like Higgs boson calculation in VH@NNLO. For a SM-like Higgs boson with all κ i = 1, the K-factor in Eq. (12) ranges between 1 and 2, with values larger than 1.1 only appearing for m φ 400 GeV. For gauge-phobic particles with κ Z = 0, K is typically of the order of 2 and can range up to 5 for masses below 10 GeV.
The coefficientsσ ab in Eqs. (8) and (9) are shown in Figs. 1 and 2 as a function of m φ . Figure 1 shows the two contributions to the pp → W φ process. Figure 2 shows the contributions to pp → Zφ grouped into pure CP-even (left) and pure CP-odd (right) contribtions. The dominant CP-even contribution isσ ZZ until the tt and Zt contributions become similarly relevant  for m φ 500 GeV. The CP-oddtt contribution is nearly identical to the tt contribution and is the largest contribution for the CP-odd case, where ZZ contributions are absent. Note that the bb andbb contributions contain the bb → Zφ subprocess in addition to the b-quark boxes of the gg → Zφ subprocess. All cross-terms proportional to the product of a CP-even and a CP-odd coupling vanish in our parametrization. These can only originate from the contribution of additional Higgs bosons, which is neglected here.
Note that our expansions in Eqs. (8) and (9) have the following limitations in their applicability to BSM Higgs models: • In all processes, the higher-order correction terms only account for virtual top-and bottom-quark loops, thus we assume that no other particles with similar quantum numbers run in the loop and give a significant contribution. Models with additional, colorcharged particles with scalar interactions, as e.g. the scalar top squarks in supersymmetric theories, may require a proper calculation of the quantum corrections of the additional particles.
• As scalar-scalar-vector interactions, g model φφ V , are not accounted for, care must be taken in models with several light scalar degrees of freedom with non-negligible interaction between these particles. These may contribute in s-channel diagrams with the additional scalar particles as propagators.
The approximation described here is a very significant improvement over the effective coupling approximation for pp → V H in previous versions of HiggsBounds. It is automatically used in the effective coupling input. In models where the assumptions above are satisfied it can also be used to substitute an explicit calculation of σ(pp → V H) for the hadronic input scheme. In this case, the approximated hadronic cross sections can be accessed through the functions in the access_effC.f90 file. Furthermore, this approximation is not only used for the inclusive σ(pp → ZH) cross section but also provides separate σ(qq → ZH) and σ(gg → ZH) cross sections that may be kinematically separable in some analyses.

Direct Charged Higgs Production
As discussed in Section 2.2, we have added many charged Higgs search channels to HiggsBounds that are or may -in the future -be probed at the LHC. The most thoroughly studied (and in many cases dominant) production channel is the production of a charged Higgs in association with a top-quark -denoted pp → H ± tb in the four-flavor scheme or bg → H ± t in the five-flavor scheme -which is typically considered in experimental searches at the LHC, see e.g. Refs. [130,134,199,244]. The cross section has been calculated including NLO-QCD corrections in the 2HDM and MSSM [240,[245][246][247][248][249].
While the 2HDM is the simplest BSM model with a charged Higgs boson H ± , the coupling structure of its charged-Higgs-quark couplings readily generalizes to a large variety of models. In the 2HDM (and also e.g. in the MSSM at tree-level) the relevant charged Higgs coupling to top and bottom quarks has the form with κ ± t = 1/ tan β, κ ± b = tan β for flavor conserving 2HDM Yukawa sectors of type II (or in the MSSM) and κ ± t = κ ± b = 1/ tan β in Yukawa type I. For a generic charged Higgs boson H ± j we expand the cross section in terms of κ j± t,b defined as in Eq. (13) and obtain We keep the interference term -even though it is suppressed by an additional mass insertion due to the helicity structure of the coupling -as its contribution becomes important for charged Higgs masses close to the top-threshold region.
For charged Higgs bosons lighter than the top-quark, the internal top-quark propagators can go on-shell introducing a dependence on the width of the top-quark, Γ t . According to Ref. [249], this can be approximately included by rescaling the cross section as Neglecting decays of the top-quark into first and second generation quarks in the SM, this ratio of widths is simply given by BR(t → W + b) in the BSM model under consideration. We therefore parameterize σ tH − j in the entire m H ± j range as where We use the results 8 of Refs. [240,249] tabulated in the 2HDM type II as a function of the charged Higgs mass and tan β and extract the mass dependent coefficientsσ tH − j ab by solving the resulting system of linear equations. In the region m H ± j < m t we use HDECAY-6.52 [250,251] to calculate the required branching ratios of the top-quark. The resulting parametrization reproduces the original results to a relative accuracy of better than 10 −4 for m H ± j > 2m t and -with deviations of at most 2 % -stays well within the theoretical uncertainties of the original calculation [249] even for m H ± j < m t .
The parametrization in Eq. (16) holds for any charged Higgs boson H − j that has a coupling structure of the form Eq. (13), and is valid as long as no other BSM effects contribute up to NLO in QCD. In particular this neglects possible contributions of the form pp → bbφ → bbH + W − that can appear in 2HDM-like models. However, these are typically only relevant for resonant φ production where they are treated separately. SUSY QCD corrections also impact these results, and the dominant ∆ b corrections can be included through a rescaling of tan β. 8 In the region m H ± j < m t the approximation relies on the assumption BR(t → W + b) + BR(t → H + j b) ≈ 1 (no sum over j). If this assumption is violated -e.g. because the top quark decays into multiple H ± j or into additional new-physics decay modes -the threshold behavior would be incorrect and a full model-specific calculation should be performed. However, heavier H ± j are insensitive to the top width, and valid cross sections for any number of H ± j heavier than the top quark can be obtained.
In HiggsBounds-5 this approximation can be accessed through the HCCS_tHc function in the access_effC.f90 file, which requires m H ± , κ j± t , κ j± b , and BR(t → H + j b) as input and assumes The function returns the cross section as the cross section is charge-symmetric, σ tH + j = σt H − j . This inclusive value then corresponds to the required HiggsBounds input quantity CS_Hpmjtb, see Section 2.2.

Exclusion Likelihoods for LHC Higgs to τ + τ − Searches
In experimental searches for additional Higgs bosons decaying into τ + τ − the ATLAS [127,223] and CMS [164,196] collaborations have released simplified exclusion likelihoods as a function of the two contributing single Higgs production modes, gg → φ and gg → bbφ, and for a wide range of narrow scalar resonance mass hypotheses. The implementation of these nearly model independent likelihoods in HiggsBounds includes an approximate scheme for treating multiple contributing Higgs bosons of similar mass. The implementation of the first analysis from LHC Run-1 [164] for which this input was provided has been described in detail in Ref. [21]. We present the implementation and validation of new LHC Run-2 analyses and discuss improvements to the derivation of exclusion limits from the provided likelihood information. More details on the underlying likelihood reconstruction method can be found in Ref. [21].
The profiled likelihood analyses underlying the experimental results use the test statistic where N is the observed data, b is the background expectation, and s(m) is the signal expectation for a given hypothesized resonance mass m and given contributions of the two sub-channels. A limit is set on the signal strength modifier µ in the presence of the globally optimized nuisance parametersθ, the globally optimized signal strengthμ and the conditionally optimized nuisance parametersθ µ for the given value of µ. The experiments provide expected, q exp µ , and observed, q obs µ , values for this test as a function of the two sub-channel contributions for different resonance masses.
Since the likelihood was parametrized in terms of a single narrow scalar resonance, in a specific model application in HiggsBounds multiple Higgs bosons potentially contributing to the signal have to be combined and mapped onto the likelihood parametrization. This is done in HiggsBounds by a clustering algorithm. For this all Higgs bosons within are combined into a cluster, their rates are incoherently summed 9 , and a signal-rate-weighted cluster mass is used to approximate the mass of the single resonance mass. The numerical coefficient ∆ res is chosen to approximately match the mass resolution of the τ + τ − channel under consideration, and is currently set to 20% for all implemented analyses. This algorithm has already successfully been applied in various analyses, including cases in which more than two Higgs bosons form a cluster [21,29,225].
In the limit of large numbers the test statistic q µ can be treated as a ∆χ 2 above minimum such that 68 % C.L. and 95 % C.L. exclusion bounds are obtained at q µ = 2.28 and 5.99, respectively, corresponding to a two-sided limit (or fit). This approach was employed by the experimental collaborations to obtain the confidence regions in the presented two-dimensional cross section planes for fixed resonance mass, and hence also used to obtain 95 % C.L. exclusion limits from the likelihood information in HiggsBounds-4 [21]. More appropriate for limit setting, however, is a one-sided upper limit on the signal cross section. Therefore, HiggsBounds-5 uses an improved approach in which CL s is directly calculated from the provided likelihood information and the 95 % C.L. allowed region is obtained at CL s > 0.05 = 1 − 95 %. This reconstruction relies on the fact that the quantity q exp µ provided by the experiments can be interpreted as 10 where σ is the expected effective Gaussian uncertainty of the signal strength modifier µ. The expected and observed CL s can then be obtained just from q exp µ and q obs µ in the asymptotic limit as using [252] CL exp for the expected limit and for the observed limit. In all of these, Φ denotes the cumulative normal distribution function. The experimental collaborations use the CL s > 0.05 criterion for model-specific limit setting, e.g. in the context of MSSM benchmark scenarios. HiggsBounds-5 uses the improved approach that directly employs CL s > 0.05 instead of ∆χ 2 < 5.99 both for determining the sensitivity of the analysis via CL exp s and for obtaining the 95 % C.L. limit via CL obs s . This improved methodology in HiggsBounds that more closely resembles the one employed by the experimental collaborations enables a reliable reconstruction of limits from the provided likelihood information. As we will demonstrate below, with HiggsBounds it is now possible to  pp → H/A → τ + τ − searches at ATLAS [223] (top panels) and CMS [196] (bottom panels) in the M h 125 scenario [29]. The solid black and white lines show the reconstructed 95 % C.L. limit in HiggsBounds using the old (black, based on ∆χ 2 < 5.99) and the improved new method (white, based on CL s > 0.05). For the limit based on CL s > 0.05 the white-dotted lines indicate the variation of the limit due to signal-model-dependent theory uncertainties. The red dashed line shows the official ATLAS or CMS 95 % C.L. limit. even reproduce and understand methodical differences between the official ATLAS and CMS model interpretations. Figure 3 shows a validation of the HiggsBounds likelihood reconstruction algorithm against the most recent 13 TeV experimental analyses by ATLAS [223] (at 139 fb −1 , top) and CMS [196] (at 36 fb −1 , bottom) in the M h 125 scenario [29] of the MSSM. The theoretical predictions for the scenario are taken from the LHCHXSWG, based on the following prescription [240,253]: FeynHiggs [254][255][256][257][258][259][260] is taken for the calculation of the MSSM masses and couplings, a combination of FeynHiggs and HDECAY [250,251] for the τ + τ − branching ratio. The gluon fusion cross section predictions are obtained with SusHi [261,262] including all available higher order corrections [263][264][265][266][267][268][269][270][271][272][273], and the bb associated channel uses matched cross section predictions [274][275][276][277][278][279][280]. The prescription for the theoretical predictions in the MSSM benchmark scenarios is such that the model-dependent theoretical uncertainties of the predicted masses, cross sections and branching ratios should be incorporated by the user through appropriate variations of the theoretical predictions. We will discuss the impact of these theoretical uncertainties on the obtained limits below.
The color code in Fig. 3 shows the expected (observed) reconstructed likelihood value −2 ln(L) from HiggsBounds on the left (right). The black contour indicates the HiggsBounds 95 % C.L. exclusion limit reconstructed using the old ∆χ 2 < 5.99 approach, while the solid white contour displays the 95 % C.L. limit from the new CL s method. The reconstructed limits displayed by the solid black and white contours do not take into account signal-model-dependent theoretical uncertainties. For the limit based on the CL s method the white-dotted lines indicate the uncertainty band around the solid white contour according to the model-dependent theoretical uncertainties in the considered M h 125 benchmark scenario of the MSSM [29], see below. For comparison, the red-dashed lines show the official 95 % C.L. limits from ATLAS (upper panels) and CMS (lower panels).
We start by comparing the limits obtained with the ∆χ 2 and the CL s method (solid black and white contours). The CL s > 0.05 criterion results in a larger excluded area compared to the ∆χ 2 < 5.99 criterion. This feature is expected and can be understood as follows. The ∆χ 2 = 5.99 limit corresponds to a 95 % C.L. limit on CL s+b in the Gaussian approximation. However, in order to prevent erroneous exclusions in regions where the search has no sensitivity [281] CL s is constructed by dividing CL s+b by CL b , see Eq. (22). Since the expectation value of CL b in the absence of any signal is CL b = 0.5, a 95 % C.L. limit on CL s approximately corresponds to a 1 − 0.05/ CL b ≈ 90 % C.L. limit on CL s+b . 11 We now compare the limit obtained with the CL s method with the results obtained by ATLAS and CMS in their analyses for the M h 125 benchmark scenario. As explained above, the reconstructed limit in HiggsBounds is based on the (nearly) model-independent likelihood provided by ATLAS and CMS which by construction does not contain any model-specific theoretical uncertainties on the cross sections and branching ratios of the signal processes. Accordingly, the solid white contour corresponds to the limit obtained with the CL s method without taking into account model-specific theoretical uncertainties. We find that the resulting limit agrees almost perfectly with the one obtained in the ATLAS analysis, both for the expected (left upper panel) and the observed limit (right upper panel). On the other hand, the benchmark analysis of CMS excludes a smaller region than in our CL s analysis (solid white contour), both for the expected and the observed limit (lower panels). As we have verified via direct communication with members of the ATLAS and CMS collaborations [282], this feature can be understood from the fact that the experimental interpretation of the benchmark scenario by CMS includes model-specific theoretical uncertainties on the signal cross-sections, while in the ATLAS analysis no such signal-model-dependent theoretical uncertainties have been taken into account.
As a final step of this comparison we now take into account signal-model-dependent theoretical uncertainties with HiggsBounds. The dotted white contours indicate the uncertainty band around the solid white contour that has been obtained by running HiggsBounds with input rates at the upper and lower end of the theoretical uncertainties and interpreting the resulting difference in the 95 % C.L. exclusion as a theoretical error band. In line with the CMS analysis, we include scale and parton distribution function uncertainties on the cross section predictions, provided by the LHCHXSWG for the M h 125 scenario. Theoretical uncertainties on the BRs are not considered since they are negligible by comparison [240]. The upper branch of this band shows how the limit is weakened by the incorporation of the signal-model-dependent theoretical uncertainties. We find that this contour agrees well with the limit that has been obtained in the CMS benchmark analysis. We expect that an ATLAS limit incorporating the signal-model-dependent uncertanties would be less constraining and closer to the upper white-dotted limit obtained with HiggsBounds. 12 By default, the likelihood information is used to reconstruct a 95 % C.L. limit that is then treated like any other limit in HiggsBounds. HiggsBounds selects the Higgs cluster giving the largest expected exclusion likelihood as the most sensitive Higgs boson combination to be tested against the observation. The advantage of reconstructing the limit from the likelihood is that the full efficiency information on the two involved production channels is incorporated. This is especially important in the BSM φ → τ + τ − channel, since the relative contributions of gg → φ and gg → bbφ can change drastically through the model parameter space, e.g. in the MSSM or the 2HDM. The value of the likelihood can also be accessed directly through the HiggsBounds_get_likelihood subroutines (see online documentation), such that it can be included in a global likelihood analysis of BSM models, see e.g. Refs. [225][226][227][228].
The implementation of the recent ATLAS result [223] raised a different kind of issue, as it contains an observed excess with a local significance of more than 2σ for masses in the range 400 GeV to 500 GeV. Therefore, the parameter point of zero signal rates lies outside the 95 % C.L. contour of the provided likelihoods. A direct application of these results in the original form to a model would exclude parameter points that feature scalar boson(s) with small or vanishing pp → φ → τ + τ − rates in that mass range. At the same time, parameter points with no scalars in that mass range would not be excluded, as the likelihood would not 12 In Fig. 3 the black contour indicating the HiggsBounds limit using the ∆χ 2 method in all cases happens to be close to the limit that incorporates the signal-model-dependent uncertainties (and thus close to the official CMS results). We stress that this is a scenario-dependent coincidence and that the theoretical uncertainties are not captured by the ∆χ 2 limit. be evaluated if no scalar boson is within the mass range. We therefore use an approximate approach to avoid this inconsistent behavior while keeping the overall features of the likelihood profile intact. This approach is applied to the mass range with an excess of more than 1σ, i.e. for the mass planes at 350, 400 and 500 GeV provided by ATLAS. It restores the property q obs µ = 0 forμ > µ required of a test statistic used in limit setting [252], whereμ is the best fit rate.
We use the M φ = 400 GeV mass plane, shown in Fig. 4, to illustrate the approach. The grey contours are the 95 % C.L., ∆χ 2 = 5.99 contours 13 of the original (processed, i.e. based on our approach for avoiding the exclusion of too small BSM rates) observed likelihood profile q obs µ shown on the left (right). We first fit an ellipse centered at the best fit point (grey cross) to the 95 % C.L. contour of the original q obs µ . In Fig. 4 (left) this ellipse is shown in black. We then construct the ellipse shown in red, which is centered at the origin with axes parallel to the coordinate axes. We fix the eccentricity by requiring this ellipse to be tangential to the long axis of the black ellipse in the best fit point (the black-dashed line). We consider all signal rates on this red ellipse to be equal to the best fit rate, µ ∼μ. 14 Accordingly, the 13 We use the simpler ∆χ 2 approach only to construct the processed likelihood tables. The CLs method described above is always used to obtain 95 % C.L. limits in HiggsBounds-5. 14 To first approximation, the ratio between the long and short axes of this ellipse resembles the ratio of signal efficiencies in the two production channels. Therefore, our construction indeed approximately determines likelihood inside the red ellipse is set to zero, since lower predicted rates than at the best-fit point should not be disfavored in a limit setting procedure. To obtain a smooth transition, we introduce polar coordinates (r, θ) and, for each angle θ, approximate the likelihood profile by where r E (r 95 ) is the radial component of the red ellipse (the grey 95 % C.L., ∆χ 2 = 5.99 contour) for the given θ. Outside the grey contour the likelihood remains unchanged from the original. This leads to the processed likelihood profile shown in the right panel of Fig. 4.

User Operating Instructions
For HiggsBounds-5 we have made substantial changes on the structure and online platform of the HiggsBounds source code. The code has moved to a GitLab repository and is now available at https://gitlab.com/higgsbounds/higgsbounds. This modernization effort also included moving to CMake as the build system. HiggsBounds is now compiled by running: mkdir build && cd build cmake .. make The only requirements are CMake and a Fortran compiler. This will compile the library, the main executable and a number of example programs that illustrate different use cases. More detailed information on building and linking HiggsBounds can be found on the abovementioned website.
the parameter region where the total fiducial signal rate is equal to the one at the best fit point.

Fortran Subroutines
The Fortran subroutines are the most powerful and versatile way of using HiggsBounds. Upto-date descriptions of the various Fortran subroutines and functions can be found in the online documentation at https://higgsbounds.gitlab.io/higgsbounds.
Compared to HiggsBounds-4, all of the input subroutines have been extended to include all of the quantities discussed in Section 2. The input subroutine arguments are named as in Tables 1 to 5 and are either double precision or integer arrays with dimensions given by the number of neutral Higgs bosons and/or the number of charged Higgs bosons.
As a further improvement, HiggsBounds-5 includes a C interface to all of the Fortran subroutines to facilitate the use of HiggsBounds from C or C++ codes. This interface automatically handles type conversion and accounts for the different storage orders of multidimensional arrays between C and Fortran. The C interface is included in the online documentation.

Command-line Version
Compiling HiggsBounds generates a main executable that can be run as

./HiggsBounds <whichanalyses> <whichinput> <nHzero> <nHplus> <prefix>
where the arguments specify the following: <whichanalyses> specifies which experimental data is selected for the model test -'LandH' for all implemented results, 'onlyL' for LEP results only, 'onlyH' for hadron-collider results only, and 'onlyP' for published results only; <whichinput> specifies whether the model input on the production and decay rates is provided in the effective couplings approximation ('effC'), at the cross section level ('hadr'), or via an SLHA input file ('SLHA'). The arguments <nHzero> and <nHplus> specify the number of neutral and charged Higgs bosons, respectively. The argument <prefix> denotes the path to the input files including any part of the filename that is common to all input files.
Additionally an executable called AllAnalyses is generated. It prints a

HiggsBounds Data Files Input
If HiggsBounds is run from the command line with the option <whichinput>=effC or hadr the model input needs to be specified via HiggsBounds specific input files. These are whitespace separated tabular text files containing the input quantities for one datapoint per row. With respect to HiggsBounds-4, the input files have been adjusted to the changes in the input quantities detailed in Section 2.
An overview of all data input files and their data structure is given in Table 7. For some higher-dimensional arrays only some elements have to be specified, as will be explained below. Table 7 also specifies whether the data file is required for either of the two HiggsBounds input schemes (effC or hadr), or used as optional input. If a required file is not provided as input, HiggsBounds warns the user but proceeds to run while setting the unspecified input quantities to zero. In Table 7 should be specified in the input file in the order    Higgs production cross sections (in pb, arbitarity values). The cross sections for charged Higgs production in association with one or two light flavor quarks (u, d, s) are generally combined to inclusive production processes containing one or two untagged jets. For the vector boson fusion process we set both quark PDG numbers to 1 in order to differentiate it from the other quark-associated production processes. All cross sections correspond to the sum of H + and H − production, hence, all PDG numbers are taken to be positive (except for H + H − production).
of the decaying Higgs boson). The N 2 h 0 (N h 0 − 1)/2 non-redundant elements can be specified in the following way: For every k ∈ {1, N h 0 } we specify the lower left triangle (including the diagonal), but with the kth column and kth row removed, e.g. for N h 0 = 3 The input arrays BR_hjHpiW, BR_HpjhiW and CS_Hpjhi are not reducible and should be specified row by row in the input files.

SLHA
In HiggsBounds-4 the squared SM-normalized effective Higgs couplings to bosons and third generation fermions were provided in the two SLHA blocks HiggsBoundsCouplingInputBosons and HiggsBoundsCouplingInputFermions, respectively. Since HiggsBounds-5 requires the sign information for the effective couplings, we have replaced these blocks by very similar blocks named HiggsCouplingsBosons and HiggsCouplingsFermions containing the non-squared, sign sensitive effective Higgs couplings as described in Section 2. In case only the old blocks are specified in the SLHA input file for HiggsBounds-5, the effective Higgs couplings are taken to be the positive square-root of the given values.
In addition, we introduced an SLHA input block containing the hadronic cross sections for direct charged Higgs boson production. In absence of a corresponding SLHA convention, we call these input blocks ChargedHiggsLHC8 and ChargedHiggsLHC13 for the predictions for the LHC at 8 and 13 TeV, respectively. 15 The first three columns specify the final state particle PDG numbers in increasing order (modulo a sign in case of anti-particles). In case of a twobody final state the first column is filled by a zero. The fourth column gives the cross section in pb. An example (employing the particle spectrum of a 2HDM Higgs sector) for one of these SLHA blocks is given in Table 8.

Summary
This paper documents a major update of the public Fortran code HiggsBounds which tests general BSM models against exclusion limits from LEP, Tevatron and LHC Higgs searches. We have presented the theoretical input framework of HiggsBounds-5, which has been significantly extended to allow predictions for all current and many potential future Higgs search channels. In particular this extension adds sub-channels for several production modes -such as qqand gg-induced Zh production -that may be kinematically separable, and incorporates flavor-violating decay modes and decays into BSM particles. The charged Higgs boson input framework has also been extended by many different direct production processes, some of which are already probed at the LHC.
We discussed the main experimental input for HiggsBounds -the (nearly) model-independent upper cross section limits -and the possible limitations of their application to BSM models. In fact, many of these limitations in current search results can be overcome if more detailed information, in particular on the signal composition in terms of Higgs production processes, is released publicly by the experimental collaborations. Therefore we suggested guidelines for the publication of experimental search results that we deem essential for a proper reinterpretation in terms of BSM models. These recommendations are in line with and partly extend those presented in Ref. [224].
In many BSM models precise calculations for "exotic" production cross sections are in many cases not readily available. Therefore, for two of the important production modes -neutral Higgs production in association with a massive gauge boson and top-associated charged Higgs production -we have added model-independent parameterizations of existing calculations. The effective coupling (or scale factor) approximations of the Zh and W ± h production cross sections are based on results obtained with the code VH@NNLO [242,243] and include CPsensitive contributions. The tH + cross section parametrizes the precise calculations in the 2HDM [240,[245][246][247][248][249] through model-independent coupling scale factors.
In most searches for additional Higgs bosons the final result provided by the experimental collaborations is a -potentially multidimensional -upper limit on the cross section at 95 % C.L. as a function of the relevant kinematic variables, i.e. the masses and widths of the involved particles. However, for the neutral Higgs boson searches in the τ + τ − final state a simplified exclusion likelihood was provided by both the ATLAS and CMS collaborations. These likelihoods are implemented in HiggsBounds, and the resulting likelihood value is made available for use in model fits. We encourage the release of such likelihoods also for other experimental search channels in the future [224]. We have improved the derivation of 95 % C.L. limits from provided likelihood information by using the CL s method and presented a procedure to prevent overexclusion in case of excesses in the searches. A validation of the likelihood implementation in the M 125 h scenario of the MSSM using these new techniques found very good agreement with the official results from CMS (ATLAS) taking into account (ignoring) the model-dependent theoretical uncertainties on the signal rates.
HiggsBounds-5 also involves substantial technical changes. The code is now available in a public git repository at https://gitlab.com/higgsbounds/higgsbounds with updates being released whenever new analyses have been implemented. We have modernized the build system to use CMake and added a C interface to make it easier for other codes to link to HiggsBounds. A technical description of the user subroutines and further details on the code are given in the online documentation at https://higgsbounds.gitlab.io/higgsbounds .