Testing Exotic Scalars with HiggsBounds

The program HiggsBounds is a well-established tool for testing Beyond-the-Standard Model (BSM) theories with an extended Higgs sector against experimental limits from collider searches at LEP, Tevatron and LHC. Thus far, it could be applied to any neutral or charged Higgs bosons originating from the modified Higgs sector. Implicitly, these particles were assumed to exhibit a somewhat hierarchical Yukawa structure as present in the Standard Model, where in particular the couplings to first generation fermions could be neglected. In this work, we extend the HiggsBounds functionalities to go beyond these restrictions, thus making the code applicable to any neutral or charged BSM scalars. Moreover, we develop a new approach to implement experimental searches whose kinematic acceptance depends significantly on the values of the involved couplings. We achieve this by recasting the searches to general scalar models. Using this approach we incorporate relevant current experimental limits from LHC searches for exotic scalars, and present the implications of these limits for a dark matter scalar mediator model, a flipped Two-Higgs-Doublet Model and a supersymmetric model with R-parity violation.


Introduction
The Higgs boson discovered at the LHC in 2012 is the first observed potentially elementary particle with spin-0, i.e., a scalar boson. In the Standard Model (SM) of particle physics, this scalar boson arises from augmenting the theory with a scalar SU(2) L doublet field and associated renormalizable and gauge-invariant terms in the scalar potential. For a specific constellation of the scalar potential parameters the Higgs mechanism leads to a spontaneous breaking of the electroweak (EW) symmetry and provides masses to the W ± and Z bosons. The detection of the Higgs boson is thus a crucial sign that this theory is at least approximately realized at the currently probed energy scales.
While the Higgs mechanism in the SM provides a plausible explanation for the broken EW gauge symmetry and the massiveness of its associated gauge bosons, both observational and theoretical puzzles -for instance, the presence of dark matter (DM) in our Universe, the quantum description of gravity, and the naturalness (or hierarchy) problem -remain unsolved in the SM framework, and thus motivate the postulation and experimental search for new physics beyond the SM (BSM).
Many of such BSM theories postulate the existence of additional scalar particles. A prime example is Supersymmetry (SUSY) [1][2][3] which introduces a new scalar boson for every SM fermion, leading to a plethora of new spin-0 states. Other models simply extend the SM scalar sector by additional scalar SU(2) L doublet or singlet fields, leading, for example, to the Two-Higgs-Doublet-Model (2HDM) [4,5]. Yet another approach are simplified BSM models, in which only a minimal particle content is introduced to address a particular problem. For instance, addressing the DM problem, one may only postulate the existence of a stable DM candidate particle and a mediator particle that interacts both with the DM and SM particles [6]. This mediator could e.g. be a neutral gauge boson ("Z boson") of a new gauge symmetry, or simply a new scalar particle.
The ATLAS and CMS experiments at the LHC follow various strategies to search for these additional scalar particles. In SUSY -under the additional assumption of conserved R-parity or proton hexality [7] -the new supersymmetric particles can only be produced pairwise at the collider, and would then decay (possibly via cascades) into the lightest supersymmetric particle (LSP), which is electrically neutral and stable, thus escaping the detector unseen. SUSY searches, therefore, typically look for inclusive production of highly energetic objects (leptons, jets) associated with a fairly large amount of missing transverse energy (MET) (see Refs. [8][9][10] for an overview of current searches). LHC searches targeting BSM models with an extended Higgs sector usually look for single production of a new scalar boson which either decays directly to SM particles (tt, bb, τ + τ − , µ + µ − , W + W − , ZZ, γγ, etc.) or to final states containing also other scalars (often including the discovered Higgs boson). The third class of BSM models mentioned above -the simplified DM modelsare often targeted by searches for a single highly-energetic object (jet, Z-boson, or Higgs boson, h 125 ) which recoils against the produced DM particles. The DM particles escape the detector and lead to MET [11].
Thus far, all these searches have not found any significant deviations from the SM background expectation. The results are, therefore, presented as upper limits on the signal cross section -which rely on fairly modest model assumptions, e.g. the spin and CP property of the involved new particle(s), specific coupling relations, or the decoupling of other BSM particles -or as excluded parameter regions within specific BSM benchmark models. While the cross section limits can be applied directly to models that feature the same process with particles that fulfil the same assumptions, exclusion regions in BSM benchmark models can generally not be re-interpreted within different models. Therefore, if the model under study does not strictly fulfil all the assumptions underlying the presented search limit, a painstaking recasting analysis based on detailed Monte-Carlo (MC) simulations is often needed in order to estimate the sensitivity of the search and to derive a corresponding upper cross section limit on the alternative signal process. Naturally, such a recasting analysis cannot be as accurate as if the alternative signal process was directly analysed by the experimental collaboration. Ideally, the experimental collaborations would provide limits for all possible signal processes and parameter constellations their search is sensitive to, or, alternatively, provide additional information on how efficiencies and signal acceptances change under different model assumptions. Of course, resources and manpower of the experiments are limited, and the number of potentially relevant alternative signal processes can be large, rendering a complete coverage of all model-interpretations unrealistic in most cases. 1 Several tools have been developed over the years to facilitate the assessment of the experimental validity of a new physics model, given the latest experimental limits from LHC searches. One such tool is HiggsBounds [13][14][15][16][17][18] allowing to test extended Higgs sectors against upper limits from BSM Higgs searches. In this work, we present an extension of HiggsBounds that enables the program to test scalar particles that have quite different properties (to be quantified in detail later) than normally expected of BSM Higgs bosons. While for most BSM Higgs models the additional scalar particles couple only weakly to first and second generation quarks, this is not a necessity in generic scalar extensions of the SM. For specific high-energy models accommodating scalars with large first and second generation quark couplings see e.g. Refs. [19,20]. Moreover, we test models involving quark-flavor and lepton-flavor violation.
First, this requires an extension of the HiggsBounds framework that handles the model predictions, in which hadronic cross sections are approximated from internal fit functions and effective couplings provided by the user. Second, new types of LHC search results had to be implemented. In particular, we considered limits from di-jet resonance searches [21][22][23][24][25][26][27], searches for di-jet resonances in association with an energetic photon [28] or lepton [29], and di-lepton resonance searches [30][31][32][33][34].
The ATLAS search for di-jet resonances in association with a photon from initial state radiation (ISR) is targeted to models with a Z boson. In order to make the limit applicable to scalar resonances we perform a detailed MC recasting analysis and derived cross section and signal acceptance functions for various constellations for the scalar-quark-quark couplings. All these results are incorporated in HiggsBounds in the form of simple fit formulae allowing the user to derive accurate and fast limits without running any MC generation himself. We describe our recasting analysis in detail in order to highlight important parameter dependencies and to motivate and guide the experimental collaborations to provide a corresponding signal interpretation in upcoming searches.
We demonstrate the extended HiggsBounds features and the impact of the newly implemented search limits in three exemplary model applications: first, a simplified DM model with a scalar mediator particle ("scalar DM portal model"); second, a 2HDM with large Higgs-b-quark couplings; third, we discuss the case of resonant scalar lepton (slepton) production and decay in SUSY models with R-parity violation.
This paper is organized as follows: In Section 2, we discuss the generic scalar model we use for the implementation of all relevant di-jet and di-lepton searches. The implementation of these searches is then discussed in Section 3 with a special focus on the di-jet plus photon search of Ref. [28] in Section 3.3. In Section 4, we discuss a simplified DM portal model, the 2HDM, and the R-parity violating MSSM as exemplary model applications. The conclusions can be found in Section 5. Appendix A provides details on the validation of our implementation of the ATLAS di-jet plus photon search. The derived fit formulas for the cross section and the acceptance can be found in Appendix B. Appendix C documents the new HiggsBounds routines implemented in the course of this work.

Generic scalar models
For the implementation of the di-jet and di-lepton limits, we employ generic scalar models. While these models are not meant to be complete BSM models, they are designed to allow for all relevant interactions.
For a neutral scalar S, we use the following generic scalar model, where the second and third line encode the couplings of S to up-type and down-type quarks. The couplings g q,ij are the CP-even Yukawa couplings to quarks; the couplingsg q,ij are the CP-odd Yukawa couplings to quarks. Note that we also allow for quark-flavor violating interactions. Similar to the quark couplings, we write down the couplings to leptons in the third line of Eq. (1) parameterized by the couplings g ,ij andg ,ij allowing for lepton-flavor violation.
The analogous generic scalar model for a charged scalar S ± reads where P L,R are the left-and right-handed chirality projection operators. The second line of Eq.
(2) encodes the interaction of the charged scalar with quarks allowing again for quark-flavor violation. The ellipsis in the last line of Eq. (2) stands for interactions of S ± with other SM particles (e.g. with a photon) or other BSM particles. The general scalar model of Eqs. (1) and (2) is available as FeynRules model file [35][36][37]. This implementation is based on the DMsimp_s_spin1 UFO model [38]. The model file is available as ancillary file accompanying the present paper.

Implemented searches
In this Section, we discuss current experimental searches for BSM resonances decaying to di-lepton or di-jet final states. We implement all searches that either directly provided cross section limits or at least gave efficiencies for scalar resonances. Additionally, we use MC simulations in the generic models of Section 2 to obtain and tabulate the required efficiencies for scalar resonances in several models where they were not provided by the experiments.

Di-lepton final states
Higgs bosons decaying into τ + τ − are among the most sensitive signatures in searches for models with extended Higgs sectors. As such, the existing searches in this channel, including the most recent CMS [39] and ATLAS [31] limits were previously implemented in HiggsBounds. 2 In contrast, decays into e + e − and µ + µ − pairs are strongly suppressed by their small masses for Higgs-like particles. While some searches involving µ + µ − final states [40][41][42][43][44][45] were already implemented in HiggsBounds, e + e − signatures were not previously included at all. To allow implementing e + e − searches in HiggsBounds, we extend the input framework to include this decay mode, see Appendix C for details. This extension allows us to implement the latest ATLAS di-lepton resonance search [30] which covers a huge, previously largely uncovered mass range from 250 GeV to 3000 GeV including finite width effects of up to 10 % of the resonance mass. The limit is set on a fiducial cross section and the required selection efficiencies for spin-0 resonances are available as auxiliary material in the analysis.
Limits on lepton flavor violation had previously only been implemented as limits on BR(h 125 → eµ/eτ /µτ ) [46,47]. We extended these by the latest lepton flavor violating resonance searches by ATLAS [32] and CMS [33,34] that all provide useable cross section limits for scalar particles.
For an overview of the newly implemented di-lepton searches see Table 1.

Di-jet final states
Searches for di-jet resonances probe the highest masses accessible at the LHC. For scalars with a SM-Higgs-like coupling structure, they are typically not very sensitive, since at high masses the decays into top-quarks usually far outweigh those into light quarks and gluons. An exception are di-b-jet final states, which are commonly probed in searches for bb → H → bb [48,49]. This signature can be particularly sensitive to e.g. the high tan β regions of type II (e.g. MSSM) or flipped Yukawa sectors. However, di-b-jet searches are also carried out   independently from dedicated Higgs searches, and we have implemented additional results from both ATLAS [22] and CMS [26]. Furthermore, ATLAS performed a di-jet -including di-b-jet -resonance search for a resonance produced in association with a γ from initial state radiation [28]. While the CMS di-b-jet result includes a useable interpretation for a scalar produced in either gluon fusion or bb-associated production, the two ATLAS results instead set limits on σ · A · for a Gaussian resonance. To implement these limits, the acceptances A and efficiencies for a particular signal model need to be derived from Monte-Carlo simulations. In Section 3.3, we discuss this procedure in detail for the di-jet + γ search. The acceptances for the bb search [22] were obtained analogously to the more complicated case described in Section 3.3.2.
Most searches for di-jet resonances that target gluon and light-quark jets -with the current strongest in different mass ranges being [21,[23][24][25]27] -only provide limits on a Gaussian resonance and no useable signal acceptance information for scalar signals. In principle, all of those searches could be implemented in HiggsBounds by deriving acceptances from Monte-Carlo simulations. However, this would be a considerable amount of work, in particular since most of those searches include finite-width effects. We, therefore, have not implemented those searches at this time, but strongly encourage the experimental collaborations to provide scalar acceptances -in a similar way as derived in Section 3.3.2 -for di-jet searches in the future.
Finally, ATLAS has also set a limit in the di-jet+ final state [29] that includes an interpretation as pp → tb(H ± → tb). We have implemented this limit. It is, however, strictly weaker than the dedicated charged Higgs search [50] in the same channel and using the same dataset.
All newly implemented searches in di-jet final states are summarized in Table 2. Figure 1: Feynman diagrams for di-jet production via a scalar resonance in association with initial-state photon radiation.

ATLAS search for di-jet final state in association with initialstate photon radiation
In this Section, we discuss the implementation of the ATLAS search for a BSM resonance decaying into a di-jet final state with initial-state photon radiation [28]. While the implementation of most other searches listed in Tables 1 and 2 is comparably straightforward, the implementation of the di-jet plus photon search requires a dedicated effort. 3 In Ref. [28] cross sections, acceptance, and efficiency values are only given for a simplified Z model, in which the Z couples to five lightest quarks with equal strength.
In order to make this experimental search usable within HiggsBounds (i.e. within the generic scalar model presented Section 2, see also the corresponding Feynman diagrams in Fig. 1), we need to derive the following dependencies of the cross section, the acceptance, and the efficiency: • Dependence on the different scalar-quark-quark couplings.
While in Ref. [28] the resonance with universal couplings to all quarks (except the top quark 4 ) was considered, we want to allow for a non-universal coupling structure and also take into account quark-flavor-violating couplings.
• Dependence on the electric charge of the resonance.
While in Ref. [28] the considered resonance was assumed to have no electric charge, we want to consider also the case of an electrically charged resonance allowing for the radiation of the photon from the resonance itself (see right diagram of Fig. 1).
• Dependence on the spin character.
While in Ref. [28] only the example of a spin-1 resonance is discussed, we want to implement the case of a spin-0 resonance.
• Dependence on the CP character of the decaying resonance.
While the spin-1 resonance considered in Ref. [28] was considered to have only CP-odd couplings, we want to allow for CP-even and CP-odd couplings of our spin-0 resonance.
In principle, a separate Monte Carlo simulation and a subsequent event analysis, reproducing the experimental analysis, has to be run for every parameter point and every possible quantum numbers of the resonance. This is not feasible for large parameter scans as they are often performed with HiggsBounds. While for the detector efficiency, we employ the values given in the auxiliary material of Ref. [28], 5 , more work is needed to work out fit functions for the cross section and the acceptance encoding the dependencies on all relevant parameters.

Fitting the cross section
For evaluating the cross section for di-jet production in association with a photon, we employ MadGraph5_aMC@NLO 2.7.0 [51] using the NNPDF3.0 LO PDF set [52]. As model file, we employ the UFO version [36] of the models presented in Section 2, which we generated using FeynRules.
Using this setup, we obtain the cross section for scalar production in association with a photon at leading oder (LO) 6 as a function of the resonance mass and its couplings to quarks. We calculate this cross section for multiple parameter points deriving simple fit formulas (listed in Appendix B) which we then implement into HiggsBounds. Within HiggsBounds, this cross section is then multiplied by the branching ratio of the scalar resonance to a di-jet final state.
In Fig. 2, we show exemplary cross-section results. 7 In the left panel, the cross-section for quark-initiated neutral scalar production in association with a photon is shown in dependence of the resonance mass. The scalar particle is assumed to have only one non-zero coupling to quarks. The non-zero coupling, which we set to one, is indicated in the plot legend. The cross section is highest if the quarks to which the resonance couples are abundant in the initial-state protons -i.e., if no PDF suppression occurs -and if up-type quarks, which have a higher electrical charge, are involved. Consequently, the cross section is highest if the resonance couples to up quarks (blue) reaching values of up to ∼ 60 pb for M S = 225 GeV, closely followed by the cross section for a resonance coupled to up and charm quarks (purple). All other cross sections are significantly lower. E.g., the cross section for a resonance coupled to down quarks (orange) is reduced by a factor of ∼ 5 in comparison to the up-quark-induced and up-charm-quark-induced processes.
The behavior is quite similar for the production of a charged scalar in association with a photon (see right panel of Fig. 2). The overall size of the cross sections is comparable to

Criterion
Single-photon trigger Combined trigger Number of jets n jets ≥ 2 Number of photons

Fitting the acceptance
To derive the analysis acceptance, we implement the analysis cuts of Ref. [28] into the MadAnalysis5 (MA5) framework [53][54][55][56]. 8 For the event generation, we employ the same setup as in Section 3.3.1. To emulate the detector effects, we smear the invariant mass distribution of the final-state jets by the detector mass resolution as given in the auxiliary material of Ref. [28]. The analysis has been performed using two distinct data sets recorded using two different triggers: • a single-photon trigger requiring at least one photon candidate with E γ T,trig > 140 GeV, and • a combined trigger requiring at least one photon candidate with E γ T,trig > 85 GeV and two jet candidates, each with p T > 50 GeV. Table 4: Cuts applied to distinguish inclusive and b-tagged signal regions.
The single-photon trigger is used for signal masses from 225 GeV to 450 GeV; the combined trigger from 450 GeV to 1100 GeV. The photons and jets used within the actual analysis are then defined using stricter requirements. The photons are required to fulfill • E γ T > 95 GeV, and • tight identification and tight isolation criteria [58][59][60].
Jets are reconstructed using the anti-k t algorithm [61] with radius parameter R = 0.4. After reconstruction, jets are required to fulfill • |η jet | < 2.8, • p jet T > 25 GeV, and • a jet within ∆R = 0.4 of a photon is removed.
Based upon these objects a series of selection cuts is applied to the events passing the trigger selection. These are summarized in Table 3. In addition, the selections based on the singlephoton and the combined trigger are further subdivided into an inclusive and a b-tagged signal region by applying the cuts displayed in Table 4.
We validated our analysis by comparing to the example numbers given in Ref. [28] for the simplified Z model. More details can be found in Appendix A.
For our present paper, it is, however, more appropriate to employ the limits for a generic Gaussian resonance given in Ref. [28] instead of the limits given for the simplified Z model. Correspondingly, we follow the procedure outlined in App. A of Ref. [62]: Instead of retaining the full tail of the m jj distribution (see Table 3), all events with a m jj value between 0.8M and 1.2M are selected where M is the signal mass. The mean mass for the events within this interval is then used to read off the limit from the results given in Ref. [28].
Using this analysis setup, we can derive the acceptance as a function of the resonance mass and the resonance's couplings to quarks. Since the simplified models as defined in Eqs. (1) and (2) involve a relatively high number of free parameters, it is quite costly to run the MC generation and analysis chain for a sufficient number of parameter points in order to derive a fit formula for the acceptance in dependence of all relevant parameters. To avoid this, we first checked whether the acceptance depends on the CP character of a neutral scalar resonance, or in case of a charged resonance on the chirality of its couplings to quarks. We find that these dependencies can be neglected (finding variations of 0.1%).
As a further simplification, we make use of the following observation: For almost all analysis cuts, the dependence on the flavors of the final-state quarks can be neglected. Therefore, the cut efficiencies can be approximated by a function depending only on the initial-state quark flavors. The only exceptions are the m jj cut as well as requirement of two b-tagged jets in the case of the b-tagged signal region. These cuts can be approximated by a function depending only on the flavors of the final-state quarks while the dependence on the flavors of the initial-state quarks can be neglected.
This approximate factorization can be written down in the form where A incl-SR is the acceptance of the inclusive signal region and A b-SR the acceptance of the b-tagged signal region. A initial is the cut efficiency of the initial cuts (excluding the m jj and the b-tagging cuts); A m jj , the cut efficiency of the m jj cut; and, A b-tag is the cut efficiency of the b-tagging cut. q 1 and q 2 are the initial-state quarks; q 3 and q 4 , the final-state quarks (see also Fig. 1). We cross-checked the approximation of Eqs. (3) and (4) for several example points (including cases with different initial and final state quark flavours) by directly evaluating the values for A incl-SR and A b-SR finding absolute deviations of 0.3%. We do not expect larger deviations for other parameter points.
The functions A incl-SR , A m jj , and A b-tag are then fitted by performing MC event generation and analysis for parameter points with different mass and coupling values. Since each of the functions depends only on two quark flavors -or equivalently on one scalar-quark-quark coupling (see Eqs. (1) and (2)) -far less parameter points have to evaluated in comparison to fitting the functions A incl-SR and A b-SR directly.
We display some of the resulting acceptance fits in Fig. 3. For this Figure, we assume that the resonance couples to only one type of quark pairs (as shown in the plot legend) corresponding to only one coupling in Eq. (1) being non-zero. Note, however, that our implementation also fully supports the situation in which the initial-state quark pair is different from the final-state quark pair. The vertical dashed lines at M S = 450 GeV mark the transition from the single-photon to the combined trigger. The markers denote points for which we produced MC samples. The uncertainty bands denote the corresponding statistical uncertainty. The colored lines represent the derived fit functions.
In the upper left panel of Fig. 3, the acceptances of a neutral scalar resonance in the inclusive signal region are shown in dependence of the resonance mass. While the acceptance is relatively low (∼ 1 − 2%) in the single-photon-trigger region, it increase to up to ∼ 8% in the combined-trigger region. As a consequence of the different PDFs for different quark flavors, the acceptance forūu-initiated neutral scalar production (blue curve) is up to ∼ 3% higher than e.g. the acceptance of thebb-initiated channel (gray curve).
In the upper right panel of Fig. 3, the acceptance of a neutral scalar resonance in the b-tagged signal region is shown. As a consequence of the b tagging, the acceptance is nonzero only for a bb final state. Due to the tighter cuts on the pseudo-rapidity of the jets, the acceptance is slightly reduced in comparison to the one of thebb-initiated channel in the upper left plot of Fig. 3.
The acceptances for a charged scalar resonance in the inclusive signal region are shown in the bottom panels of Fig. 3. While the overall behavior closely resembles the behavior of the acceptances for a neutral scalar resonance, the overall acceptance values are reduced by ∼ 1 − 2%. This a consequence of the cuts being optimized for a neutral resonance not taking into account the possibility of the photon being radiated by the resonance itself (see right panel of Fig. 1). Since the charged resonance can not decay to two bottom quarks, it does not contribute to the b-tagged signal region. Due to charge asymmetry of the proton's quark content, the acceptances for a negatively charged scalar (see bottom right panel of Fig. 3) are slightly smaller than the acceptances of for a positively charged scalar (see bottom left panel of Fig. 3). As already mentioned above, we require the mean invariant mass for all events with a m jj value between 0.8M and 1.2M , where M is the signal mass, as an additional input for applying the limits for a Gaussian resonance derived in Ref. [28]. Following the same procedure as for the acceptance, we fit this mean mass as a function of the initial-state quarks and the signal mass. The derived formulas are implemented into HiggsBounds and used to apply the Gaussian limits presented in Ref. [28]. For reference, we list them in Appendix B.
For a particle with different channels contributing, we add the cross-section values for the individual channels multiplied by the respective acceptances.

Model applications
In this Section, we discuss four model applications for the HiggsBounds extensions discussed in Section 3: constraints on a scalar particle mediating between SM and DM particles, constraints on BSM Higgs bosons with a non-SM like Yukawa coupling structure within a 2HDM framework, and constraints on the sneutrino and slepton sectors in the R-parity violating MSSM. We want to emphasize that the presented applications are exemplary applications and that the implemented searches can also be used to constrain other models.

Scalar Dark Matter portal model
As a first exemplary application for the extensions of HiggsBounds discussed in Section 3, we consider a scalar DM portal model, in which a neutral scalar mediates interactions between the SM particles and the DM particles. More concretely, we focus on a model which was used in Ref. [63] for the interpretation of mono-jet plus missing transverse energy searches. In this model, a pseudo-scalar mediator directly couples to all quarks (except the top quark) with an universal interaction strength g q . In addition, it couples to fermionic dark matter with the interaction strength g χ . For simplicity, we assume that the mediator does not couple to any other particles. We choose this model as an example in order to highlight the interplay of mono-jet searches -as performed in Ref. [63] -and di-jet searches -as discussed in Section 3. Since a recasting of the mono-jet search of Ref. [63] would be beyond the scope of the current paper, we choose the same DM mass as in Fig. 7a of Ref. [63] (i.e. 1 GeV).
This choice allows us to display the mono-jet limits of Ref. [63] and the di-jet limits as implemented in HiggsBounds in the same plot (see Fig. 4). In this plot, the collider constraints are shown in the (m S , g χ ) parameter plane, where m S is the mediator mass. The coupling of the mediator to quarks, g q , is set to 0.5. While the mediator is always produced 10 −1 10 0 Scalar DM portal model (g q = 0.5) Figure 4: Collider constraints on scalar dark matter portal model in the (m S , g χ ) parameter plane for g q = 0.5. The blue coloured areas indicate which di-jet search is expected to be most sensitive. The green colored area is excluded by the ATLAS mono-jet plus missing energy search of Ref. [63]. The area below the red curve is excluded by the ATLAS di-jet plus photon search [28]. The gray contours indicate the branching ratio of the scalar mediator S into invisible particles.
via a di-quark initial state controlled by the coupling g q , the relative strength of its decays into two quarks and two DM particles depends on the relative size of g χ and g q . The resulting branching ratio into DM particles is shown as gray contours in Fig. 4. While this branching ratio is below 10% for g χ 0.65, it lies above 90% for g χ 6.
Correspondingly, the mono-jet plus missing transverse momentum search of Ref. [63] limit (light green area) excludes the area of g χ 2 and m S 370 GeV. The sensitivity of this searches decreases rapidly for higher mediator masses. As a complementary constraint, the di-jet plus photon search of Ref. [28] excludes the area below the red curve -i.e., the region of g χ 3 for mediator masses between 225 GeV and 1100 GeV. The light (dark) blue regions indicate the area in which the inclusive (b-tagging) signal region of Ref. [28] is most sensitive. 9 The results shown in Fig. 4 demonstrate the complementary of mono-and di-jet searches for new resonances and the importance to provide a unified interpretation framework allowing to explore this complementarity.

Two-Higgs-Doublet model
if the Higgs at 125 GeV becomes SM-like. These effective couplings are independent of the fermion generation.
To demonstrate the impact of di-jet searches, we choose a simplified 2HDM parameter plane, where we only consider the heavy CP-even neutral scalar H. We assume that sin(β − α) = 1, which puts the light CP-even neutral scalar h into the exact alignment limit so that it can easily accommodate the SM-like properties measured for the Higgs at 125 GeV. We then vary tan β and m H through a large range and use 2HDMC [64] linked to HiggsBounds to obtain and test the model predictions for the different parameter values. 10 The results are shown in Fig. 5. For each parameter point, the color code indicates which limit is most sensitive to the heavy BSM scalar H. The analyses indicated in blue are previously implemented Higgs searches targeting the bb → H → bb channel [48,49,65], while the green analyses are the newly implemented di-b-jet [22,26] and di-b-jet+γ [28] searches. The red lines indicate the overall HiggsBounds 95 % C.L. upper limit on tan β. The solid line includes all Higgs and di-jet searches, while the dashed line only includes dedicated Higgs searches. Finally, the hatched region indicates where Γ tot H > m H /4, such that the narrow width approximation is no longer applicable.
Let us now consider the impact of the di-jet searches in Fig. 5 in detail. At m H < 200 GeV the low-mass di-b-jet search by CMS [26] leads to a huge difference in the excluded parameter region. The only dedicated Higgs search in this region was Ref. [65], which only covers m H ≥ 100 GeV. Below that mass, the only existing limits at large tan β came from LEP [66] and the Tevatron [67], leading to the weak exclusion of the dashed line. The newly implemented CMS di-b-jet limit, on the other hand, puts a very stringent upper bound on tan β, especially for very low m H 70 GeV, where it excludes the region down to tan β 6.
For intermediate mass values, the exclusion is dominated by the newly implemented dib-jet+γ limit [28]. For m H 500 GeV the limit is comparable to the dedicated CMS Higgs search [49]. Which of the two is the most sensitive changes based on the precise mass and -due to the width dependence of the bb + γ limit -tan β value. For larger masses, the bb + γ limit becomes more dominant, with a significantly stronger resulting limit on tan β. On the one-hand, this behavior is not unexpected, since the bb + γ analysis uses ∼ 80 fb −1 of data, while all other analyses use at most ∼ 36 fb −1 . However, especially for larger masses, the bb + γ search also benefits from the H + γ production cross section, which falls off less rapidly with increasing Higgs mass.
Out of the dedicated Higgs searches, Ref. [48] covers the highest masses, with a limit set up until m H = 1.4 TeV. The limit arising from the di-jet search Ref. [22] is substantially stronger already from m H 1.3 TeV. However, both limits only exclude parameter points in the hatched region, where the narrow width approximation is no longer applicable. Therefore, at very high masses even with the newly implemented searches no physically meaningful limit can be set.

R-parity violating MSSM
In the MSSM, the couplings of the scalar fermions are fixed by a small set of parameters. As a consequence of R parity, these bosons always decay into a final state containing missing transverse energy. If R parity is violated, however, the R parity violating couplings of scalar fermions are essentially free parameters. Within the superpotential of RPV SUSY one can find the terms, Thus, the couplings of the sneutrinos and the sleptons are controlled by the parameters λ ijk -coupling a left-handed slepton to two leptons -, and λ ijk -coupling a left-handed slepton to two quarks, where i, j, and k are generation indices. For a non-zero λ coupling, the respective sneutrino and slepton can be produced at the LHC via a di-quark initial state. Moreover, the presence of this coupling allows the decay into a di-jet final state. For a non-zero λ coupling, the sneutrino can decay into two charged leptons, whereas the slepton can decay to a neutrino and a charged lepton.
In Fig. 6, we display four exemplary parameter scans in which the two couplings λ ijk = −λ jik and λ ijk are varied. We choose four scenarios that illustrate different production channels and phenomenological signatures that could be covered by existing searches. The top left plot in Fig. 6 we performed a scan over the two couplings λ 323 = −λ 233 and λ 311 . This   means that the τ sneutrino, whose mass we choose to be at ∼ 1 TeV, 11 can be produced via a dd initial state and decays either to two jets or a muon and a tau lepton. The associated left-handed stau can be produced via a ud initial state and decays into two jets or the leptonic channel, µν τ . 12 Two of the newly implemented searches discussed in Section 3 are sensitive to this signature: the ATLAS pp → S → µτ resonance search [32] and the ATLAS di-jet plus photon search [28]. While the µτ search excludes the region of |λ 311 | 10 −2 and |λ 323 | 2·10 −2 , the di-jet plus photon search allows to set an upper limit on |λ 311 | of ∼ 0.6 showing nicely the complementarity between the different search channels. Note that here and in all the other scenarios, both the sneutrino and the -almost mass degenerate -associated slepton contribute to the di-jet rates. This is what makes the limits from di-jet final states competitive to the leptonic ones.
The top right plot of Fig. 6 shows the result of the scan of the two couplings λ 311 = −λ 131 and λ 321 . In this case the tau sneutrino can decay into a pair of jets or a pair of electrons while its production is caused by a sd (sd) initial state. The stau for this configuration is produced mainly byūd (ud) initial state and its decay consists of a pair of jets or the pair e + e − . In this case the two most significant searches are the ATLAS di-lepton resonance search [30] and the ATLAS di-jet plus photon search [28]. The highly sensitive e + e − search leads to stringent constraints down to |λ 321 | 3·10 −3 and |λ 311 | 5·10 −3 . As a consequence, the complementary di-jet plus photon search is only significant for large values of |λ 321 | and small values of |λ 311 |.
The third scenario is shown in the bottom left plot of Fig. 6, in which we scan over λ 121 = −λ 211 and λ 122 . In this case the lightest scalar fermions are the electron sneutrino and the selectron. The electron sneutrino can be produced in this case by ass initial state and decays into a pair of jets or a muon and an electron. In this case the selectron, that is produced by the cs (cs) initial state, decays into a pair of jets or the pair eµ. The two most sensitive searches in this scenario are the ATLAS di-jet plus photon search [28] and the CMS pp → S → eµ resonance search [33]. In this case the electron muon search can set the next bounds on the parameters |λ 122 | 1.5 · 10 −2 and |λ 121 | 3 · 10 −3 .
The last scenario (bottom right plot of Fig. 6) corresponds to a scan of the parameters λ 122 = −λ 211 and λ 133 . In this case the electron sneutrino decays into a pair of b-jets or two muons while it could be produced via abb pair. The selectron that decays into a pair of jets or the pair µν µ , can be produced by the initial state cb. In this scenario, the three most sensitive searches are: the newly implemented ATLAS two b-jets plus photon search [28] and ATLAS pp → S → µµ search [30], and the previously implemented CMSbb → S → µµ search [42]. The CMS muon search is able to constrain |λ 122 | 1.5 · 10 −1 while the ATLAS di-muon search constrains |λ 133 | 1.3 · 10 −1 . It is interesting to see here why both searches cover different areas in the (λ, λ ) parameter space. The CMSbb → S → µµ search [42] that is sensitive to small values of |λ 122 | sets bounds on resonances up to a mass of 1 TeV. Higher values of |λ 122 | make the physical mass of the sneutrino slightly greater than 1 TeV which makes this search insensitive. When we reach this threshold the ATLAS pp → S → µµ search [30], which covers higher resonance masses, becomes the most sensitive search. In this scenario we can see how the newly implemented searches are complementary to the previous ones.
As we have shown in the four cases appearing in Fig. 6 the newly implemented searches are quite sensitive and cover substantial parts of the parameter space of couplings. Our extension of HiggsBounds thus enables a deep exploration of the sneutrino/slepton sector of the RPV MSSM.

Conclusions
BSM scalars are a well-motivated target for LHC searches. Normally, new scalar states are anticipated to have coupling structures similar to the SM Higgs boson. Many examples for BSM theories containing new scalar states exists, however, for which this is not the case. In the present article, we presented an extension of the code HiggsBounds allowing to test such scalars with a general coupling structure.
Besides extending the HiggsBounds input routines, we also implemented ATLAS and CMS di-lepton and di-jet searches into HiggsBounds based upon generic scalar models. While the implementation of the di-lepton searches is straightforward, we have developed a new approach to implement limits with a non-trivial dependence on the involved couplings (as it is the case for the implemented di-jet limits): By recasting the analysis results and tabulating the acceptances and efficiencies as a function of all involved couplings in a generic scalar model, we are able to evaluate limits for arbitrary BSM scalar models without the need to run Monte-Carlo simulations for every parameter point. While this approach can be applied to existing searches by the means of recasting, we encourage the experimental collaborations to directly provide the relevant acceptance and efficiency functions in a sufficiently general simplified model, e.g. in our generic scalar models. This would significantly improve the applicability of the experimental results without the need to re-implement the analysis chain. Moreover, this would also allow for more precise limit setting, since it is normally only possible to recreate the experimental analysis chain in an approximate manner.
Evaluating the acceptances and efficiencies in our generic scalar models also allows us to judge the impact of some commonly used approximations when no applicable model interpretations are available and recasting every parameter point is not feasible. As expected, we have found that the acceptances for a scalar can differ significantly from those of a spin-1 boson, such that applying Z or W limits to scalar particles is an extremely crude approximation. Additionally, we observed that the acceptances for a scalar can vary significantly depending on the specific initial and final state, such that scalar interpretations assuming flavor-universal couplings only give a rough estimate of the true acceptances if flavor universality is not fulfilled. While we only investigated scalar particles in this study, this conclusion should be equally valid for spin-1 particles -e.g. in non-flavor-universal Z models -and we advise caution when generalizing flavor universal limits to such scenarios.
We used the extended HiggsBounds version for several example applications. First, we demonstrated the complementarity of di-jet and mono-jet limits for constraining simplified DM models. As a second example, we investigated the impact of generic di-b-jet resonance searches on the flipped 2HDM excluding parameter regions which are not yet constrained by dedicated heavy Higgs searches. As a third example, we discussed constraints on sneutrino and slepton couplings in the R-parity violating MSSM. While the newly implemented dilepton resonance searches set an upper limit on the product of sneutrino-lepton-lepton and sneutrino-quark-quark couplings, the di-jet limits allow to set an upper bound on the sleptonand sneutrino-quark-quark couplings regardless of the sneutrino-lepton-lepton couplings. The presented extension of HiggsBounds provides a useful tool for constraining BSM model containing scalars with an "exotic" coupling structure. All of the features described in this paper are available from HiggsBounds-5.10.0 onwards, which is available at https://gitlab.com/higgsbounds/higgsbounds .
In the future, we plan to apply the presented approach for implementing more complex analyses into HiggsBounds also to other search targets like mono-jet plus missing energy final states.
A ATLAS search for di-jet final state in association with initial-state photon radiation -validation In Section 3.3, we describe the implementation of the ATLAS search for a di-jet resonance produced in association with an initial-state photon into the MA5 framework. Here, we validate this implementation.
For the validation, we use the same model as employed in Ref. [28]. I.e., we generate MC events for a leptophobic Z resonance using the DMsimp_s_spin1 UFO model [38]. As in Ref. [28], the Z resonance is assumed to have only axial-vector couplings to quarks with a universal coupling strength g q (the coupling to top-quarks is, however, set to zero). All couplings to other particles are set to zero.
Using the same setup as described in Section 3.3, we generate MC samples consisting out of   Table 5: Search for a di-jet resonance produced in association with a photon: cut flows for BP1 comparing the MA5 event analysis as described in Section 3.3 and the numbers provided by ATLAS. For BP1, the single-photon trigger is employed.
For these two benchmark points, we can compare our analysis to the ATLAS cut flow tables which are given in the auxiliary material of Ref. [28]. The comparison of the cut flow tables can be found in Tables 5 and 6. This table shows the number of events remaining after each cut and the associated statistical uncertainty (corresponding to the uncertainty of the MC integration). We reconstructed the initial number of events by dividing the final number of events by the total acceptance values given in the auxiliary material of Ref. [28]. 13 The columns marked by denote the cut efficiencies.
The column denoted by δ shows the deviation between the cut efficiencies of Ref. [28] with respect to the MA5 implementation.
For BP1 employing the single-photon trigger (see Table 5), the cut efficiencies of our MA5 implementation are quite close to the numbers from Ref. [28]. Only the cut efficiency of the initial trigger and cleaning step as well as the b-tagging cut deviate by more than 10%. The numbers are, however, still well below the agreement level of under 30% recommended in Ref. [53]. Also note that the total acceptance values (before and after the b tagging), obtained by dividing the final number of events by the initial number of events, are very close with absolute deviations smaller than 0.1%.
The agreement level is very similar for BP2 employing the combined trigger (see Table 6). While the cut efficiency of the number of photons cut deviates by more than 10%, all deviations in the cut efficiencies are still well below 30%. Also the total acceptance values deviate by less than 0.4%.
As an additional validation, we compare the total acceptance values for different M Z masses given in Ref. [28] to our implementation. This comparison is shown in Fig. 7. In the left panel -showing the acceptances for the inclusive signal region -the agreement with the numbers given in Ref. [28] is very good for the single-photon trigger used for M Z ≤ 450 GeV. For heavier Z masses, for which the combined trigger is used, the agreement is slightly worse  Table 6: Search for a di-jet resonance produced in association with a photon: cut flows for BP1 comparing the MA5 event analysis as described in Section 3.3 and the numbers provided by ATLAS. For BP2, the combined trigger is employed. with a maximum deviation of ∼ 1%. For the b-tagging signal region, shown in the right panel of Fig. 7, a very similar behavior is visible. The overall size of the acceptances is, however, reduced by a factor ∼ 5. The maximum deviation between the numbers given in Ref. [28] and our implementation amounts to ∼ 0.2% for M Z = 950 GeV. Given this level of agreement, we consider our MA5 implementation as validated.

B Fit formulas cross sections and efficiencies
In this section we list the fit formulas used in the implementation of the di-jet + photon search [28].

B.1 Cross section fits
As described in Section 3.3.1, we evaluated the pp → φγ cross sections for each initial flavour separately and performed fits to this data. We use the functional form for the cross section fits, where g is the respective coupling as defined in Eqs. (1) and (2) and m S is the particle mass in GeV. The values of the fit coefficients for the different initial state flavours and final state charges are given in Table 7.

B.2 Acceptances and mass corrections
As described in Section 3.3.2, the acceptances and mass corrections are similarly fit to simulations for each initial or final state flavour and electric charge. The functional form used is For the preselection acceptance A initial this is a function of the particle mass m = m S and depends on the initial state flavours with the coefficients for the two trigger categories shown in Table 8. The final acceptance A m jj is a function of m = m S -the mean mass after the m jj cut -and depends on the final state flavours with the coefficients shown in Table 9. The additional acceptances of the b-tagged category are given by The mass transformation itself is parametrized as where the corresponding coefficients are given in Table 10.     C Applying HiggsBounds to exotic scalars: A User's Guide In this Appendix we provide a brief documentation of the new functions and subroutines that we added to HiggsBounds in HiggsBounds-5.10.0 to enable the application to neutral and charged scalars with an exotic, non-Higgs-boson-like coupling structure. We refer to Ref. [18] for the full documentation of the HiggsBounds-5 code.

C.1 Extension of the Theoretical Input Framework
We extended the standard HiggsBounds theoretical input framework both for the neutral and charged scalars. The cross sections for exotic production channels at the LHC at 13 TeV are parametrized as fit functions that take effective couplings as input, as described above.
For neutral scalars, we therefore introduce two new input subroutine for the neutral scalar couplings: First, the subroutine HiggsBounds_neutral_input_effC_firstgen, which takes the couplings to first-generation fermions as input, which are defined as SM-Higgs-bosonnormalized quantities, to be in accordance with the standard HiggsBounds effective coupling input: with y f = m f /v and v = 174.1 GeV. The up-quark, down-quark and electron masses are set to m u = 2.16 MeV, m d = 4.67 MeV and m e = 0.511 MeV, respectively. Second, the subroutine HiggsBounds_neutral_input_effC_FV takes the input for the flavor-violating (FV) couplings to quarks. 14 As these couplings vanish for the SM Higgs boson, they are not normalized and defined as g,g from Eq. (1).
In the standard HiggsBounds framework the branching ratios for neutral Higgs bosons can be approximated from the provided effective coupling through the rescaling of the corresponding partial width predictions from the LHC Higgs Cross Section Working Group (LHC HXSWG). As these predictions do neither exist for first generation fermionic final states, nor for flavor-violating fermionic final states, the corresponding branching ratios (BR) have to be set directly as input. The input for the branching ratios of neutral scalar boson two-body decays to first generation fermions and FV quark final states can be provided via the new subroutines HiggsBounds_neutral_input_firstgenBR and HiggsBounds_neutral_input_FVBR, respectively, see Table 11 for details.
For charged scalar bosons the effective couplings can be specified via the subroutine HiggsBounds_charged_input_effC_fermions. These are defined as g qL , g qR from Eq. (2). The branching ratios for charged scalar boson decays to fermionic final states ud, us, cd, ub, eν and µν can be set via the subroutine HiggsBounds_charged_input_firstgenBR.