Di-Higgs signatures from R-parity violating supersymmetry as the origin of neutrino mass

Motivated by the naturalness and neutrino mass generation, we study a bilinear R-parity violating supersymmetric scenario with a light Higgsino-like lightest super-symmetric particle (LSP). We observe that the LSP can have substantial decay branching ratio to νh in a large part of the parameter space, and thus study the pair production of electroweakinos followed by the decays χ˜1±→χ˜10W±∗\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\tilde{\chi}}_1^{\pm}\to {\tilde{\chi}}_1^0{W}^{\pm \left(\ast \right)} $$\end{document} and χ˜10→νh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\tilde{\chi}}_1^0\to \nu h $$\end{document}. This leads to an interesting signature of Higgs boson pair production associated with significantly large missing transverse energy which is grossly distinct from the di-Higgs production in the Standard Model. We investigate the perspective of probing such signatures by performing a detector level simulation using a toy calorimeter of both the signal and corresponding backgrounds for the high-luminosity high energy phase of the Large Hadron Collider (LHC). We also advocate some observables based on kinematical features to provide an excellent handle to suppress the backgrounds.


Introduction
Supersymmetry (SUSY) has been considered as one of the theoretically well motivated framework to protect the electroweak scale from the quadratic divergence caused by a certain ultra-violet (UV) physics. As a bonus, supersymmetry provides a viable candidate for dark matter of the Universe: the lightest supersymmetric particle (LSP) which is neutral and stablized by assuming R-parity conservation. 1 On the other hand, R-parity violation (RpV) brings another interesting possibility of generating tiny neutrino masses  in the context of the Minimal Supersymmegtric Standard Model (MSSM). The observed neutrino masses and mixings determine the lepton flavor structure of R-parity violating couplings which typically leads to clean signature of same-sign dileptons and predicts specific leptonic branching ratios of the LSP decayχ 0 1 → l ± W ∓ [28][29][30][31][32][33][34][35][36][37]. As no hint for supersymmetry appeared yet at the Large Hadron Collider (LHC), the naturalness argument for the TeV scale SUSY is in question due to a severe fine-tuning which turns out to be much more than expected. The electroweak symmetry breaking in SUSY requires a potential minimization condition: where m H u,d are the soft masses of the two Higgs doublets, tan β ≡ v u /v d is the ratio of their vacuum expectation values, and µ is the Higgs bilinear parameter in the superpotential. As the LHC pushes up the soft mass scale above TeV range, the condition (1.1) requires
In this paper, we investigate the LHC signatures of the light Higgsino in association with bilinear R-parity violation (BRpV) as the origin of the observed neutrino masses and mixings. Contrary to the conventional studies on BRpV predicting a peculiar signature of same-sign dileptons from pp →χ 0 1χ 0 1 → l ± l ± W ∓ W ∓ , we will focus on the unusual case of the LSP decay dominated by the Higgs channelχ 0 1 → νh which will be shown to occur in a large region parameter space of the scenario under consideration. As a consequence, it leads to an interesting LHC signature of di-Higgs bosons with missing transverse energy. Measurement of Higgs-pair production cross-section will be one of the main focuses of the high energy and high luminosity LHC run. It is also an important step towards our understanding of the electroweak symmetry breaking mechanism. At LHC energies, Higgs boson pair production occurs dominantly through the gluon fusion in the SM [50,51]. Other processes, such as weak boson fusion qq ( ) → qq ( ) hh, associated productions qq ( ) → W hh, Zhh or associated production with top quarks, gg, qq → tthh also occur albeit with cross-sections which are 10-30 times smaller than the gluon fusion [52,77]. Di-Higgs production at the LHC has been studied in the context of triliniear Higgs self coupling measurements by various authors [53][54][55][56][57][58][59][60] and references there in. While the SM crosssection is small and a lot of efforts have been put in to observe it at the LHC, it also plays a crucial role in the context of new physics searches at the LHC since physics beyond the SM can lead to an enhancement of the observable cross-sections and/or different event kinematics. We demonstrate that di-Higgs in association with a large / E T can be very useful to probe the BRpV SUSY where conventional channels fail to be sensitive. Assuming only the electroweak production of the electroweakino pairs, we analyze the di-Higgs signal in the channel of γγbb/ E T at the LHC14 with the integrated luminosity ranging from 1-3 ab −1 . Search for Higgs-pair production as a window to probe new physics is one of the major activities in the context of LHC and future collider, e.g., resonant Higgs-pair production in the context of singlet and doublet extension of the SM [61][62][63][64], double Higgs production via gluon fusion in the effective field theory framework [65][66][67], Higgs pair production in the context of SUSY extension of the SM [68][69][70][71][72] and various other extensions of it [73][74][75][76].
This paper is organized as follows: in section 2 summarizing the results of refs. [31-33, 36, 37], we provide a brief review on the bilinear RpV couplings constrained by the resulting tree-level neutrino mass matrix. In section 3, we analyze the LSP decay modes for various choices of the Higgsino (µ) and gaugino mass parameters (M 1,2 ) In section 4, we compute the Higgsino pair production cross-section for some benchmark points and perform a detector-level simulation using a toy calorimeter and PYTHIA hadronisation and showering algorithm for both signal and backgrounds in the di-Higgs decay channel of hh → γγbb and obtain the LHC14 perspective to probe our scenario. Finally, we summarize our results and conclude in section 5. In appendix A, we collect the effective R-parity violating couplings relevant for the Higgsino decays and appendix B shows the decay widths of the neutralinos induced by the BRpV couplings.

JHEP12(2016)062 2 Bilinear RpV and neutrino mass matrix
Allowing lepton number violation in the supersymmetric standard model, the superpotential is composed of the R-parity conserving W 0 and violating W 1 part; Among soft supersymmetry breaking terms, let us write R-parity violating bilinear terms; It is clear that the electroweak symmetry breaking gives rise to nonzero vacuum expectation values of sneutrino fields,ν i , as follows: Given the BRpV couplings i and a i , the neutrino-neutralino sector form a 7 × 7 mass matrix whose 3 × 4 (Dirac) neutrino-neutralino mass matrix takes the form of where s W ≡ sin θ W is the weak mixing angle, and the index i runs for three neutrino flavors (ν e , ν µ , ν τ ), and j runs for the neutralino states (B,W 3 ,H 0 1 ,H 0 2 ) which has the usual 4 × 4 mass matrix M N containing the bino, wino and Higgsino masses denoted by M 1 , M 2 and µ, respectively.
As is well-known, a seesaw diagonalization rotating away M D [see appendix A for details] generates the "tree-level" neutrino mass matrix M ν = −M D M N −1 M D T whose components are given by This makes massive only one neutrino, ν 3 , in the direction of ξ. The other two get masses from finite one-loop corrections and thus ν 3 is usually the heaviest component. We fix the value of m ν 3 from the atmospheric neutrino data and thus the overall size of ξ ≡ | ξ| is determined to be (2.6) Furthermore, among three neutrino mixing angles defined by the mixing matrix with c ij = cos θ ij and s ij = sin θ ij , etc., two angles are determined by the tree-level mass matrix (2.5) as follows: These two angles define the atmospheric and reactor neutrino oscillation angles, respectively, and thus one has sin 2 2θ 23 ≈ 1 and sin 2 2θ 13 ≈ 0.09 [78]. This implies that the sizes of ξ i should follow the relation: The other angle θ 12 can be determined only after including one-loop corrections which are assumed to be smaller than the tree-level contribution (2.5) and thus irrelevant for our discussion.
The BRpV terms induce mixing between neutrinos (charged leptons) and neutralino (charginos) as well as their scalar partners. Rotating them away, one gets the effective RpV vertices of neutralinos and charginos which are summarized in appendix A.

Light Higgsino decays
A distinct feature of the RpV SUSY models is that the LSP,χ 0 1 , is not stable. In our analysis, neutralino decays via sfermions are highly suppressed in the limit of heavy sfermion masses and small trilinear RpV couplings responsible for one-loop neutrino mass generation. Therefore, we discuss the branching ratios (BR) of the LSP for the following decay processes:χ where is any of the three charged-leptons.
To explore the phenomenological features, we consider specific scenarios in this work where |µ| < M 2 < M 1 and |µ| < M 1 < M 2 leading to a situation of the Higgsino LSP, with the lightest statesχ 0 1,2 andχ ± 1 being Higgsino-like. In the following, we fix M 1 = 500 GeV, 700 GeV and 1 TeV, and vary both µ and M 2 . Given M 1 and M 2 , the masses of electroweakinos are determined by the value of µ. We also vary tan β to see the variation of branching ratios with respect to µ, M 2 and tan β ( figure 1 and 2). Both the figures 1 and 2 show that the branching ratio for the W or νh decay mode is quite sensitive to the values of µ, M 2 and tan β.
The main contribution of LSP decay into νh mode comes from the neutrino-Higgsino mixing introduced in the R-parity conserving vertex of type: Higgs-Higgsino-gaugino interaction ((φ * T a ψ)λ a ) and also from sneutrino-Higgs mixing introduced in the sneutrinonetrino-guagino vertex. All of the RpV vertices which we have obtained in eqs. ξ 1 : ξ 2 : ξ 3 = 0.1 : 1 : 1 (eq. (2.9)). Thus, BR(χ 0 1 → eW ) is quite suppressed and major contribution to lW mode comes from µW and τ W channels. An important point to note is that the decay rate (χ 0 1 → W ) also depends on charged-lepton masses through c R 2 (eq. (B.3)). Although this mass dependence is suppressed by m /F C , its contribution can be substantial for τ -lepton at large tan β as c R 2 is proportional to tan β (eq. (A.4)). Thus, BR(χ 0 1 → W ) grows at large tan β through the enhancement of τ W decay width. Another important aspect to note from eq. (A.4) is that both c L i and c R i are inversely proportional to µ. Thus for small values of µ, the decay rate Γ(χ 0 1 → W ) is much larger than νZ and νh decay modes. The W -decay modes becomes subdominant, as expected, for increasing value of µ. For µ = 300 GeV, the νh decay mode is substantial (> 60%) at low tan β and M 2 ∼ 700 GeV. For µ > 300 GeV and M 2 < 800 GeV we always have at least 50% branching ratio in the νh-mode irrespective of tan β choices.
From figure 1 we conclude that for large values of tan β and M 2 > 800 GeV, the χ 0 1 → νh is suppressed. This can be understood from the expression of C R j (eqn 16) which starts dominating when tan β and M 2 are large. To have better understanding of the behavior of BR(χ 0 1 → νh), we plot the BR in the plane of (M 2 , µ) for three different values of tan β = 5 (left), 15 (center) and 25 (right) in figure 2. One can see from the figure 2, that the largest possible values of the BR(χ 0 1 → νh) can be achieved when tan β is small. Also when M 1 ∼ 1 TeV, it is difficult to achieve BR (χ 0 1 → νh) ∼ 80% for our choice of Higgsino mass parameter (200GeV ≤ |µ| ≤ 600GeV) even for small tan β. Figure 3. Feynman diagrams representing pair production of electroweakinos. Similar diagram for chargino pair production mediated via γ * and Z boson also exist.

JHEP12(2016)062
Finally, we also note that our choice of ξ, µ, M 1,2 and tan β which predict correct neutrino mass are also consistent with the constraints coming from flavor violating Z and W -decay, lepton flavor violating decay of muon and neutrinoless double beta decay [79,80]. The last two gives the most stringent constraint among others.

Light Higgsino production at the LHC
From the analysis of previous section, we conclude that the decay channelχ 0 1 → νh is significant in large part of the parameter space and, in fact, close to ∼ 70-75% in certain regions. Motivated by this observation, we consider a very interesting signature of Higgs boson pair-production at the LHC in RpV SUSY. To probe this yet unexplored parameter space through di-Higgs, we consider pair production of Higgsinos at the 14 TeV LHC run: The main contributions to these processes come from s-channel mediation of γ, Z and W ± bosons (see figure 3). The contributions from t-channel squark mediated processes are suppressed by heavy squark masses and can be ignored. Hence, the crosssections only depend on the masses of the electroweak gauginos and their couplings with γ, Z and W ± bosons. We have worked with leading order (LO) cross-sections. The effects of next-to-leading order (NLO) QCD corrections can be achieved by including multiplicative K factors which can give a 10-20% enhancement [81]. However, here we present a conservative estimate without multiplying by any K-factors. In figure 4, we show the cross-sections for these production channels as a function of µ at 14 TeV LHC. The cross-sections for the associated production are the largest followed by the chargino and neutralino pair productions. In ref. [82], the authors have studied the impact of NLO corrections on various kinematical distributions and concluded that the NLO corrections do not change features of distributions rather they tend to change the overall scale of distributions. As in our analysis, we only present normalized distributions, they are robust under any NLO corrections.
The chargino mainly decays via R-parity conserving coupling toχ ±  decays ( h, νW, Z) are suppressed compared to the R-parity conserving one due to the tiny RPV couplings. Thus, in most of the events, one ends up with a pair of LSP's which subsequently decays toχ 0 1 → νh via R-parity violating interactions. This naturally leads to pair production of Higgs and two invisible neutrinos which contribute to missing transverse energy (/ E T ). The signal we investigate at the LHC, therefore, consists of pp → hh+ / E T +X, where X stands for additional jets and/or unclustered particles but no isolated leptons.
In table 1, we show the branching ratios of di-Higgs decays in various channels. The dominant decay of the SM-like Higss to a pair of b quarks occurs with a branching ratio of 60%. Thus, the dominant branching ratio for di-Higgs decay is to hh → bbbb, which is ∼ 36%. But the backgrounds for 4b and 2b2τ final states are overwhelming. While the backgrounds for 4τ decay mode are moderate, it suffers from small τ detection efficiencies. As far as decays to electroweak gauge bosons are concerned, bbW W and bbZZ decaying to bb + − νν suffer from huge QCD tt pair backgrounds. The bbZZ * → bb + 4 leptons and JHEP12(2016)062 Table 2. Photon and b jet identification efficiencies and misidentification probabilities for charm quarks, τ , light jets to b jets and photons at the LHC. [89] bbZγ channels suffer from too small rates. So, this leaves us a very few options to choose from. Among all the di-Higgs decay channels, we focus on the possibility where one Higgs decays to a bb pair and the other to a di-photon pair. Thus the final signature for di-Higgs search of interest in this work is γγbb/ E T + X at the LHC. We generate the spectrum for the BRpV SUSY model using SARAH [83] and SPheno [84,85], 2 and then, use these spectrum files in the PYTHIA [86] to generate the events. We also use PYTHIA for the parton showering and hadronisation. We have used simple cone algorithm of PYCELL inside PYTHIA to form jets out of clusters of hadrons with a cone size of R = 0.4. For simulating a generic detector, we also include the appropriate Gaussian smearing of the energies of each objects (i.e., jets and photon) in an event using the following resolution function, For jets, we take a = 1 GeV, b = 0.8 GeV 1/2 , c = 0.05 while for photons the values are a = 0.35 GeV, b = 0.07 GeV 1/2 , c = 0.007 [87,88]. Here, E is in the units of GeV. The identification efficiencies for a true b-jet, photon and their respective mistagging probabilities have been given in table 2. The tagging and mistagging efficiencies are more or less consistent with the ATLAS experiment and have been taken from [89].
In the following, we will perform a detailed collider simulation for the di-Higgs production process and the corresponding backgrounds at the 14 TeV high luminosity LHC with the luminosity of 3 ab −1 in the BRpV model. We will also discuss the strategy and cuts required to suppress the backgrounds and enhance the signal-to-background ratio. In our simulation, we consider different values of Higgsino masses, ranging from 200 GeV to 600 GeV, for the signal reconstruction and analysis. However for the purpose of illustration, we consider the following benchmark point for the model where BR(χ 0 1 → νh) is maximized: µ = 300 GeV, tan β

γγbb channel
We perform the signal calculation pp → hh / E T +X → bbγγ / E T +X. Though the branching ratio of the γγbb channel is very small, it can reach the similar sensitivity of the bbW W channel in probing the di-Higgs signals because of the precise resolution of two photons JHEP12(2016)062 µ (GeV) Cut 1 (fb) Cut 2 (fb) Cut 3 (fb) Cut 4 (fb) Cut 5 (fb) Table 3. Effects of the cut flow (discussed in section 4.1) on the signal events. The BR(χ 0 1 → νh) is assumed to be 100% for all benchmark points. which leads to a very prominent peak around Higgs mass in M γγ distribution. The background processes for the γγbb channel are QCD bbγγ, bbh(→ γγ), γγh(→ bb) and multijet QCD backgrounds resulting from jets faking either as b-jets or photons like jjγγ with two fake b jets; bbjγ with j faking photon; bbjj with two fake photons; jjjγ with two fake b jets and one fake photon; jjjj with two fake b-jets and two photons; hjj with either two fake b jets or two fake photons; and hjγ with one fake photon. While estimating and generating multijet backgrounds, we consider the misidentified charm quarks separately from the light flavour jets because of the very different mistagging factors as given in table 2.
In the following we present our strategy to suppress the background and enhance the signal-to-background ratio. The acceptance and selection cuts applied in our analysis are as follows: • Identification cuts (Cut 1): 1. Accept events with two photons, 2 b-jets and missing energy, 2. Photons must have transverse momentum p γ T > 10 GeV and rapidity |η | < 2.5, 3. All b-jets must have following p T and η requirements: 4. All pairs of jets, photons and photon plus jets should be well separated with each other by: ∆R jj,jb,bb,γj,γb,γγ ≥ 0.4 where ∆R = (∆φ) 2 + (∆η) 2 .
• Selection requirements: when an event satisfies above requirements, it is further processed for the signal reconstruction and background reduction as follows: -Cut on p T of photons (Cut 2): to get rid of soft photons coming from the decay of mesons or radiations, we further put the following p T cuts on the two photons: p γ 1 T > 30 GeV and p γ 2 T > 20 GeV. In figure 5, we show the transverse momentum distribution of the two photons for the signal as well as for the backgrounds before the cuts. For the hardest photon, the distribution peaks at around 100 GeV while most backgrounds peak at low p T .

JHEP12(2016)062
-Invariant mass cuts (Cut 3): in figure 6, we show the invariant mass distribution of two photons, M γγ and two b jets, M bb for both the signal and backgrounds. In the case of the signal and some of the backgrounds, the two photons come from a resonant Higgs. Thus, the diphoton invariant mass distributions for the resonant Higgs have a very narrow peak around M h = 125 GeV. On the other hand, for other backgrounds, the distribution is continuum with no localised events around M h . Also, the fact that the photons are measured at a very high degree of precision at the LHC, allows a very sharp distribution even at the detector level. On the other hand, due to the large jet energy uncertainties and energy resolution of the jets, the invariant mass distributions of the resonant Higgs decaying to two b-jets show a relatively-wider peak. The peak is also shifted towards lower value (around 112.5 GeV) than the actual Higgs mass due to invisible neutrinos coming from b-decays and missing particles outside the jet cone. Both the diphoton and di-b-jet invariant mass distributions are quite distinct from the background and have prominent peaks. We utilize this to suppress the backgrounds and further cuts on invariant masses are imposed: -Missing transverse Energy, / E T distribution (Cut 4): in this analysis, the di-Higgs signal arise from the decays of two lightest neutralinos to Higgs and neutrino. The neutrinos coming from heavy particles are expected to have a large transverse momentum contributing to / E T . In figure 7, we show the / E T distribution for the signal as well as backgrounds. For the M χ 0 1 = 300 GeV, the / E T distribution is expected to peak at around 200 GeV. For the non-top backgrounds, the only source of / E T comes from uncertainties in the measurements of p T of the various objects at the detector, dominant being jets. Thus, the peaks in / E T distributions for such backgrounds are at very low / E T . On the other hand, the backgrounds, which contain top-pairs, produce neutrinos when they decay semi-leptonically and thus may have large / E T . From figure 7 we see that the / E T distribution for a non-top background almost vanish at (≤ 80 GeV) while the backgrounds containing a top-pair have significant / E T until / E T < 200 GeV. For the signal, the peak is at 200 GeV and the distribution remains significant until a very large value of / E T ∼ 800 GeV. Thus, we further put a cut of / E T > 100 GeV to further suppress the backgrounds. With this cut, the bbγγ background is completely eliminated and the total background is suppressed by a factor of 30 while the signal events are affected only by 20%.
-∆R separation (Cut 5): we find two interesting angular correlations which are significantly different between signal and backgrounds. The ∆R separation of the photon-photon and bb pairs are small for the signals and larger for the backgrounds. The shape of the signal distribution stems from the fact that the  photon pair and b pair arise from the highly boosted resonant Higgs bosons. These Higgs bosons are boosted as they are produced from the decays of heavy neutralinos/charginos. On the other hand, diphotons and di-b jets for the backgrounds do not come from any heavy resonance and thus are expected to be farther spaced. This fact can be seen from figure 8 where we show the ∆R distribution for γγ and bb pairs. A cut on ∆R can effectively suppress the background relative to the signal. Based on these observations, we put following additional cuts: ∆R bb < 2.0, ∆R γγ < 2.0.
These cuts cuts reduce the background by an order of magnitude while the signal events are reduced by only 20%.  In table 3, we show the cut flow of the cross-sections for the signal for different values of the Higgsino masses. One can observe that the cross-section rapidly decreases with the increase in masses but the efficiency of the cuts is better for the heavier Higgsinos. In table 4, we show the cut flow of the cross-sections for the various background processes considered in the analysis. One can see from table 4, the dominant contributions to the background come from the bbγγ continuum, ttγγ and tth processes, each having crosssections of 44 fb, 1.2 fb and 0.11 fb respectively after cut 1. On the other hand, the signal cross-section after cut 1 is 0.096 fb for Mχ0 = 300 GeV. Thus the signal-to-background ratio at this stage of cut flow is only 5 × 10    ttγγ and tth are reduced to 10 −4 and 10 −3 respectively. After the cut 4, the S/B ratio is 6 which is a tremendous improvement over previous value after cut 1. The cut 5 (on ∆R) further suppresses the background contribution and thus enhancing the S/B ratio. We estimate the statistical significance of the signal using following formula [90] Sig = 2 (S + B) ln 1 + S B − S (4.1) where S and B are the number of signal and background events, respectively. We will now comment on the possibility of detecting the di-Higgs signal for different masses of the Higgsinos. In table 5, we show the cut efficiency of the signal for different masses of Higgsinos. As the masses increase, the total production cross-section goes down rapidly, but the handle over the background improves significantly. Thus, even though the cross-section for the signal is decreased with the mass, the efficiency for the detection of the signal is increased.
Let us now look at the various distributions discussed earlier, for the different values of Higgsino masses. In figure 9, we show the missing transverse energy ( / E T ) distribution for the di-Higgs signal for different values of the Higgsino masses ranging from 200 GeV to 600 GeV. We notice that the / E T distribution peaks at 60 GeV, 100 GeV, 150 GeV, 200 GeV and 250 GeV for M χ ± = 200 GeV, 300 GeV, 400 GeV, 500 GeV and 600 GeV, respectively. Thus, in order to further improve the cut efficiency, optimizing of the cuts on / E T should be performed for different Higgsino masses. However, one should keep in mind that as the backgrounds are already low, any further stringent cuts would also result in reducing the signal events.

JHEP12(2016)062
In figure 10, the ∆R separation between the photon-photon (left) and bb pairs (right) for the di-Higgs+/ E T signal at the 14 TeV LHC have been shown. For very heavy Higgsinos, the Higgs coming from the decay will be highly boosted. Thus the ∆R separation between bb/γγ becomes very small. This fact determines the shape of the ∆R distributions for different chargino masses. For the mass of 600 GeV, the distribution is peaked at a very low value of ∆R while for a light chargino of 200 GeV, the diphotons and di-b jets are much farther spaced. Based on these findings, the cuts on ∆R need to be judiciously chosen so as to enhance the signal-to-background ratio for the different values of Higgsino masses.
We now discuss the exclusion/discovery limits for the Higgsino LSP in di-Higgs production at the LHC14 for the different values of integrated luminosities. Based on the cut-flow analysis for the signal and backgrounds, we estimate the significance of the signal, using eq. (4.1), for various values of µ and for three different values of integrated luminosities 1, 2 and 3 ab −1 . The numbers for the significance have been presented in table 6. For a light LSP of Mχ0 1 ∼ 400 GeV, it would be possible to discover them with more than 5σ significance at 1 ab −1 of integrated luminosity. On the other hand, for a heavy LSP Mχ0 1 ∼ 600 GeV, it will need around 3 ab −1 of data to have sufficient discovery prospects.
In figure 11 we show the statistical significance for the signal in the plane of (µ, tan β) for M 1 = M 2 = 700 at integrated luminosity of 3 ab −1 . For these values of M 1 and M 2 , we find that the BR is maximized and di-Higgs production has better prospects. As mentioned earlier, Br(χ 0 1 → νh) also depends on parameters tan β and µ. This is also evident from figure 2. The dark regions in the figure denote the regions with high signal JHEP12(2016)062 significance denoting better prospects. For µ = 200 GeV, we find that it would be possible to discover the LSP in di-Higgs production process at the LHC with only 500 fb −1 of integrated luminosity only for small values of tan β < 10. For small µ, the BR to νh drops quickly at larger values of tan β while for large µ ∼ 500 GeV, the BR fairly remains constant in the range %. Thus we find that the signal significance for µ = 500 GeV is considerably large even for large tan β up to 40.

Summary and conclusions
Motivated by the naturalness, we study a Bilinear RpV SUSY scenario with µ below TeV scale and the LSP has substantial Higgsino component, and BRpV couplings determine the tree-level neutrino mass matrix. We investigate the parameter space of this scenario and study the decay patterns of the LSP. We find that in a large part of the parameter space, the LSP decays to νh with branching fraction larger than 0.7. A large region of this yet unexplored parameter space can still be insensitive to the existing constraints coming from the LHC searches. We then study the pair production of the electroweakinos followed by the decaysχ ± 1 →χ 0 1 W ± * andχ 0 1 → νh. This leads to a very interesting signature of Higgs boson pair production at the LHC accompanied with significant missing transverse energy. This di-Higgs production in BRpV model, occurring in association with a large missing transverse energy, is in distinct contrast to the Higgs boson pair production in the SM. This fact makes this signal quite feasible to search at the luminosity improved version of the 14 TeV LHC despite having a very small cross-section.
Among the various decay channels for di-Higgs, we focus on the scenario where one of the Higgses decays to diphoton and the other decays to a bottom pair. This particular decay has the advantage of manageable SM backgrounds. Thus the signal which we are looking for includes 2 photons, 2 b jets and large missing transverse energy. We perform a realistic detector level simulation for the signal γγbb + / E T taking some benchmark points in the parameter space at the 14 TeV LHC. We also perform a full systematic study of all the background processes. It is found that the cuts on / E T and ∆R are instrumental in eliminating the QCD multijet backgrounds and suppressing the total background. We also notice that even though the cross-sections for the signal decrease as the masses of the Higgsinos get heavy, the increased efficiency of the / E T and ∆R cuts helps to compensate the overall signal to background ratio. Finally we conclude that the LSP of mass 300-500 GeV would be amenable to discovery in the early LHC14 data in the di-Higgs channel. On the other hand for the heavy LSP ∼ 600 GeV, it would require full data set (3 ab −1 ) of 14 TeV run to have reasonable discovery prospects.

A Effective R-parity violating vertices
Let us summarize all the relevant RpV diagonalization matrices and vertices confirming the results in refs. [31-33, 36, 37].
Rotating away the neutrino-neutralino mixing mass terms (by θ N ) can be made by the following redefinition of neutrinos and neutralinos: where (ν i ) and (χ 0 j ) represent three neutrinos (ν e , ν µ , ν τ ) and four neutralinos (B,W 3 , H 0 1 ,H 0 2 ) in the flavor basis, respectively. The rotation elements θ N ij are given by Here s W = sin θ W and c W = cos θ W with the weak mixing angle θ W .

Charged lepton/chargino diagonalization.
Defining θ L and θ R as the two rotation matrices corresponding to the left-handed negatively and positively charged fermions, we have where e i and e c i denote the left-handed charged leptons and anti-leptons, (χ − j ) = (W − ,H − 1 ) and (χ + j ) = (W + ,H + 2 ). The rotation elements θ L,R ij are given by We find that the above expression for θ R differs from the result by Nowakowski et al., in refs.  and in refs. [34,35].

Sneutrino/neutral Higgs boson diagonalization.
Denoting the rotation matrix by θ S i = a i , we get (A.5)

JHEP12(2016)062
With the expressions for the rotation matrices, we can obtain the effective R-parity violating vertices from the usual R-parity conserving interaction vertices, which are relevant to the LSP decays. We list them below by taking only the linear terms in θ's which are enough for our purpose.

JHEP12(2016)062
Let us remark that the RpV vertices of the neutralinos and charginos depend only on ξ i although the RpV rotation matrices depend also on i or a i . Contrary to eq. (19), eqs. (18), (20) disagree with the results of refs. [34,35].

B Decay widths of neutralinos
For a generic decay processχ i → L j V where L is either ν or ± and V is either Z or W ± , the decay width can be written as: where C L i and C R i are the left-and right-handed couplings, r V is (m 2 V /m 2 χ ) and I 2 (r V ) = (1 − r V ) 2 × (1 + 2 r V ).

Forχ
Similarly, for the decayχ 0 i → ν j h, the decay width can be written as where r h is m 2 h /m 2 χ .
Here, N is the 4 × 4 matrix and diagonalize the neutralino mass matrices.
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.