Higgs Production in Neutralino Decays in the MSSM - The LHC and a Future e+e- Collider

The search for the production of weakly-interacting SUSY particles at the LHC is crucial for testing supersymmetry in relation to dark matter. Decays of neutralinos into Higgs bosons occur over some significant part of the SUSY parameter space and represent the most important source of $h$ boson production in SUSY decay chains in the MSSM. We study h production in neutralino decays using scans of the phenomenological MSSM. Whilst in constrained MSSM scenarios the decay chi^0_2 ->h chi^0_1 is the dominant channel, this does not hold in more general MSSM scenarios. On the other hand, the chi^0_2,3 ->h chi^0_1 decays remain important and are highly complementary to multi-lepton final states in the LHC searches. The perspectives for the LHC analyses at 8 and 14 TeV as well as the reach of an e+e- collider at 0.5, 1, 1.5 and 3 TeV are discussed.


Introduction
The discovery of a Higgs-like boson by the ATLAS and CMS experiments [1,2] at the CERN LHC opens up an era of precision studies of its production and decay properties. In particular, establishing if the newly discovered particle is the Standard Model (SM) Higgs boson or the manifestation of an extended Higgs sector is a key question. If the observed Higgs-like particle is the lightest Higgs boson, h, of a Supersymmetric extension of the SM (SUSY), it could be significantly produced also in the decay chains of supersymmetric particles, in particular those of neutralinos,χ 0 j →χ 0 i h. There have been several phenomenological studies of h production in neutralino decays in various constrained supersymmetric models [3][4][5][6][7][8][9][10][11], including some detailed assessments of the LHC potential [7,12]. This paper discusses the regions of the MSSM parameters where these decays are relevant and the perspectives for their observation and study at the LHC and a future lepton collider in the framework of the phenomenological minimal supersymmetric extension of the SM (pMSSM) with 19 free parameters. This provides us with sufficient freedom to the masses and couplings to explore the supersymmetric parameter space in a largely unbiased way. The study of the pMSSM parameter space by high statistics, flat scans of its parameters is briefly described in section 2. Results from measurements and searches at LEP and LEP2, flavour physics, dark matter experiments provide already powerful constraints. These are discussed in section 3 together with the implications of the LHC results in the search for SUSY particles and in the first determi-nation of the properties of the newly discovered boson. Section 4 presents the regions of the viable pMSSM parameter space probed by the contribution of the h boson in neutralino decays, while section 5 discusses the expected event rates, experimental signatures and relevant kinematics of these processes at the LHC and a lepton collider, for various centre of mass energies. Section 6 has the conclusions.

MSSM Scans and Tools
The study of two-body neutralino decays into a h boson is based on the analysis of a large sample of MSSM points obtained through a flat scan of the pMSSM parameters, which are varied in an uncorrelated way within the following ranges: 1 ≤ tan β ≤ 60 , 50 GeV ≤ M A ≤ 2 TeV , −10 TeV ≤ A f ≤ 10 TeV , 50 GeV ≤ Mf L , Mf R , M 3 ≤ 3.5 TeV , 50 GeV ≤ M 1 , M 2 , |µ| ≤ 3.0 TeV (1) to generate a total of 1 × 10 8 pMSSM points. We also perform a scan of the CMSSM, in order to contrast the results for the pMSSM with those of a highly constrained model. For this we vary the CMSSM parameters within the following ranges: In all scans, we select the generated points which fulfil a set of constraints derived from flavour physics and lower energy searches at LEP2 and the Tevatron, as summarised below and discussed in more details in Ref. [13], to which we refer also for further details on the pMSSM scans. These selected points are referred to as "accepted points" in the following. The pMSSM scans employ a number of programs and software tools. Only those most relevant to this study, are mentioned here, while further details are given in [13]. SUSY mass spectra are generated SOFTSUSY 3.2.3 [14]. Decay branching fractions are calculated using SDecay [15,16]. Higgs decay branching fractions are calculated with HDECAY (5.0) [17]. The flavour observables and dark matter relic density are calculated with SuperIso Relic v3.2 [18]. Cross sections are computed with Pythia 6.424 [19], also used for event generation, while the NLO k-factors are evaluated using Prospino v2 [20]. The physics object response for the LHC analyses is simulated using the Delphes fast simulation package [21], tuned on the ATLAS detector performance and validated for the ATLAS multi-lepton neutralino and chargino analyses, as discussed below.

Constraints
We apply constraints from flavour physics, dark matter and SUSY searches at LEP, Tevatron and LHC. These have been discussed in details in Ref. [13]. In particular, we consider the decay B s → µ + µ − , which can receive significant SUSY contributions at large values of tan β. Here, we adopt the latest result from LHCb with the measurement of a branching fraction of BR(B 0 s → µ + µ − ) = (3.2 +1.5 −1.2 ) × 10 −9 with 2.2 fb −1 of data [22], from which we derive the constraints at 90% C.L. after accounting for theoretical uncertainties. Dark matter also provides significant bounds for this analysis. We impose the constraint on the neutralino relic density of 10 −4 < Ω χ h 2 < 0.155, derived from the WMAP-7 result [23], accounting for theoretical and cosmological uncertainties and allowing the neutralino be responsible for only part of the observed dark matter. In addition, we impose the new limit on dark matter direct detection for spin-independentχ-nucleon scattering recently obtained from 225 live days of the XENON-100 data [24].
Searches for SUSY particles at the LHC in channels with missing E T [25,26] have already provided a number of constraints relevant to this study, by excluding part of the pMSSM phase space, corresponding to gluino masses below ∼ 600 GeV to 1 TeV and scalar quarks of the first two generations with masses below ∼ 400 to 800 GeV. These constraints are included here using the same analysis discussed in Ref. [13], extended to an integrated luminosity of 12 fb −1 . In addition, we also include the results of the searches for events with large missing transverse momentum and two b-tagged jets, conducted by ATLAS [27]. Finally, we apply the constraints in the Higgs sector by selecting only points having a mass of the lightest Higgs boson in the range 123 < M h < 129 GeV and compatible with the constraints on the heavier MSSM Higgs bosons in the channel H/A → τ + τ − [28,29].

Neutralino decays to h: MSSM parameter regions
In the MSSM the couplings of the Higgs bosons to charginos and neutralinos are maximal for higgsino-gaugino mixed states, while the gauge boson couplings to neutralinos are maximal for higgsino-like states. The h couplings to neutralino are suppressed by powers of |µ|/M 2 or M 2 /|µ| in the gaugino-like and higgsino-like regions, respectively. The branching fraction for a set of our accepted pMSSM points is shown as a function of |µ|/M 1 and |µ|/M 2 in Figure 1. In the case where the charginos and neutralinos are gaugino-like (i.e. when the higgsino mass parameter is much larger than the gaugino mass parameter, |µ| M 2 ) or higgsino-like (i.e. in the opposite situation |µ| M 2 ), this results into the dominance of the decays of the heavier charginos and neutralinos into the lighter states and Higgs bosons, over the same decays with gauge boson final states [6,7]. This is also the case for the so-called "little cascade" in the gaugino region M 2 > ∼ M 1 |µ| where the branching fraction for the decay of theχ 0 2 into the LSP neutralino and the lighter h bosonχ 0 2 → hχ 0 1 , is in general larger than that for the decayχ 0 2 → Zχ 0 1 , when kinematically accessible in the two-body channel 1 . Thus, the decaỹ χ 0 2 → hχ 0 1 is important in these regions, unless sleptons are light and the additional channels whereχ 0 2 decays into a lepton and its super-partner are open. These decays are most important if the neutralino is gaugino-like, since the coupling to the higgsino component is suppressed due to the small lepton mass. It is important to observe that for small values of ||µ| − M 2 |, where the contribution of the decay into h is enhanced, the mass difference M χ 0 3 − M χ 0 2 is small, typically of the order of 50 GeV or less. In this regime, the production cross sections forχ 0 2χ ± 1 andχ 0 3χ ± 1 in pp collisions are comparable. Due to the nature of theχ 0 2 andχ 0 3 states, the branching fractions into h and Z of these two neutralinos are complementary, i.e. the yield into h for thẽ χ 0 2 is approximately equal to that into Z of theχ 0 3 states, and viceversa, as shown in Figure 2. This highlights the complementarity between the search forχ 0 2 andχ 0 3 with decays into h and Z and the interest in pursuing both of them in the LHC analyses.
We study the regions explorable by theχ 0 2,3 → hχ 0 1 decay using our set of accepted pMSSM points. These depend only on a subset of the pMSSM parameters, mostly . In order to characterise the occurrence of the various decays in different regions of the parameter space, we first study the fraction of the accepted pMSSM points into each bin in the [µ, M 1 ] and [µ, M 2 ] planes where hχ 0 1 , Zχ 0 1 and the sum of˜ andτ τ is the dominant two-body neutralino decay mode. We define the dominant decay mode as that having a branching fraction larger than both 0.2 and any of the other two-body decay channels. Results are shown forχ 0 2 in Figures 3 and  4, where we also show the region where the sum of branching fractions for the three-body decays exceeds both 0.2 and any of the individual two-body processes. While the values for the fraction of points obtained depend on the parameter ranges adopted for our scans given in Eq.(1), the picture emerging from Figure 4 is of general validity.
As we can see from Figure 3, the hχ 0 2χ 0 1 couplings are larger than the Zχ 0 2χ 0 1 whenever M 1 > ∼ |µ| and thẽ χ 0 1 is mostly bino. The region where |µ| ≈ M 1 is excluded as the decayχ 0 2 → hχ 0 1 is not allowed kinematically. Very small values of |µ| are excluded from the searches for charginos at LEP2, which imply that |µ| > ∼ 90 GeV. Theχ 0 2 → Zχ 0 1 mode is dominant for low values of µ where the LSP is higgsino-like and, thus, there is a strong coupling of the Z boson toχ 2 1χ 0 1 states. There is a line 1 is kinematically close and the decayχ 0 2 → Zχ 0 1 becomes dominant, provided that the sleptons are not light enough for the˜ channel to contribute. Finally, at low µ values,χ 0 1 andχ 0 2 are higgsino-like, but the splitting is typically small, so the two-body decays through h or Z are forbidden and the three-body modes dominate.

Neutralino decays to h: pMSSM to CMSSM comparison
W h production fromχ ± 1χ 0 2 decays at the LHC has been recently discussed in the context of the CMSSM in Ref. [30], where the authors find that the decayχ 0 2 → hχ 0 1 is dominant for typical choices of the model parameters. In our analysis we find indeed thatχ 0 2 → hχ 0 1 is the dominant decay, with branching fraction in excess of 0.20, for 71% of the accepted CMSSM points before applying the relic density constraint. This is no longer the case when imposing the full set of constraints. Of the accepted points withχ 0 2 → hχ 0 1 as the dominant decay, only 10 −4 fulfil the constraint 10 −4 < Ω χ h 2 < 0.155 and 8×10 −5 the tighter requirement of 0.068 < Ω χ h 2 < 0.155 set by the WMAP bound. In fact, the region where theχ 0 2 → hχ 0 1 is dominant corresponds to values of neutralino relic density, Ω χ h 2 , which largely exceed the upper limit of the WMAP result [23], as shown in Figure 5. This is due to the peculiar nature of the CMSSM, which has large µ values, with |µ| M 2 making the neutralinos to be gaugino-like, as observed above. The general prevalence of theχ 0 2 → hχ 0 1 decay, observed in the CMSSM before the Ω χ h 2 constraint, cannot be generalised to the MSSM. But in turn the greater flexibility of the pMSSM allows to reconcile the relic density constraint to large coupling of the h to the neutralinos (see Figure 5). For comparison, in the pMSSM the decay intoχ 0 1 h is dominant in only 10% of the accepted points before the relic density constraint. But now 28% of these fulfil also the loose relic density constrain of 10 −4 < Ω χ h 2 < 0.155 and still 2.2% make also the tighter requirement of 0.068 < Ω χ h 2 < 0.155. It is therefore important to reconsider the apparent dominance of theχ 0 2 →χ 0 1 h decay in the CMSSM, when examining the SUSY phenomenology in the context of more general models, or even in the CMSSM itself after imposing the dark matter relic density constraint. Still, the process of h production in neutralino decays is of crucial importance in specific regions of the parameter space, which are largely complementary to those yielding decays into sleptons or Z bosons, probed by the multi-lepton analyses at the LHC.

Neutralino decays to h: Experimental Signatures
In the following we discuss the perspectives for detecting and studying these decays at the LHC and an e + e − linear collider (LC). We consider for the LHC an integrated luminosity of 25 fb −1 at 8 TeV and 300 fb −1 at 14 TeV. For the LC we consider 500 fb −1 of integrated luminosity at √ s = 0.5 TeV and 1 ab −1 at 1, 1.5 and 3 TeV. For each of these scenarios, we combine the decay branching fractions with the relevant neutralino production cross sections to study the number of expected events. For the case of the LHC, we also characterise the kinematics of the detectable final states and study a possible event selection to isolate neutralino decays into hχ 0 1 (and Zχ 0 1 ). We compare these results to multi-lepton final states.

LHC
As a first step, we combine the results on the branching fractions discussed in the previous section to the expected relevant cross sections for the associated chargino neutralino production. Through the decayχ ± 1 → W ±χ0 1 , W ± → ± ν, this process produces a high p T lepton, which is very valuable for the event trigger and subsequent selection. The cross sections for chargino neutralino production at 8 and 14 TeV, followed byχ 0 2,3 → hχ 0 1 ; h → bb, τ τ and χ ± 1 → W ±χ0 1 , W ± → ± ν are given in Figures 6 and 7, respectively.
Then, we study the fraction of accepted pMSSM points in bins of the M 1 , M 2 and µ parameters, which yield more than 10 signal events forχ 0 2 → hχ 0 1 andχ 0 3 → hχ 0 1 with 25 fb −1 at 8 TeV and 300 fb −1 at 14 TeV. Results are summarised in Figures 8 and 9, respectively. Again the values for the fractions of points given in these figures depend on the the ranges adopted in the scan at Eq.(1). However, the location of the regions highlighted by the large occurence of points with enough signal events are of general validity. Due to the factor of ∼3 increase of the pp →χ ± 1χ 0 2,3 production cross section and the factor of ∼10 increase in the assumed statistics, the coverage of the parameter space expands significantly when moving from 25 fb −1 at 8 TeV to 300 fb −1 at 14 TeV. Still, the 8 TeV data provide interesting limits, when comparing to the other channels sensitive to SUSY weak production at the LHC. Results from searches for neutralino and chargino production have been reported for the two-and tri-lepton channels by both ATLAS and CMS [31][32][33]. These channels are sensitive to the contribution of either light sleptons in the decays χ 0 2 →˜ ;˜ → χ 0 1 , or Z bosons in the decayχ 0 2 → Zχ 0 1 ; Z → . In the unconstrained MSSM, slepton masses can be pushed well above theχ 0 2 mass, so that the only remaining dominant two-bodyχ 0 2 decays are either hχ 0 1 or Zχ 0 1 . This motivates us to pursue the study of the challenging hχ 0 1 channel, through the bb + E missing T and τ + τ − + E missing T topologies, which complement the Zχ 0 1 channel. In addition, the decaysχ 0 2 →τ τ are of special importance in the region of lightτ 1 , which is highlighted by a possible enhancement of h → γγ yields in the ATLAS and CMS data. This process leads to the same τ + τ − + E missing T final state as the decays hχ 0 1 ; h → τ τ and Zχ 0 1 ; Z → τ τ . If a signal is observed in the bb + E missing T and τ + τ − + E missing T modes it would become essential to attempt to identify the contributing channels, possibly through a fit to the bb invariant mass and the τ τ transverse mass. In order to compare the coverage of the space of the relevant pMSSM parameters for the Higgs and the multilepton channels, we compare the results obtained above to the fraction of the accepted pMSSM points for which a signal can be obtained in the two-and three-lepton channels across the [µ − M 1 ] and [µ − M 2 ] planes at 8 and 14 TeV, based on the signal selection of the ATLAS analyses. Accepted pMSSM scan points fulfilling the selection criteria discussed above and having an inclusive chargino and neutralino production cross section larger than 0.01 fb at 8 TeV collision energy are studied in details. A set of 20k SUSY events with inclusive neutralino and chargino production is generated for each point at both 8 and 14 TeV. Simulated events are passed through a fast parametric simulation using the Delphes package which simulates the physics objects used in the subsequent analysis. The selection criteria of the two-and three-lepton + missing E T ATLAS analyses [31,32] are applied. The observability of the pMSSM point is assessed based on the corresponding Cl value, where the number of expected signal events is obtained from the result of the selection on the generated events rescaled to a luminosity of 25 fb −1 and that  of background events is taken from that estimated for the ATLAS analyses, rescaled again to the assumed luminosity. The results of this simulation is validated by comparing the 95 % C.L. exclusion contours obtained in the M 2 -µ plane with those expected for the ATLAS analyses, under the same assumptions. These contours agree within  decay, based on the Delphes fast simulation response. In this analysis, jets are reconstructed using the anti-kt algorithm [34] with a cut 0.4. Leptons are selected requiring p T > 10 GeV and pseudo-rapidity below 2.40 for electrons and 2.47 for muons and to be isolated. We select events with only one lepton and two jets with E T > 25 or 45 GeV. We further require the reconstructed jet to be associated to a generated b hadron. These selection criteria are inspired to those adopted in the W H Higgs analysis with H → bb, leading to a similar final state, but with smaller values of E missing T . Here, we do not investigate the backgrounds, for which the tt and W + jets processes are expected to give the dominant contributions. In our signal analysis, we require significant missing transverse energy, E missing T in the event, due to the two escapingχ 0 1 s. The average value of the transverse missing energy and of the average energy of the reconstructed bb system iñ χ 0 2,3 → hχ 0 1 , h → bb for signal events from the pMSSM points which pass the cuts above, is shown in Figure 11. We perform the analysis for values of the cuts on the b jets p T of 25 and 45 GeV and on E missing T of 50, 100 and 150 GeV. The efficiencies at 14 TeV for b jet p T 25 GeV are typically in the range 15 to 30% depending on the requirement on E missing T . Figure 12 shows the fraction of the pMSSM points in the [µ, M 1 ] parameter space where the analysis selects more than 20 signal events. In the case thẽ χ 0 2 andχ 0 3 states are close in mass, they would be both produced and they may decay into hχ 0 1 and Zχ 0 1 , respectively, as discussed above. In this case the bb + E missing T final state will receive contribution from both h → bb and Z → bb, which can be resolved at least on a statistical basis from an analysis of the invariant mass of the bb system.

e + e − Linear Collider
An e + e − collider of sufficient energy and luminosity is very well suited for reconstructing two-body decays of pair-produced neutralinos, e + e − →χ 0 iχ 0 j , including those decaying into h and Z bosons. In particular, since the energy of the producedχ 0 i is known, apart from beam-related effects, the energy distribution of the reconstructed boson can be used to precisely determine theχ 0 i mass. The decayχ 0 2 → Zχ 0 1 was studied as a benchmark reactions at √ s = 0.5 TeV using full simulation and reconstruction as part of the ILC LoIs [35,36]. In these analyses theχ 0 2 mass of 217 GeV was determined with a statistical relative accuracy of 1%. The decaysχ 0 2 → hχ 0 1 andχ 0 2 → Zχ 0 1 have been studied for CLIC at √ s = 3.0 TeV [37]. Thẽ χ 0 2 → hχ 0 1 decay, followed by h → bb, leads to the signature bbbb+ missing energy final state. Theχ 0 2 decay can be isolated in an inclusive SUSY samples, where its mass and h yield are determined with a statistical accuracy of 4% and 5%, respectively. The pair production cross sections multiplied by thẽ χ 0 2,3 → hχ 0 1 branching fractions are shown in Figure 13 as a function of the neutralino mass, for the four √ s energy values (0.5, 1.0, 1.5 and 3.0 TeV) chosen for this study.
The production of neutralino pairs in e + e − collisions is mediated by s-channel Z boson exchanges as well as t and u-channel left-and right-handedẽ exchange. If the sleptons are assumed to be very heavy or the produced neutralinos are higgsino-like (and thus the couplingχ 0ẽ e is negligibly small) only the s-channel Z boson exchange is relevant. In fact, the neutralino cross sections are rather small, even in the case ofχ 0 1χ 0 2 production where the phase space is more favourable. The only exception is when both χ 0 1 andχ 0 2 are higgsino-like and have maximal couplings to the Z boson. Also the cross sections with identical neutralinos are extremely small. Due to Fermi statistics, the neutralinos are produced in p-waves and therefore suppressed at threshold and, in the case of gauginos or hig-  4 are significant, except in the case where one of the neutralinos is a pure gaugino, which leads to suppressed Zχχ couplings. It is important to observe that, in the region of parameters where theχ 0 2 → hχ 0 1 process is the dominant neutralino decay, theχ 0 2 tends to decouple from the Z. This suppresses the e + e − →χ 0 2χ 0 2 pair production, leaving instead the e + e − →χ 0 2χ 0 3 process as the only significant neutralino production mechanism, when kinematically available. This correlation between the pair production cross section and the decay branching fractions is demonstrated in Figure 14 for a sample of accepted pMSSM points where the sum of the neutralino masses is more than 100 GeV below the assumed collision energy, to avoid threshold effects. With an expected integrated lu- minosity of 0.5 -1.0 ab −1 a high energy linear collider can observe these decays over a significant part of the kinematically allowed region of the parameter space. This is quantified here by the fraction of the pMSSM points yielding at least 20 signal events for an integrated luminosity of 0.5 ab −1 for √ s = 0.5 TeV and of 1 ab −1 for √ s = 1, 1.5 and 3 TeV, over the [µ, M 1 ] and [µ, M 2 ] parameter space, as shown in Figure 15.
Finally, we study the determination of the yield of h bosons in an inclusive SUSY sample for the specific case of a MSSM model havingχ 0 1 with mass of 340 GeV, χ 0 2 with mass of 643 GeV decaying predominantly into hχ 0 1 andχ ± 1 with mass of 643 GeV decaying exclusively into W ±χ0 1 and M h = 126.1 GeV, in 3 TeV e + e − collisions. An inclusive sample of SUSY events, corresponding to 0.5 fb −1 of integrated luminosity, is generated using Pythia and events are fully simulated using Geant-4 [38] and reconstructed following the same analysis procedure discussed in [37]. Events are pre-selected requiring a visible energy 0.08 √ s < E tot < 0.6 √ s, an energy in charged particles larger than 150 GeV, transverse energy larger than 200 GeV and jet multiplicity 2≤ N jets <5. Jet clustering is performed using the Durham jet algorithm [39], with y cut = 0.0025, on the reconstructed particle flow objects of the Pandora particle flow package [40]. Selected events are clustered into four jets and the di-jet invariant mass for all the three possible pairings is computed.
The jet pairing minimising the difference between the dijet invariant masses is selected, provided the mass difference is below 20 GeV. The fraction of background nonresonant events is obtained from the di-jet mass sidebands 20 < E jj < 60 GeV and 140 < E jj < 200 GeV. The fraction of W + W − , Z 0 Z 0 and h 0 h 0 events in the inclusive 4-jet SUSY events is extracted by a χ 2 fit to this di-jet mass distribution. The W ± and Z 0 mass peaks are parametrised as Breit-Wigner functions convoluted with a Gaussian term describing the experimental resolution. The mass and width values of the Breit-Wigner functions are fixed to their generated values, while the total area and the width of the Gaussian resolution terms are left free in the fit. The h 0 peak, which has negligible natural width, is modelled as the sum of two Gaussian curves, one representing the correctly reconstructed signal events, centred at the nominal M h value, the second describing decays where the mass has a lower reconstructed value due to semi-leptonic b decays. The central value, width and fraction of events in this second Gaussian is extracted by a fit to a pure sample of decays into h 0 bosons and fixed in the fit, while the Gaussian width of the main peak is kept free. The yield of h bosons, with mass 126.1 GeV, is extracted by a χ 2 fit to the di-jet mass distribution of events in 4jet final states, shown in Figure 16 for all di-jets and for those passing di-jet b-tagging based on the ZVTOP algorithm [41] implemented in the LCFIVertex package [42]. We measure the fraction of h bosons to be 0.269±0.013 and that of Z bosons to be 0.037±0.016, after the nonresonant background subtraction, which compares well to the original values of 0.290 and 0.025, respectively, of the generated events. This specific example shows the feasibility to accurately determine the yield of h bosons in SUSY particle decays in the data of an e + e − collider operating at sufficient energy.

Conclusions
The observation of a Higgs-like particle with mass 126 GeV promotes its study in the decay of supersymmetric particles. Neutralino decays offer a prime opportunity, since the decayχ 0 2,3 → hχ 0 1 is highly complementary to the channels leading to multi-lepton final states, such as Z and˜ . The typical products of cross sections and decay branching fractions to the bb final state is of the order 0.01 -1 fb in 8 TeV pp collisions and a factor of ∼ three larger at 14 TeV for M χ 0 2,3 ≤ 600 GeV. This should make possible to search for this decay in the bb + MET over a significant region of the parameter space at 14 TeV, complementary to that covered by searches with multi-leptons, provided the tt and W + jets backgrounds can be properly suppressed. An e + e − linear collider of sufficient energy, √ s ≥ 1 TeV, and luminosity, yielding ≥0.5 ab −1 of data, can study h production in neutralino decays for masses up to the kinematical limit for pair production, obtaining O(1%) precision on the mass and O(5%) on the yield of Higgs bosons from neutralino decays in inclusive multi-jet SUSY events.