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 supersymmetric particle (LSP). We observe that the LSP dominantly decays to $\nu h$ in a large part of the parameter space, and thus study the pair production of electroweakinos followed by the decays $\tilde\chi^\pm_1\to \tilde\chi^0_1 W^{\pm*}$ and $\tilde\chi^0_1\to \nu h$. 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 realistic detector level simulation 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 the best framework to protect the electroweak scale from the quadratic divergence caused by a certain ultra-violet (UV) physics. As a bonus, supersymmetry provides a good candidate of dark matter of the Universe: the lightest supersymmetric particle (LSP) which is neutral and stablized by assuming R-parity conservation. On the other hand, R-parity violation (RpV) brings another interesting possibility of generating tiny neutrino masses [1,2] 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 ∓ [3,4,5,6]. 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) requires a fine cancellation among different terms. Barring too huge cancellation, one may arrange m H u,d and µ not too larger than m Z [7], which still remains a viable option for SUSY. It is a challenge for the LHC and future colliders to probe such a degenerate electroweakino [8].
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 Higgsino 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 [9]. 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 [10,19]. Di-Higgs production at the LHC has been studied in the context of triliniear Higgs self coupling measurements by various authors [11] and references there in. While the SM cross sections are too small to allow observation 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 γγbbE / 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 [12,13], double Higgs production via gluon fusion in the effective field theory framework [14], Higgs pair production in the context of SUSY extension of the SM [15,16,17] and various other extensions of it [18].
This paper is organized as follows: In Section 2 summarizing the results of Refs. [4,6], we provide a brief review on the bilinear RpV couplings constrained by the resulting tree-level neutrino mass matrix. In Section 3, the Higgsino-like LSP decay modes are analyzed for the cases of µ < M 1 ≤ M 2 where M 1,2 are the masses of the bino and wino components, respectively. In Section 4, we compute the Higgsino pair production cross-section for some benchmark points and perform a realistic detector-level simulation for signals 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.

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 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 (6) 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 [20]. 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 (6) 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 a specific scenario in this work where µ < M 1 , M 2 leading to a situation of the Higgsino LSP, with the lightest statesχ 0 1,2 andχ ± 1 being Higgsino-like. If M 1 < M 2 ,χ 0 3 is bino-like whileχ ± 2 andχ 0 4 are wino-like; and for M 2 < M 1 ,χ 0 4 is bino-like whileχ ± 2 andχ 0 3 are wino-like. In the following, we fix M 1 = 1.0 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 β (Fig. 1). Fig. 1 shows that the branching ratio for the W or νh decay mode is quite sensitive to the values of µ, M 2 and tan β.
All of the RpV vertices which we have obtained in Eqs. 18-23 of Appendix A depend only on one type of BRpV variables ξ i which have taken to satisfy ξ 1 : ξ 2 : ξ 3 = 0.1 : 1 : 1 (Eq. 10). 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. 26). 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. 16). Thus, BR(χ 0 1 → W ) grows at large tan β through the enhancement of τ W decay width.
Another important aspect to note from Eq. 16 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. This feature can be seen from the top row of Fig. 1 where we have displayed the BRs of W , νh and νZ decay modes for µ = 200 GeV. In these figures, it can be seen that BR( W ) dominates at large tan β. Decay mode νh is only significant at low tan β. For larger value of µ = 500 GeV (see bottom panel of Fig. 1), BR( W ) becomes subdominant as expected. Thus, the dominant decay channel with the branching ratio varying between (0.6-0.9) is the νh mode for all values of tan β. BR(νh) decreases by about 20% because of the increase in the decay rate of χ 0 1 → W at large tan β.
From the Fig. 1, we conclude that for large M 2 = 3 TeV, the BR(χ 0 1 → νh) always dominates over all otherχ 0 1 decays. For moderate M 2 = 2 TeV, the BR(χ 0 1 → νh) still dominates over W and νZ modes but for low values of tan β while for small M 2 = 1 TeV, we find that BR(χ 0 1 → νh) dominates but only for large value of tan β and for large values of µ. All these facts can easily be deduced from the effective RpV couplings which are derived in Appendix A. 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)

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 ∼ 90 − 95% 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: 3. Associated neutralino and chargino production : The main contributions to these processes come from s-channel mediation of γ, Z and W ± bosons (see Fig. 3). The contributions from t-channel squark mediated processes are suppressed by heavy squark masses and can be ignored. Hence, the cross sections 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. However, here we present a conservative estimate without multiplying by any K-factors. In Fig. 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.q Figure 3: Feynman diagrams representing pair production of electroweakinos. Similar diagram for chargino pair production mediated via γ * and Z boson also exist.
The chargino mainly decays via R-parity conserving coupling toχ ± 1 →χ 0 1 W ± * → f f χ 0 1 almost 90% of time in all regions of the parameter space. The R-parity violating 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.  Figure 4: Production cross sections of the charged and neutral Higgsinos in terms of µ at the LHC14.
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 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 γγbbE / T + X at the LHC.
We generate the spectrum for the BRpV SUSY model using SARAH [21] and SPheno [22], and then, use these spectrum files in the PYTHIA [23] 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 a realistic detector simulation, 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 [24,25]. 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 [26].
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 β = 15, M 2 = 3 TeV. For this benchmark point, we have Mχ0 1 = 300 GeV, Mχ± 1 = 310 GeV, BR(χ 0 1 → νh) = 0.9, BR(χ + 1 →χ 0 1 W + ) = 0.95 and σ(pp → χχ) = 232 fb. We have shown all the figures for the various distribution for this benchmark point. Nevertheless we also present our analysis for various Higgsino masses in Table 3.

γγ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 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 > 30GeV and p γ 2 T > 20GeV. In Fig. 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 .
-Invariant mass cuts (Cut 3): In Fig. 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 relativelywider 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 Fig. 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 Fig. 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 > 100GeV 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 Fig. 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: 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 cross sections 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 −4 . The cut on missing transverse energy E / T > 100 GeV (cut 4) almost eliminates the bbγγ background while the contributions of 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 [27] Cut 1 (fb) Cut 2 (fb) Cut 3 (fb) Cut 4 (fb) Cut 5 (fb) 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 Fig. 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. In Fig. 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. 12, 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 Fig. 11 we show the statistical significance for the signal in the plane of (M 2 , tan β) for two different values of µ = 300 (left) and 500 GeV (right). As mentioned earlier, Br(χ 0 1 → νh) also depends on parameters tan β and M 2 other than µ. The partial width in individual decay mode of the LSP depends on ratio of M 2 /M 1 and tan β on such a way that the branching ratio to νh is minimum when M 2 /M 1 ≈ 1. Strictly speaking this depends on tan β as well. This is also evident from Fig. 2. The dark regions in the figures denote the regions for low discovery prospects while the brighter region has better prospects to observe a Higgsino LSP in di-Higgs process. For µ = 300 GeV, we find that it would be possible to discover the Higgsino-like LSP in di-Higgs production process at the LHC with only 300 fb −1 of integrated luminosity and a large portion of the parameter can be excluded at 2σ with 1 ab −1 . The region of exclusion has been shown by red curves in the figures. The 5σ discovery prospect has been denoted by green curves in both the plots. For µ = 500 GeV, the 5σ discovery can only be obtained with full data set of the 14 TeV LHC viz. 3 ab −1 for large values of M 2 .

Summary and Conclusions
Motivated by the naturalness, we study a Bilinear RpV SUSY scenario where the LSP is Higgsino-like, 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 Higgsino LSP. We find that in a large part of the parameter space, the LSP decays to νh with branching fraction larger than 0.9. 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.
Acknowledgment: EJC is supported by the NRF grant funded by the Korea government (MSIP) (No. 2009-0083526) through KNRC at Seoul National University.

Neutrino-neutralino diagonalization
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 and F C = M 2 + M 2 W s 2β /µ. Sneutrino/neutral Higgs boson diagonalization Denoting the rotation matrix by θ S i = a i , we get 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.

Appendix 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.