Higgs, di-Higgs and tri-Higgs production via SUSY processes at the LHC with 14 TeV

We have systematically investigated the production of a Higgs boson with a mass of about 125 GeV in the decays of supersymmetric particles within the phenomenological MSSM (pMSSM). We find regions of parameter space that are consistent with all world data and that predict a sizeable rate of anomalous Higgs, di-Higgs and even tri-Higgs events at the 14 TeV LHC. All relevant SUSY production processes are investigated. We find that Higgs bosons can be produced in a large variety of SUSY processes, resulting in a large range of different detector signatures containing missing transverse momentum. Such Higgs events are outstanding signatures for new physics already for the early 14 TeV LHC data. SUSY processes are also important to interprete deviations found in upcoming Standard Model Higgs and di-Higgs production measurements.


Introduction
The Higgs-boson discovery at the Large Hadron Collider (LHC) [1,2] marks the beginning of a new era in particle physics. It gives us exciting new possibilities to study the physics of the Standard Model (SM) of particle physics. In this paper we investigate the next level of Higgs-boson searches, namely the possibility that Higgs bosons with a mass of about 125 GeV are produced by processes involving physics beyond the SM.
Supersymmetry (SUSY) [3][4][5][6][7][8][9][10][11][12][13][14][15][16] is one of the conceivable extensions of the SM. It could provide a natural candidate for cold dark matter if R-parity is conserved [17,18] and it allows for a stabilization of the electroweak scale by reducing the fine tuning of higherorder corrections to the Higgs mass [16,[19][20][21][22][23]. In its minimal version, i.e. the Minimal Supersymmetric Standard Model (MSSM), SUSY predicts superpartners for the existing JHEP05(2015)044 SM particles and two Higgs doublets instead of one. On top of that, R-parity is assumed to be conserved in the MSSM, which results in the existence of a lightest supersymmetric particle (LSP). If the LSP is a neutralino, i.e. a Majorana-fermion superpartner associated with the neutral SM bosons in the electroweak sector, it is only weakly interacting and stable. It escapes detection, which results in missing transverse momentum in the detector.
In the present study we investigate systematically the possibilities to produce Higgs bosons with a mass m h 0 ≈ 125 GeV 1 in the decay of SUSY particles. This analysis is based on the phenomenological MSSM (pMSSM) [24,25]. The pMSSM is scanned for parameter regions where the SUSY particles have a viable branching ratio to Higgs bosons. Only those models are selected that fulfil the current constraints on SUSY. The relevance for the upcoming LHC runs at 14 TeV is discussed in detail and the most relevant Higgs production processes are identified. Higgs production via particular SUSY processes has been studied e.g. in [26][27][28][29]. We calculate the allowed production rates for anomalous Higgs, di-Higgs and tri-Higgs events. Subsequently, LHC events are simulated for each interesting model. These events are classified into topologies according to the SM particles produced in association with the Higgs boson(s) and the Higgs kinematics is studied. We identify topologies that are interesting for extending the current SUSY searches. Experimentally the events might be best detectable by explicitly "tagging" the Higgs boson(s) in SUSY searches. Since the invariant mass of the (lightest) Higgs boson is known and well reconstructable in many decay modes, and since we know that the SM rate to produce Higgs events with large missing transverse momentum (and maybe other SM particles) is small, a "Higgs-tag" can provide a unique signature for new physics. A special "Higgs-tag" for a boosted Higgs has also been suggested [30][31][32][33].
A few analyses have already searched for such events in ATLAS and CMS data. Higgs production via χ 0 2 χ ± 1 neutralino-chargino production has been investigated in ATLAS [34] and CMS [35,36]. In addition, searches have been pursued by CMS for a simplified model with a Higgs produced in top squark decays [37,38]. The present study aims to systematically investigate the possibility to produce Higgs bosons within the current constraints on SUSY by considering all relevant SUSY processes and decays. This paper is organized as follows. In section 2 the most important supersymmetric decay mechanisms for producing light Higgs bosons are discussed. In section 3 the pMSSM parameter space is scanned for models that are consistent with all current experimental constraints on SUSY and that have the potential to produce sizeable Higgs-boson event rates. Finally, in section 4 the surviving pMSSM models are studied with regard to the expected Higgs-boson event rates at the early stages of the upcoming LHC run and with regard to special kinematical features, such as boosts and missing transverse momentum.

Supersymmetric decays into the lightest Higgs
In view of its important role in producing Higgs bosons, we start with a detailed discussion of the neutralino/chargino sector. In the MSSM some of the superfields mix as a result of SUSY breaking to form new mass eigenstates. Let's first consider the neutral SM bosons in JHEP05(2015)044 the electroweak sector, i.e. the hypercharge B boson, neutral weak W 3 boson and neutral components of the two Higgs doublets. The associated Majorana-fermion superpartners, i.e. the Bino B, neutral Wino W 3 and neutral Higgsinos H 0 d and H 0 u , mix to form neutral mass eigenstates called neutralinos ( χ 0 1,2,3,4 , numbered in increasing mass order). This mixing is caused by off-diagonal terms in the neutralino mass matrix, which acts on the Bino, Wino and Higgsino fields [39]: Here s α ≡ sin α and c α ≡ cos α. The parameters M 1 and M 2 are the SUSY-breaking mass parameters for the Bino and Winos, µ is the SUSY version of the SM Higgs-mass parameter, cos θ W = m W /m Z is the ratio of the SM W -boson and Z-boson masses, and tan β is the ratio of the two Higgs vacuum expectation values. A similar mixing phenomenon occurs in the charged sector, belonging to the charged weak bosons W ± and the charged components of the Higgs doublets. The associated Dirac-fermion superpartners, i.e. the charged Winos W ± and Higgsinos H ± u/d , mix to form charged mass eigenstates called charginos ( χ ± 1,2 , numbered in increasing mass order) as a result of the mixing in the chargino mass matrix [39]: The mixing in the neutralino and chargino mass matrices stems from terms that go with the Z-boson mass. However, in the case that M 1 , M 2 and |µ| largely exceed the mass of the Z-boson, the mixing terms are relatively small. If we neglect the mixing terms, the neutralinos are either a Bino, a Wino or a symmetric/antisymmetric mix of both Higgsino The charginos are in that case either a charged Wino or a charged Higgsino. The composition for all possible regimes is shown in table 1. In this simplified case the mass of the Bino neutralino is M 1 , the masses of the Wino neutralino and charginos are M 2 , and the masses of the Higgsino neutralinos and charginos are |µ|.
In fact some of the eigenvalues of the mass matrices will turn out to be negative. For instance, H 0 S corresponds to the eigenvalue −µ, whereas H 0 A corresponds to the oppositesign eigenvalue +µ. In order to arrive at a proper (non-negative) definition of the mass of JHEP05(2015)044

Regime
Composition neutralinos Composition charginos all particles, an extra factor γ 5 will have to be absorbed into the definition of the negativemass eigenstates, which flips the sign of the corresponding mass eigenvalue. As we will see, this extra factor γ 5 has important consequences for the decay properties of the neutralinos. When we switch on the mixing again, mixed neutralino states consisting of Binos, Winos and Higgsinos will exist. However, since the mixing is small, there will always be a part that dominates the state, which we then refer to as Binolike, Winolike or Higgsinolike. The true masses of all the neutralinos and charginos behave as in the previously discussed simplified case, which is governed by the three mass parameters M 1 , M 2 and µ.

Neutralino and chargino decays into the lightest Higgs
If we choose the lightest neutralino to be the LSP, all supersymmetric particles will eventually decay into a lightest neutralino. The branching ratios of the most important direct decay channels of neutralinos into the lightest Higgs boson h 0 accompanied by a LSP are shown in figure 1. The lightest chargino plays an important role if it is of almost the same mass as the lightest neutralino. Therefore the branching ratio of figure 1d is also included.
Some of the features of these decay processes can be explained very well kinematically with the previously discussed simplified case. For example, the decay χ 0 2 → χ 0 1 + h 0 is very unlikely in the case that M 1 , M 2 > |µ| or when the smallest two parameters of the set M 1 , M 2 , |µ| are relatively close (i.e. less than m h 0 apart), as can be seen in figure 1a. This is because both neutralinos have more or less the same mass in that case, which means that the decay χ 0 2 → χ 0 1 + h 0 is kinematically not allowed. For the same reason the decay χ ± 2 → χ ± 1 + h 0 is greatly suppressed in the region around M 2 ≈ |µ|, as can be seen in figure 1d. In figure 1b we see that a similar thing holds for the decay χ 0 3 → χ 0 1 + h 0 for M 2 ≈ |µ| < M 1 or M 1 ≈ |µ| < M 2 , since in that case the lightest three neutralinos have more or less the same mass.
For some of the features of these decay processes, such as the apparent complementarity of BR( χ 0 2 → χ 0 1 + h 0 ) and BR( χ 0 3 → χ 0 1 + h 0 ) for M 2 > |µ| > M 1 , we have to dig a little bit deeper. In order to facilitate the discussion we first list in table 2 the possible interactions between the Binos, Winos and Higgsinos, from which the neutralinos and charginos inherit their decay properties. In order to identify the interactions that involve the light Higgs boson, the two Higgs doublets are represented by the associated five Higgs mass eigenstates JHEP05(2015)044  Figure 1. Prominent direct branching ratios for the decay of neutralinos/charginos into the lightest Higgs boson and the lightest neutralino/chargino in the case that µ = 500 GeV, tan β = 50 and all other parameters scaled up to very high values.

JHEP05(2015)044
Besides the kinematical observations mentioned before, we observe the following features in the various mass domains. These features mainly involve the competition between the decay modes into the h 0 and alternative decay modes involving W or Z bosons.
The Binolike-Winolike-Higgsinolike mass domain M 1 < M 2 < |µ|: • The decay χ 0 2 → χ 0 1 + h 0 tends to dominate the branching ratio of χ 0 2 if it is kinematically allowed, resulting in values for the branching ratio BR( χ 0 2 → χ 0 1 + h 0 ) that can get close to unity (see figure 1a). As can be read off from table 2, this is caused by the fact that the decay into the h 0 involves the (suppressed) Higgsino component of either the Binolike χ 0 1 or Winolike χ 0 2 , whereas the decay into a Z boson involves the (double suppressed) Higgsino components of both these neutralinos. At the same time the decay χ 0 can only reach values that are substantially smaller (see figures 1b and 1c), since the alternative decay modes χ 0 3,4 → χ ± 1 + W ∓ and χ 0 3,4 → χ 0 2 + h 0 /Z seriously reduce the maximum branching ratio for direct decays into the LSP. Note, however, that the χ 0 3,4 neutralinos might be of interest for di-Higgs decay modes in view of the possible two-step decays χ 0 3,4 → χ 0 2 + h 0 followed by χ 0 The Winolike-Binolike-Higgsinolike mass domain M 2 < M 1 < |µ|: the same arguments in principle apply to this mass domain. However, in this case the (single suppressed) decay mode χ 0 2 → χ ± 1 + W ∓ cannot be avoided since m χ ± 1 ≈ m χ 0 1 . As a result BR( χ 0 2 → χ 0 1 + h 0 ) will at most reach 0.3 in this mass regime.
The Binolike-Higgsinolike-Winolike mass domain M 1 < |µ| < M 2 : • We observe both very large (almost unity) and very small branching ratios for the decay into the h 0 , with as additional striking feature the apparent complementarity of BR( χ 0 2 → χ 0 1 + h 0 ) and BR( χ 0 3 → χ 0 1 + h 0 ) (see figures 1a and 1b). This has to do with the occurrence of negative-mass eigenstates in the Higgsino sector and the associated factor γ 5 that is introduced in order to flip the sign of the mass eigenvalue. If χ 0 2 corresponds to a genuine positive-mass eigentstate, then χ 0 2 → χ 0 1 + h 0 is an unsuppressed scalar decay mode (see table 2) that tends to dominate the single suppressed decay into a Z boson. If χ 0 2 corresponds to a negative-mass eigenstate, then χ 0 2 → χ 0 1 + h 0 is a velocity-suppressed pseudoscalar decay mode and this time the decay mode into a Z boson dominates. The observed complementarity follows from the fact that χ 0 2 and χ 0 3 correspond to opposite-sign mass eigenvalues, while at the same time the suppression factors are such that the role of h 0 and Z are effectively interchanged in the two cases [40].

JHEP05(2015)044
The Winolike-Higgsinolike-Binolike mass domain M 2 < |µ| < M 1 : the previous arguments in principle apply to this mass domain as well. However, in this case the decay modes χ 0 2,3 → χ ± 1 + W ∓ cannot be avoided. This reduces the maximum combined branching ratio for the decays into h 0 and Z to roughly 0.3. For this maximum combined branching ratio again a complementarity phenomenon is observed (see figures 1a and 1b).
The Higgsinolike LSP mass domain |µ| < M 1,2 : • As mentioned before, the decay χ 0 2 → χ 0 1 + h 0 is not allowed kinematically since both neutralinos have more or less the same mass in the Higgsinolike LSP case.
The chargino decays: • The lightest chargino can never decay directly into the lightest Higgs boson. This is due to the fact that R-parity conservation forbids a decay to a neutral Higgs boson. For this a lighter charged supersymmetric fermion is needed, which is not present in the case of the lightest chargino. It is possible for a lightest chargino to decay into three particles, but such three-particle chargino decay modes with a Higgs boson featuring in the final state are rather rare.
• For M 2 < |µ| the branching ratio BR( χ ± 2 → χ ± 1 +h 0 ) can reach maximum values of up to 0.35 as a result of the competition from the unavoidable decay mode χ ± 2 → χ 0 1 +W ± as well as the decay mode χ ± 2 → χ ± 1 + Z. The total branching ratio for multi-step decays into the lightest Higgs boson can substantially exceed 0.35 in view of the possibility of two-step decays of the form χ ± 2 → χ 0 2 + W ± followed by χ 0 • For |µ| < M 2 the branching ratio BR( χ ± 2 → χ ± 1 + h 0 ) gets additionally reduced by the alternative decay mode to the second Higgsinolike neutralino. Again multi-step decay modes can substantially enhance the branching ratio for the decay into the h 0 .
In conclusion, the branching ratios for direct decays of neutralinos/charginos into the LSP and the lightest Higgs boson can be pretty large, reaching maximum values close to one for χ 0 2,3 . For χ ± 1 there is effectively no decay into the lightest Higgs boson. For the heavier states χ 0 4 and χ ± 2 the direct-decay branching ratios can reach 0.35 at best. However, for these heavy SUSY particles the total branching ratio for multi-step decays into the LSP and the lightest Higgs boson can be substantially larger if the non-Higgs decay step gives rise to χ 0 2,3 , which can subsequently decay into the lightest Higgs boson with high probability.

Sfermion decays into the lightest Higgs
Next we give a brief summary of the other supersymmetric decay channels that can produce a lightest Higgs boson, starting with the sfermions (squarks and sleptons), the scalar superpartners of the SM fermions (quarks and leptons). Such decay modes will play a role later on when the masses of the sfermions are not artificially scaled up to very high values. In this context it should also be noted that, apart from the interactions listed in table 2, the Binos, Winos and Higgsinos can also decay into fermion-sfermion pairs, involving the Yukawa interactions.
Since the Wino couples only to left-handed sfermions, the decays of left-handed ( f L ) and right-handed ( f R ) sfermions are different. In addition, the couplings to Higgsinos are Yukawa suppressed. This results in a profound difference between the decays of 1st/2nd generation sfermions and 3rd generation sfermions, since only the latter may have a large coupling to the Higgsinos.
Direct decays of sfermions into the lightest Higgs boson: • First of all there is the possibility to have a mass difference between left-and righthanded sfermions. As a result, there is the possibility for f L,R → f R,L + h 0 decay modes if the mass difference between the left-and right-handed sfermions exceeds 125 GeV. The couplings involved in this decay mode are Yukawa suppressed in the pMSSM. Therefore, the direct decay is mostly relevant for 3rd generation sfermions.
• The sfermions of the 3rd generation are mixtures of left-and right-handed states, indicated by f 1,2 (numbered in increasing mass order). Therefore, there is an additional possibility for a direct decay via f 2 → f 1 + h 0 . This decay can involve a non-Yukawa-suppressed (gauge) coupling between two left-or two right-handed components of the sfermion mass eigenstates. For 3rd generation squarks the gauge coupling and the Yukawa coupling can be of the same order of magnitude, which can lead to unexpected cancellations between both direct h 0 production mechanisms in that case.
Indirect decays of 1st/2nd generation sfermions into the lightest Higgs boson: sfermions can decay to heavy neutralinos or the heavy chargino, which can subsequently decay into lighter neutralinos or charginos and the lightest Higgs boson. The decay pattern differs for the left-and right-handed sfermions, depending on the composition of the LSP.
• Winolike LSP: the direct decay of the right-handed sfermions to the LSP is suppressed. If kinematically allowed the right-handed sfermions will decay to the Binolike neutralino, with the decay to the Higgsinolike states being Yukawa suppressed. As explained above, this Binolike neutralino can decay with a moderately large branching ratio to the h 0 , since the decay to the Z boson is double suppressed. The left-handed sfermions predominantly decay to the LSP, which strongly reduces indirect decays into the lightest Higgs boson.
• Higgsinolike LSP: the decay of the right and left-handed 1st/2nd generation sfermions to the LSP is Yukawa suppressed. If possible, these sfermions will decay to the heavier

JHEP05(2015)044
Bino-or Winolike states. As explained above, these states can decay with reduced branching fraction to the h 0 .
• Binolike LSP: if the χ 0 1 is Binolike, the right-handed sfermions predominantly decay to the LSP. However, the left-handed sfermions still prefer to decay (if kinematically allowed) to the heavier Winolike neutralino/chargino. This is caused by an intrinsic c θ W /s θ W ≈ 1.9 enhancement factor of the weak coupling of sleptons compared to the hypercharge coupling, with an additional factor 3 enhancement for squarks. As explained above, in these models the Winolike neutralino can have a large branching ratio to Higgs bosons. If it is a χ 0 4 , then decays to charginos are also possible. If it is a χ 0 2 , then its branching ratio to h 0 bosons is potentially very large and can be close to unity.
For 3rd generation sfermions the couplings to the Higgsinolike states are not Yukawa suppressed anymore and can even become large for top squarks. This results in a richer structure of possible decay modes and a more prominent role of Higgsinolike states as decay products. In that case indirect Higgs production can also become important in scenarios where Higgsinolike states have a large branching ratio into the lightest Higgs boson, as described in the previous subsection.

Heavy Higgs-boson decays into the lightest Higgs
Also the heavy Higgs particles can decay into the h 0 . These particles are a consequence of SUSY, which requires more than one Higgs doublet, but as far as R-parity is concerned they qualify as "SM" particles. Consequently, these particles do not necessarily have to decay into the LSP and therefore do not necessarily give rise to large missing transverse momentum in their decay chains. A comprehensive overview of the decays of the heavy Higgs bosons is given in [41].
• Direct decays: as will be discussed later, all surviving MSSM models have M A values exceeding 300 GeV. This is known as the decoupling limit (large M A ) and consequently all heavy Higgs bosons have similar masses, which blocks decays among heavy Higgs bosons. The heavy CP-even Higgs boson H 0 can directly decay to two h 0 bosons. The CP-odd Higgs boson A 0 can decay to h 0 Z. The charged Higgs bosons H ± can decay to W ± h 0 . The corresponding branching ratios tend to be rather small, because in most surviving MSSM models tan β is relatively large ( > 10) and consequently the decays of the heavy Higgs bosons to b-quarks are dominant. Later on we will encounter a couple of exceptional models that have the lowest values for M A and at the same time a relatively small value for tan β in order to survive the experimental constraints. Such models have noticeable branching ratios for direct heavy Higgs-boson decays into the lightest Higgs. More details can be found in ref. [41].
• Indirect decays: the heavy Higgs bosons also have the possibility to decay into heavy neutralinos or charginos (if kinematically allowed), especially if one of these decay states is Higgsinolike and the other Bino-or Winolike (see table 2). Those states can subsequently decay to h 0 , sometimes with high branching ratios. This can, for instance, result in di-Higgs production from an A 0 decay.

JHEP05(2015)044
3 Finding candidate pMSSM models: simulation and constraints The MSSM has more than 100 free parameters. Most of those parameters are not relevant for LHC physics. In the pMSSM the free parameters are reduced to 19 by demanding CP-conservation, minimal flavour violation and degenerate mass spectra for the 1st and 2nd generations of sfermions. The LSP is required to be the neutralino χ 0 1 in order to have a viable dark-matter candidate. This reduced model should cover a large fraction of the relevant SUSY phase space for h 0 production. The 19 remaining parameters are 10 sfermion masses, 2 3 gaugino masses M 1,2,3 , the ratio of the Higgs vacuum expectation values tan β, the Higgsino mixing parameter µ, the mass m A of the CP-odd Higgs-boson A 0 and 3 trilinear scalar couplings A b,t,τ .

Generation and pre-selection of pMSSM model-sets
SUSY-HIT [42] is used to generate the particle spectra of the 19-parameter pMSSM models. Only models are selected with a neutralino as LSP. The Higgs mass has been precisely determined by ATLAS and CMS to be 125.4 (ATLAS [43]) and 125.0 GeV (CMS [44]) with uncertainties of 0.3 − 0.4 GeV for each experiment. We select only models with a lightest Higgs boson h 0 within the range: In addition we produce two statistically independent sets of models: • Set A: Higgs production via direct decay of an arbitrary SUSY particle or a heavy Higgs boson. As described in the previous section, Higgs production can occur via various different decays of SUSY particles. In addition, h 0 bosons can be produced in the decay of heavy Higgs bosons. For this set we require in the preselection that at least one SUSY particle or heavy Higgs boson has a direct branching ratio to h 0 exceeding 20%.
• Set B: Higgs production via direct decays of charginos or neutralinos. Since Higgs production via neutralino or chargino decays is most important, a second set of models dedicated to these decays is produced. For this set we required that at least one of the following direct branching ratios exceeds 20%: The advantage of set B is that less model points are needed to study the most relevant Higgs production modes, since Higgs production predominantly originates from the decay of a heavy neutralino or chargino. Those neutralinos and charginos can be directly produced JHEP05(2015)044 or they are produced in cascade decays of predominantly squarks or gluinos, since these coloured SUSY particles can have a large cross section. The advantage of set A is its larger coverage of possible 3rd-generation and heavy-Higgs decay modes.

Parameter space coverage with a particle filter
This study has not the objective to provide a statistical interpretation like a "coverage" or a "likelihood" for a given parameter region. The objective is to find regions in the parameter space that are consistent with the global constraints on SUSY and where in addition the production of h 0 bosons is large (or close to maximal) in order to determine possible rates and topologies for SUSY Higgs production at the LHC. Each of our parameter sets represents a viable model point that could be realized in nature.
We use a simplified two-step particle filter algorithm [45] to find model points in the pMSSM parameter space.
1. First the 19 parameters of the pMSSM (3 gaugino masses, 6 squark masses, 4 slepton masses, 3 trilineair couplings, M A , µ and tan β) are chosen randomly from a flat prior distribution. The squark and slepton masses and M A have ranges between 100 GeV and 3000 GeV. The Higgsino mixing term µ, which in principle can be negative, ranges between -3000 GeV and 3000 GeV. This is also the case for the trilineair couplings, although we choose the ranges in that case between -5000 GeV and 5000 GeV to be sure that the trilinear couplings will not restrict the simulation too much. The lower bound on the gaugino masses is chosen to be 10 GeV to ensure that neutralinos, charginos and gluinos with very low masses are also evaluated. Finally, the ratio tan β of the Higgs vacuum expectation values is chosen between 1 and 50. For each set of pMSSM parameters SUSY-HIT [42] is used to generate the SUSY particle spectra and mixing matrices. Subsequently the preselection criteria of the previous subsection are checked. Model-sets are generated randomly within the given parameter range until we find 10000 model-sets fulfilling the preselection requirements.
2. These model-sets are then used as seeds (or particles) to build a posterior probability distribution from which further model-sets are generated. The posterior probability distribution is chosen as a sum of multi-dimensional Gaussian distributions centered around the parameter values S of each seed point. The multi-dimensional width of the Gaussian distributions is set to 10%, 25% and 40% of L d , where L d is the extent of the parameter space in dimension d out of 19. Around each seed further models are generated. The three sets are compared in order to evaluate the dependence on the width of the sampling. A comparison of the width dependence and a comparison of sets A and B is shown in figure 18 in the appendix. Since no significant difference is found all sets have been merged. This simulation process continues until in total at least 250000 models survive all preselection criteria.

Experimental constraints
The code micrOMEGAs [46] is used to calculate specific observables for each model-set in order to compare them with the experimental restrictions. The following constraints

JHEP05(2015)044
impact especially the neutralino and chargino mixing and can also influence their decay to the lightest Higgs boson.
• From the WMAP and the Planck data we adopt the cold dark-matter (DM) relic density in the universe [47,48]. We select a region corresponding to the last Planck central value 0.1186 ± 0.0031 including an 10% (upper) and 20% (lower) theoretical uncertainty: • The limits from the 85.3 days Large Underground Xenon (LUX) data [49] are taken into account. To compare the calculated proton/neutron cross sections σ (p/n) to the experimental limits, we use a normalized cross section for a point-like nucleus [50]: with A and Z the mass number and atomic number of the target. In our case the target is xenon with A = 131 and Z = 54.
• We implement the LHCb and CMS measurements of the branching ratio of the strange B meson to two muons [51,52] by demanding (3.5) • We impose the LEP limits on the invisible width of the Z boson and on the SUSY particle masses as implemented in micrOmegas.
The WMAP/Planck results place severe constraints on the models as can be seen in figure 2. The LSP's of the surviving models turn out to be mostly Binolike, with a relatively low mass, and to a lesser extent Higgsinolike or Winolike, with a relatively high mass. This is caused by the possibility of coannihilation of the LSP together with the lightest chargino or next-to-lightest neutralino, which is mostly absent for Binolike LSP's. In order to reduce the efficiency of the coannihilation and have a DM relic density that is not too low, Higgsinolike and Winolike LSP's are substantially heavier than Binolike LSP's. The occurrence of Winolike LSP's is suppressed within the simulated parameter space, since in that case the coannihilation turns out to be very efficient. Among the useful models that survive the WMAP/Planck constraint we have found only a few with a Winolike LSP.
Having checked the impact of the WMAP/Planck constraint, we now impose the LUX limits on the surviving models. The additional impact of the LUX limits is much smaller, as can be seen in figure 3 where the LUX experimental limits are imposed. Given the WMAP/Planck and LUX constraints, the B s and LEP constraints have little additional impact on the number of viable models. Notable exceptions are the surviving models in figure 3 with a very light Higgsinolike LSP, which are removed by the LEP constraints on the lightest chargino mass.

JHEP05(2015)044
After a first iteration about 250 models survived. These models were again used to construct a posterior for a second particle-filter iteration, producing in total about 430 models that fulfil all constraints. As expected the success rate with this posterior is increased.

ATLAS constraints: event generation, fast simulation and analysis
The remaining model-sets are compared with recent constraints from LHC SUSY searches at 7 TeV and 8 TeV centre-of-mass energies. Most important are constraints from searches for squarks and/or gluinos and chargino-neutralino production. Searches for heavy Higgs production had no influence on the remaining models. The mass constraint on the h 0 boson demands either large M A or very heavy top squarks. In fact, the model with the lightest A 0 boson in our sample has M A ≈ 330 GeV. Since tan β = 6.9 is relatively small for this model, it is not excluded by the ATLAS and CMS searches for heavy Higgs production [53,54].
The limits of the ATLAS experiment on light squarks, gluinos and chargino-neutralino production are implemented by emulating the ATLAS analysis chain. Events from LHC collisions are generated for each pMSSM model and the detector response is simulated by a fast detector simulation. The acceptance and efficiency is determined by applying the most important ATLAS analysis cuts on the simulated events. Finally, these numbers are used to calculate the expected number of signal events for each signal region and analysis. Subsequently, these expected yields are compared to the model-independent 95% C.L. limits provided by ATLAS. PYTHIA 6.4 [55] is used for the event simulation of proton-proton (pp) collisions at a 7 TeV and 8 TeV centre-of-mass energy. All SUSY production processes are enabled. For every model point and each centre-of-mass energy 10000 events are generated, which we found to be enough even for the models with the smallest selection efficiencies. To get as close as possible to the ATLAS analysis we use DELPHES 3.0 [56] as a fast detector simulation with the default ATLAS detector card, modified by setting the jet cone radius to 0.4. The PYTHIA output is read in by DELPHES in HepMC format, which is produced by HepMC 2.04.02 [57]. The object reconstruction is done by DELPHES, which uses the same anti-k T jet algorithm [58] as ATLAS. Also included in the reconstruction are isolation criteria for electrons and muons. We do not emulate pile-up events.
The 7 TeV analysis implementation is identical to ref. [59]. The selection efficiencies of our own implementation were compared to ATLAS in ref. [59] and were found to agree within uncertainties. For this study the implementation used in ref. [59] was updated using the recent 8 TeV selection criteria. For the chargino-neutralino searches the SR0τ a selection with all 20 bins was implemented as described in ref. [60]. For the squark and gluino searches all 13 signal regions without explicit W selection of ref. [61] are considered. In order to check constraints from multi-b-jet searches we included also signal region SR-0 -A from [62].
Preliminary direct searches for decays into h 0 bosons from neutralinos do not influence the remaining models. The mass of the lightest neutralino with a sizeable branching ratio to Higgs bosons is about 185 GeV and the mass of the LSP is at least 40 GeV. This is well beyond the exclusion reach of the ATLAS and CMS searches in these channels [34][35][36]. After the event selection, the event counts are scaled to the luminosities considered in the analyses with leading order cross sections as provided by Pythia. The limits on the effective cross sections given by the ATLAS analyses are used to calculate a limit on the number of signal events passing the cuts. No attempt was made to include theoretical uncertainties. In the studied SUSY mass range these uncertainties are small compared to the differences of the ATLAS and DELPHES setups and would not change drastically any conclusion of this work.

JHEP05(2015)044
In the end, 252 of the models passed all selection criteria. Figure 4 shows the excluded and non-excluded models as a function of the gluino mass and the minimal mass of the first and second generation squarks m miñ q . Most excluded model points are due to limits on squarks and gluinos and have squark or gluino masses below about 1500 GeV, in agreement with current LHC limits. All models with a gluino mass below 750 GeV are excluded. Remarkably, a large fraction of models with low squark masses is still allowed. One wellknown reason for this is that the lightest squark can be compressed with the χ 0 1 as shown in figure 5. This leads to very soft jets from squark decays. The squarks might only be visible via mono-jet signatures.
The enhancement of Higgs production in the studied models leads to a second interesting feature that causes the fraction of non-excluded models in this study to be larger than previously found in other scans (e.g. in [63]). In many non-excluded models the lightest squarks are compressed with a heavy neutralino/chargino. To illustrate this we indicate the minimal mass of all first and second generation squarks and the gluino by A. Figure 6 shows the smallest difference min(∆A) between A and the masses of all neutralinos and charginos as a function of A (given that the neutralino or chargino mass is smaller than A). In contrast to figure 5, all non-excluded models with A < 800 GeV have a mass difference ∆A below 300 GeV, which implies that many squarks are compressed with χ 0 2,3,4 or a chargino.  As discussed in section 2, in many cases the squarks do not directly decay to the LSP, especially in the model points selected for this study. If e.g. the LSP is Binolike, the lightest q L prefers to decay (if kinematically allowed) into the heavier Winolike neutralino or chargino. This is caused by a 3c w /s w ≈ 5.5 enhancement factor of the weak coupling with respect to the hypercharge coupling. If the squark happens to be compressed with the

JHEP05(2015)044
Winolike neutralino, jets from the squark decays are also soft and the remaining signature is determined by the branching ratios of the heavier Winolike neutralino. In these models the Winolike neutralino can have a large branching ratio to Higgs bosons. If it is a χ 0 4 , then decays to charginos are also possible. If the chargino decay is dominant, SUSY searches with leptons might be sensitive to these points. Searches asking for one lepton in the final state typically exclude simplified models with degenerate squarks decaying to charginos if m q < 800 GeV and m χ 0 1 < 300 GeV (see e.g. [64]). After applying all other search constraints we find no model that fulfills these requirements. If the Winolike neutralino is a χ 0 2 , then explicit searches for Higgs production from squarks might give a unique discovery possibility. Similar multi-step decays are possible in other cases as outlined in section 2.

Branching ratios
The branching ratios of all SUSY particles for decays into the lightest Higgs boson h 0 have been determined for all surviving models. These branching ratios include also multi-step decays to the Higgs boson, i.e. particles can have a non-zero branching ratio even if they do not couple to the Higgs boson directly. We show in figure 7 the branching ratio for all MSSM particles to the light Higgs boson h 0 . All models are shown in grey in order to indicate the ranges of these branching ratios. The possible decay processes have been described in more detail in section 2.
The sfermions can have decay branching ratios of up to 0.4, with the values for left-and right-handed sfermions strongly varying from model to model. The b 2 and t 2 squarks have a larger branching ratio due to the direct decay f 2 → f 1 + h 0 . As explained in section 2, the χ 0 2,3 neutralinos can have branching ratios close to unity. The χ ± 2 charginos can have a branching ratio that substantially exceeds 0.35 due to multi-step decays. The branching ratios of the heavy Higgs bosons range up to ∼ 0.4 due to direct as well as multi-step decays. Some models with interesting features are shown in colour. These models, which are labeled A-E, are shown in table 3 and are discussed in more detail below.
If we go one step further and take a look at all supersymmetric particles that decay into at least two Higgs bosons, the heaviest neutralino has the highest branching ratio, as can be seen in figure 8. Although it is not the preferred decay channel, it can decay into a χ 0 2,3 , which can subsequently decay into a LSP. Both decays can produce one h 0 boson. The sfermions can decay into two h 0 bosons via the intermediate decay into a heavy neutralino. The H 0 boson can directly decay to two h 0 bosons, the A 0 and H ± bosons decay via heavy neutralinos/charginos. The t 2 top squark can decay to t 1 + h 0 and t 1 can subsequently decay to one more h 0 .

Event generation, fast simulation and analysis
In order to determine the phenomenological relevance of h 0 production via SUSY processes, the LHC production rate needs to be determined. The generation of simulated events of pp collisions at a centre-of-mass energy of 14 TeV for each candidate model utilises JHEP05(2015)044

JHEP05(2015)044
PYTHIA6.4, HepMC 2.04.02 and DELPHES 3.0, just as described in section 3.4. In the simulation the branching ratio of the h 0 boson into two photons has been set to unity manually. This is done in order to prevent interference of jets or leptons originating from h 0 decays when analysing the jet multiplicities in the final states later on in this section. All other decay channels of the h 0 boson have been assigned a zero branching ratio. This does not affect the total h 0 production rates. A total number of 100000 events are generated for each of the approximately 250 candidate models that survive all previous constraints. All SUSY production processes are enabled.

Determining the expected number of events with Higgs bosons
The number of expected Higgs (h 0 ), di-Higgs and tri-Higgs events is calculated for an integrated luminosity of 10 fb −1 for each SUSY production process. For each model point the branching ratios to the h 0 as well as the SUSY production cross sections σ ISU B i for each subprocess ISU B i are considered. All cross sections are determined by PYTHIA. No attempt was made to include NLO corrections. In general these NLO corrections would further increase the production rate. So, in that sense our estimates are conservative. The number of events with at least one, two or three h 0 bosons is calculated as Figure 9 shows the rate of events with ≥ 1h 0 for all SUSY production processes normalized to an integrated luminosity of 10 fb −1 . The most important classes of production processes are squark-(anti)squark production, in particular for left-handed squarks (see figure 23 in the appendix), chargino-neutralino production (see also figure 19 in the appendix) and neutralino pair production (see also figure 20 in the appendix). Next in line are the associated production of a neutralino/chargino and a light squark (see also figures 21 and 22 in the appendix), and the production of pairs of bottom or top squarks (see also figure 24 in the appendix). Due to the nature of the mixing matrix, h 0 production via neutralino-pair and charginoneutralino processes are correlated. Large Higgs production rates are possible if the heavier neutralinos/chargino are relatively light and decay to h 0 . Examples are model D and E shown in table 3.

JHEP05(2015)044
Model point Largest h 0 production process Second largest h 0 production process Table 3. Parameters of the SUSY models with large Higgs production cross sections as discussed in the text.

JHEP05(2015)044
As can be seen in model B or C, Higgs production from squarks can still be large when at the same time Higgs production via chargino/neutralino processes is supressed. This happens when the charginos and non-LSP neutralinos are too heavy (about 600-800 GeV in model B and C) for direct production. These heavy neutralinos/charginos can, however, still be produced in the decay of a slightly heavier squark. As described in sections 2 and 3 such squark decays might be dominant. Searching for squark production might then be the only possibility to detect these models. Similar neutralino/chargino decays are important to produce light Higgs bosons in bottom squark decays.
The associated production of a chargino/neutralino and a squark can also be interesting for h 0 production. For the same mass of the produced particles the associated production cross section is in between the (electroweak) chargino/neutralino production and the (strong) squark production. This process can be important if the mass of the chargino/neutralino is similar to the mass of the squark, i.e. if one of the squarks is rather light. As explained before squarks can still be light in our models, e.g. if the squark decays via a heavy chargino/neutralino rather than directly to χ 0 1 . It is then difficult to detect the squarks in the conventional way at the LHC. An example is model C where the lefthanded squark (mass 760 GeV) decays with 65% branching ratio into the χ ± 1 and with 30% branching ratio into χ 0 2 (both with a mass of 627 GeV). The χ 0 2 decays with 85% branching ratio into a h 0 boson.
Higgs-boson production via top squarks can be enhanced for light stops. An example is model D, which has a t 1 mass of 850 GeV and a t 2 mass of 1130 GeV. Both stops decay to h 0 bosons predominantly via heavy neutralinos with branching ratios of 20 − 25%. This gives rise to special final-state topologies, involving top quarks, (possibly multiple) Higgs bosons and missing transverse momentum. Higgs-boson production via gluinos proceeds through the decay into light squarks. In the case of model D, these light squarks are top squarks, leading to spectacular event topologies where the gluino (or even both gluinos) can decay into h 0 tt χ 0 1 . The most important Higgs production processes are summarized in table 4. In some models (e.g. model D) h 0 production is significant for almost all important SUSY production processes. When all contributions to h 0 production from SUSY interactions are summed up, realistic models are found that lead to about 3000 events with at least one h 0 for 10 fb −1 of data. In almost all models a significant amount of missing tranverse momentum due to the LSP's is expected. This makes the events different from h 0 production via SM processes.

Di-Higgs production via SUSY processes
Since SUSY particles are pair produced and both particles can decay to a light Higgs boson, SUSY processes can be a significant source of events containing two h 0 bosons. We will show that this di-Higgs production rate can be significant.

JHEP05(2015)044
Di-Higgs production is of utmost importance in the SM to measure the triple-Higgs coupling. As such, the measurement of di-Higgs production is a central research objective for the high luminosity phase of the LHC. The SM di-Higgs production has an expected next-to-next-to-leading order (NNLO) cross section of roughly 40 fb [65], leading to about 400 events for 10 fb −1 . These events are very difficult to detect due to the enormous SM background rate.
In the MSSM another important source of di-Higgs events is the production of heavy Higgs bosons (see figure 25 in the appendix). Model A predicts an enourmous rate of > 2000 di-Higgs events, visible as a di-Higgs resonance. Heavy Higgs production is discussed separately in the next subsection. Di-Higgs events from processes involving SUSY particles are different due to the presence of large missing transverse momentum. The background from SM processes can be reduced by a large factor with cuts on this quantity. Figure 10 shows the di-Higgs production rate per SUSY process normalized to 10 fb −1 . Model C predicts the largest SUSY production rate for di-Higgs events with about 350 events for 10 fb −1 . This rate can also be compared with 10 and 4.2 events expected from the SM tth 0 h 0 or Zh 0 h 0 processes, which have a cross section of 1.0 fb at leading order [66,67] and 0.42 fb at NNLO [66], respectively. SUSY processes can therefore significantly enhance di-Higgs signatures in SM di-Higgs searches. Any deviation from the SM expectations in these searches needs therefore to be interpreted carefully, since deviations could be the result of SUSY decays.
The SUSY di-Higgs production is dominated by squark processes, followed by the direct production of heavy neutralinos/charginos. The most important SUSY di-Higgs production processes are summarized in table 4.

Tri-Higgs production via SUSY processes
Due to the multi-step decays of heavy neutralinos there is the possibility that one heavy neutralino can decay to two h 0 bosons. The corresponding branching fractions were discussed in section 4.1 and χ 0 4 was found to be the dominant source. This makes it possible to produce three Higgs bosons in one event. Figure 11 shows the number of tri-Higgs events per SUSY process normalized to 10fb −1 . Up to 20 tri-Higgs events can be produced, predominantly via squark production. The dominant tri-Higgs production processes are summarized in table 4. The SM tri-Higgs production cross section is 0.044 fb [68] leading to an expectation of only 0.4 events for 10 fb −1 . Tri-Higgs production might become important for large luminosities or, after a LHC discovery, for determining e.g. the neutralino mixing matrix.

The lightest Higgs boson from heavy-Higgs production processes
For the sake of completeness, simulated events with primary interaction processes involving heavy Higgs particles are also investigated briefly. This investigation utilises the calculation of events with at least one h 0 boson according to equation (4.1), but this time only the branching ratios of the heavy Higgses into one or more light Higgs boson(s) are taken into account.  As can be seen in figure 12, for most models the h 0 event rates from heavy Higgs production processes is low. This is caused by the decoupling limit. Due to the mass constraint on the lightest Higgs-boson, most models have an A 0 boson that is much heavier than the Z boson. In this decoupling limit all heavy Higgses are nearly mass degenerate and truly heavy. As a result, the heavy-Higgs production cross sections are relatively small and the h 0 event rates rather modest. The models with parameters that place them in the decoupling limit only reach a maximum of about 50 h 0 events for single heavy-Higgs production. Exceptions are a couple of models including model A, which have a smaller value for M A and which are therefore less firmly in the decoupling limit. These models also have a relatively small value for tan β, which results in a noticeable H 0 → 2h 0 branching ratio (see table 3) and substantially larger h 0 event rates beyond 1000 events. Also heavy Higgs production can have an effect on the Standard Model Di-Higgs production rate as discussed in [69].

JHEP05(2015)044
It must be kept in mind, though, that the heavy Higgs particles are not strictly speaking supersymmetric particles and are therefore not expected to lead to events with a large missing transverse momentum in the detector due to the LSP.

Kinematic distributions for Higgs events from SUSY
Boost of the h 0 boson. When a supersymmetric particle decays into a h 0 boson, the mass difference between mother and daughter (initial and final) state can lead to a boost. In hadronic pp collisions the main contribution to h 0 production by SM processes is expected to be from gluon-gluon fusion, and to a lesser extent from W W/ZZ fusion. A second JHEP05(2015)044 Figure 13. Qualitative distribution (in %) of the h 0 -boson boost in terms of β h 0 = v h 0 /c for the main h 0 hadroproduction processes in the SM (red, blue) and in pMSSM models via SUSY processes (grey). The distributions are each normalised to unity. Only those pMSSM models are presented that predict more than 100 h 0 bosons produced via SUSY processes for 10 fb −1 . relevant contribution is expected from associated tth 0 production, which is expected to lead to h 0 bosons that are more boosted in view of the larger (top-quark) mass scale in the process. Both processes are shown in figure 13 in order to compare the h 0 -boson boost (β h 0 ) distributions originating from SUSY and SM processes. Due to the larger mass scale of the SUSY processes the h 0 bosons are on average more boosted, even more than in tth 0 production. In extreme cases a heavy SUSY particle with mass > 1 TeV decays to a h 0 boson and a SUSY particle with a mass of O(100 GeV), leading to a very large boost. As an opposite extreme, we find one case where a squark with m ≈ 1.5 TeV decays to a χ 0 3 with m = 1.17 TeV, which subsequently decays to a h 0 and a χ 0 1 with m = 1.04 TeV. In such compressed scenarios the h 0 boost is even lower than expected from SM processes.
Missing transverse momentum. Figure 14 shows the missing transverse momentum distributions from SUSY processes for the selected pMSSM models. The generated events are normalized to an integrated luminosity of 10 fb −1 . All models have on average large missing transverse momentum up to several 100 GeV, permitting the introduction of selection cuts of 100-200 GeV in order to reduce backgrounds from SM processes. The production of heavy Higgs bosons is not considered. Model A has low missing transverse momentum since the h 0 boson originates from a heavy H 0 boson, which does not decay to the LSP. Final states with h 0 bosons. After the detector response is simulated with DELPHES, the final states are determined. Selection cuts are applied, requiring the leptons, i.e. electrons or muons, to have a transverse momentum of at least 20 GeV. For the jets this lower limit is chosen to be 50 GeV. B-jets and hadronic tau decays are counted as jets. Both leptons and jets are only considered if they are located within the pseudorapidity 4 region of |η e,µ,jet | < 2.5. In addition to the overlay removal that is automatically performed in DELPHES, a stricter overlap removal of ∆R a,b > 0.6 is applied. 5 The generated events are again normalised to an integrated luminosity of 10 fb −1 .

JHEP05(2015)044
The lepton and jet multiplicities for events requiring a missing transverse momentum of at least 100 GeV are shown in figures 15, 16 and 17. The most populated channels for single Higgs production are the channels that contain ≥ 1 − 4 jets, with close to 1000 events, and mono-higgs 6 production, with up to 200 events for 10 fb −1 . In some cases very high jet multiplicities can occur, as can be seen in figure 17. Channels with one lepton lead to ∼ 100 events and channels with two leptons to less than 10 events. Higher lepton multiplicities are not important for h 0 production. Di-Higgs and tri-Higgs production is dominantly found in channels with ≥ 2 jets. Another notable feature is that the production of neutralino pairs can lead to events with two Higgs bosons, missing transverse momentum and nothing else, i.e. no leptons and no jets. 4 η = − ln tan( θ 2 ) in terms of the polar angle θ w.r.t. the beam axis. 5 ∆R a,b = (∆η) 2 a,b + (∆φ) 2 a,b in terms of the pseudorapidity difference ∆η and the difference in azimuthal angle ∆φ between the objects (leptons/jets) a and b. 6 The 0-lepton, 0-jet channel, which is dominated by neutralino-pair production.

Conclusion
We have systematically investigated the possibilities to produce a 125 GeV Higgs boson (h 0 ) via SUSY processes within the phenomenological MSSM (pMSSM). We find the following interesting features: • Given global constraints on the pMSSM, it is possible to produce Higgs events with a large rate in the upcoming LHC data at the increased centre-of-mass energy. We have found valid pMSSM models that could produce more than 3000 Higgs, 300 di-Higgs and/or 20 tri-Higgs events already with an integrated luminosity of 10 fb −1 .
• A relation is observed between large Higgs-production rates via squark decays to heavy neutralinos and inherent difficulties to exclude such models in conventional (non-Higgs) LHC searches. This is caused by the fact that Higgs production requires a less compressed neutralino mass spectrum, which can bring the heavy neutralinos closer to the lowest-lying squark states, thereby reducing the available amount of energy for additional jets.
• In some models Higgs production is significant for almost all important SUSY production processes, which can have large repercussions on SM Higgs studies and SUSY searches.
• Higgs production via SUSY processes might significantly enhance the event rates for SM Higgs and di-Higgs searches, especially in final states with missing transverse momentum. The allowed SUSY production rates can be reduced by upcoming (negative) SUSY searches at higher LHC energies, especially if new dedicated searches for events with h 0 bosons and missing transverse momentum are performed.
• Higgs production processes can likewise be of importance for a SUSY discovery via "Higgs tagging". We found that the different SUSY production channels can lead to a large variety of interesting event topologies and kinematics. Of special interest are multi-jet channels with up to three Higgs bosons, "mono-Higgs" channels with up to two Higgs bosons, one-lepton channels with one Higgs boson and Higgs production in association with top or b-quarks, all with a sizeable amount of missing transverse momentum. The list is completed by searches for heavy Higgs bosons decaying directly or via neutralinos/charginos into h 0 bosons.

JHEP05(2015)044
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.