Bilinear R-parity violating supersymmetry under the light of neutrino oscillation, higgs and flavor data

In this work, we explore a well motivated beyond the Standard Model scenario, namely, R-parity violating Supersymmetry, in the context of light neutrino masses and mixing. We assume that the R-parity is only broken by the lepton number violating bilinear term. We try to fit two non-zero neutrino mass square differences and three mixing angle values obtained from the global $\chi^2$ analysis of neutrino oscillation data. We have also taken into account the updated data of the standard model (SM) Higgs mass and its coupling strengths with other SM particles from LHC Run-II along with low energy flavor violating constraints like rare b-hadron decays. We have used a Markov Chain Monte Carlo (MCMC) analysis to constrain the new physics parameter space. While doing so, we ensure that all the existing collider constraints are duly taken into account. Through our analysis, we have derived the most stringent constraints possible to date with existing data on the 9 bilinear R-parity violating parameters along with $\mu$ and $\tan\beta$. We further explore the possibility of explaining the anomalous muon~(g~-~2) measurement staying within the parameter space allowed by neutrino, Higgs and flavor data while satisfying the collider constraints as well. We find that there still remains a small sub-TeV parameter space where the required excess can be obtained.


Introduction
Neutrino oscillation is one of the most robust indications towards the existence of physics beyond the standard model (BSM).Over the years, multiple experiments have been studying the neutrino oscillation phenomena, see e.g., [1][2][3][4][5][6][7].Their measurement of two mass square differences and three mixing angles imply significant mixing among the three light neutrino states of which at least two must have non-zero masses [8].The standard model (SM) [9][10][11][12] or the R-parity 1 conserving minimal supersymmetric standard model (MSSM) [13][14][15] cannot address the neutrino oscillation phenomena.The light neutrino masses and mixing can be generated by simple seesaw extensions of the SM, which are different manifestations of the dimension- 5 Weinberg operator [16,17].In Type-I seesaw, we add right-handed singlet fermions in the model, in Type-II seesaw, we add fermionic triplets and in Type-III seesaw the objective is achieved by adding scalar triplets to SM [18][19][20][21][22].In R-parity violating (RPV) MSSM [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42] scenario, one can explain neutrino oscillation phenomena without incorporating the Weinberg operator.In the light of the updated neutrino oscillation data and other relevant constraints, it is worth revisiting the scenario to gauge their impact on the RPV couplings.R-parity conserving MSSM is more widely studied in literature because it offers a natural dark matter candidate in the form of the lightest supersymmetric particle (LSP) which cannot decay further and is therefore stable2 .However, one has to incorporate R-parity conservation by hand to achieve that.The symmetry principles to write the Lagrangian allow us to add four RPV terms in the superpotential as below: The bilinear term, ϵ i Li Ĥu , and the next two trilinear terms containing λ, λ ′ in the Eq.1.1 each violates lepton number by one unit and the last λ ′′ term violates the baryon number by one unit.Here Li ( Êk ) corresponds to the left (right) handed lepton supermultiplet and Ĥu is the up-type Higgs supermultiplet.Qj , Ûj ( Dk ) represent left-handed doublet and right-handed singlet up-type (down-type) quark supermultiplet respectively.One can generate non-zero light neutrino masses through the trilinear λ or λ ′ couplings at one-loop [45,46].Note that, with these couplings, the neutrinos are still massless at tree level.The bilinear term is capable of generating one neutrino mass at tree level [35,45,46].However, one also needs to take into account the one-loop contributions to explain the oscillation data.The trilinear couplings need not be non-zero in that case from the perspective of light neutrino mass generation.Note that, the bilinear RPV (bRPV) terms can exist even in the absence of the trilinear terms in the theory and the trilinear couplings can be generated starting from the bilinear couplings [47].On the other hand, if one starts from only trilinear RPV terms, the bilinear RPV couplings can be generated through renormalisation group evolution at a different energy scale [48,49].Understandably, neutrino masses and mixing angles lead to constraints on the trilinear couplings [50].Bilinear RPV can be assumed to be the fundamental theory and hence, in this work, we only focus on non-zero values of these couplings keeping all trilinear RPV couplings zero.Note that, in the alignment of bilinear coupling parameters (ϵ i ) and sneutrino vev (v i ), which arises naturally in the framework of horizontal symmetries, the three light mass eigenstates of L i correspond to the three light neutrinos [25,51].To achieve this, the alignment between soft coupling parameters, B i and ϵ i is also required i.e., B i ∝ ϵ i [25,51].In such cases, both the bilinear term and the soft breaking bilinear term can be rotated away by the field redefinition of L i and H u .
Although the contribution of bilinear RPV couplings towards neutrino masses and mixings has been studied in the past [52][53][54][55], a detailed statistical analysis that can highlight the allowed parameter space is missing in the existing literature.There has been some effort to constrain the bilinear RPV couplings from neutrino physics perspective [52][53][54][55][56][57][58][59][60] but they are either not very generalized or simply inadequate in the light of new oscillation data [8].Moreover, adding RPV couplings leads to a wide range of phenomenological implications [24].One needs to take into account the modified bounds on the SUSY particles which can now decay exclusively into SM particles.Because of the presence of the bilinear RPV term, the neutral and charged Higgs states now can mix with the sneutrinos and charged sleptons respectively.Therefore, one needs to carefully check the SM Higgs coupling strengths in the light of the updated dataset [61,62] which can further put constraints on the RPV parameters.In this study, we take into account all these possibilities.New measurement of muon magnetic moment at Fermilab has slightly changed the existing world average, which shows a 4.2σ deviation 3 at present from SM predic-tion [67,68].
Sneutrino-chargino and slepton-neutralino loops have been studied extensively in this context and it is evident that given the present collider constraints, it is quite difficult to achieve the required excess within the framework of the MSSM.One region of sub-TeV allowed parameter space still relevant in this context is the compressed region where LSP-NLSP mass difference is very small [69].These kinds of compressed regions are difficult to probe owing to the poor detection prospect of the final state particles.In the RPV context, however, these particles can decay further into SM particles if the RPV couplings are large enough for a prompt decay.That leads to new collider constraints obtained from direct searches at the LHC [70][71][72][73][74][75][76][77][78].Having said that, there are some additional contributions to the muon (g -2) in the RPV framework owing to the mixing between charged higgs-sleptons, neutral higgs-sneutrinos, charginos-charged leptons and neutralino-neutrinos.In the present context, it is worth a look since RPV couplings and some relevant particle masses, e.g., that of sneutrinos and neutralinos have a big role to play in both muon (g -2) calculation and generation of light neutrino masses and mixing.We, therefore, explore the possibility of explaining the muon (g -2) excess over the SM contribution within this framework while simultaneously satisfying all other experimental constraints.
The presence of the bilinear RPV term results in a lepton number violation by one unit.A sneutrino state therefore can now acquire non-zero vacuum expectation value (VEV).All three sneutrino VEVs along with the three ϵ i parameter and their corresponding soft terms are crucial in fitting the neutrino oscillation data.Moreover, because of the sneutrino and neutral Higgs mixing, the Higgs sector parameters like trilinear stop coupling (A t ), the ratio of up and down type Higgs vacuum expectation values (tanβ) and the µ parameter are also relevant to our objective.This results in a substantially bigger set of unknown input parameters and as a result, a conventional random scan does not produce the coveted results.We, therefore, adopt a Markov Chain Monte Carlo (MCMC) algorithm to sample the parameter space.The main objective of this work is to take advantage of the current observational data in order to constrain the model parameters by a robust statistical analysis.The method described here is not only able to enhance our understanding of the importance of currently available data on the modeling of R-parity violating supersymmetry (RPV SUSY), but also explores the ability to use this technique for upcoming experiments.We concentrate mostly on the neutrino and Higgs sector observables apart from some relevant flavor observables to locate the favored parameter region.While doing so, we ensure that the collider constraints on supersymmetric parameters are also taken (ETMC) [65] and the preliminary results of the CDM-3 detector [66] indicate that the discrepancy between the observed and predicted values of muon (g -2) will be smaller and less significant.
into account.We have kept the parameter space as generalized as possible albeit with some simplified assumptions on the soft masses of charged sleptons, sneutrinos and squarks.
The plan/structure of this paper is as follows.In Sec.2 we briefly discuss the generation of light neutrino masses and mixings in the context of bRPV SUSY model.In Sec.3 we first mention the constraints coming from the global analysis of neutrino oscillation data, then Higgs mass and the constraints from Higgs signal strength measurements.In addition, we have considered the flavor constraints arising from rare B decays.We also briefly discuss the prescription of the MCMC technique adopted by us to identify the favored parameter region.Sec.4 contains the results for normal and inverted hierarchy scenarios.We briefly discuss the prospects of addressing the anomalous muon (g -2) issue in Sec.5 and finally, we present our conclusion in Sec.6.

Neutrino Mass and mixing from Bilinear RPV SUSY
In the bRPV SUSY model the additional RPV term in the superpotential generates the mixing between the neutrinos and neutralinos [23,24,46] and is written as: Here ϵ i (i = 1, 2, 3) represents the bRPV mass parameters.The Lagrangian corresponding to the superpotential gives rise to mixing between up-type Higgsino ( H0 u ) and three light neutrinos (ν iL ).It also generates mixing between H+ u and left-handed leptons (l iL ).The Lagrangian of the superpotential and the soft term are given by where B i corresponds to the soft bRPV coupling parameter representing the coupling between sneutrino and the neutral Higgs bosons, and Li is the left-handed slepton multiplet.So, at the tree level, the resultant 7 × 7 neutralino-neutrino mass matrix in the basis of ψ 0 = B W3 H0 d H0 u ν e ν µ ν τ looks like [25,26,45,49] as where, B ( W3 ) denotes bino (wino) and M 1 (M 2 ) is bino (wino) mass parameter.v u , v d are the VEVs for up-type and down-type Higgs respectively.The sneutrino VEVs are represented by ⟨ν i ⟩ ≡ v i (i = 1, 2, 3).Diagonalising the above mass matrix gives rise to one non-zero light neutrino mass apart from four massive neutralinos.However, at least one other light neutrino must be massive in order to satisfy the neutrino oscillation data.This is achieved at one-loop level.
i The tree level mass of one of the neutrinos also receives some loop correction and understandably this one happens to be the heaviest of the three neutrinos.The admixture of the three massive neutrinos in the respective mass eigenstates of course depends on the choice of hierarchy.Tree level contribution, [m ν ] ϵϵ ij , involves ϵ i term, which indicates the mixing between the up-type Higgsino and the neutrino (see Fig. 1).There are two different loop contributions to neutrino masses from this model, namely, the BB loop, and the ϵB loop [27,28,45].Combining all the contributions we can write neutrino mass matrix as [24,45] [ correspond to the tree level contribution, the BB loop contribution and ϵB loop contribution respectively.Now if we look at the individual contribution, for the tree level, ζ represents the alignment between ϵ i and v i [29,49,79,80].For different basis choices, the alignment will be different 4 .In the tree level contribution, X T is defined as [45,46] where m γ ≡ cos 2 θ w M 1 + sin 2 θ w M 2 .In bRPV models, the neutral Higgses and the sneutrinos mix at the tree level via B i parameters.This leads to a finite mass splitting between the CP even and CP odd sneutrino mass eigenstates.This mass splitting of sneutrinos is responsible for the generation of Majorana neutrino mass at the one-loop level (see Fig. 2(a)).For a detailed discussion on the cancellation between different Higgs (h, H, A) mediated BB loop diagrams and the effect of sneutrino degeneracy from the BB loop see Refs [27,45,46].The ϵB loop diagram involves both the bilinear term ϵ i and B j .The combination of higgsino-neutrino mixing and Higgs-sneutrino mixing give rise to neutrino mass (see Fig. 2(b)).The ϵB loop contribution is subleading to the BB loop.Now if we consider that all the masses are at electroweak breaking scale ( m), then the approximate contributions to the neutrino mass matrix (mentioned in Eq.2.4) is given by [46]: where ϵ H or ϵ ′ H arises due to the cancellation of the different Higgs (h, H, A) diagram in BB and ϵB loop respectively and depending upon the parameter space, they can suppress the mass contribution by several orders [45].A recent article [82] nicely summarises the different contributions to the light neutrino mass matrix for different model scenarios in the RPV context.The equations 2.6 -2.8 are consistent with the ones quoted in the article with a small difference, that is the assumption of absence of sneutrino VEV (v i = 0) taken in [82].
As already mentioned, only one of the three light neutrinos becomes massive at the tree level.If we consider the Normal Hierarchy (NH) scenario where the ordering of neutrino mass is m ν 3 > m ν 2 > m ν 1 , then it is evident from Eq. 2.6 that the heaviest neutrino mass is proportional to cos 2 β(ϵ 2 1 + ϵ 2 2 + ϵ 2 3 ) i.e., tanβ acts as a suppression factor for m ν 3 .The masses of the other two neutrinos, m ν 2 and m ν 1 , are generated at the one-loop level where the dominant contributions come from the BB loop.In general, B i and ϵ i are not related to each other and the leading contributions to m ν 2 comes from the BB loop with an enhancement effect from the 1 cos 2 β part as shown in Eq 2.7.The lightest neutrino (m ν 1 ) in the NH scenario also can be massive from this same loop contribution.When the sneutrinos of different generations are non-degenerate, m ν 1 is proportional to the square of the sneutrino mass splitting between different generations.Hence to get the neutrino mass square differences and mixing angles in the existing ranges obtained from various neutrino oscillation experiments, we need tree level contribution as well as the loop contribution.This puts a restriction, among other parameters, on the choice of tan β which cannot be either very large or very small.It may be noted that for the Inverted Hierarchy (IH) scenario where m ν 2 > m ν 1 > m ν 3 , the relations will change accordingly.

Computational Set-up and Numerical Constraints
In this section, we first summarize the numerical constraints that have been used in this study.The neutrino observables obtained from the latest global fit of different neutrino oscillation data, the most updated measurement of Higgs mass and its coupling strength in different decay modes along with low energy data from rare bdecays are considered in this analysis.Relevant limits derived from the LHC Run-I and Run-II data are also summarised in this section.Finally, we discuss the range of parameter space considered for scanning and summarize briefly the likelihood analysis implemented using the Markov Chain Monte Carlo (MCMC) algorithm.

Constraints from Neutrino observables
Global analysis of neutrino oscillation data provides us with two mass-square differences and three mixing angles.The mass-squared differences are defined as where, the m i (i = 1, 2, 3) represent the physical masses of the three light neutrinos.The sign of ∆m 2 31 (or ∆m 2  32 ) remains unknown to date, which gives rise to the two hierarchial (NH and IH) scenarios.The three relevant mixing angles between different generations are represented by θ 12 , θ 13 and θ 23 .The mixing angles can be calculated from the 3 × 3 light neutrino mass matrix, also known as PMNS matrix [83][84][85] [8,86].For this analysis, we have used the best-fit values and 1σ ranges of neutrino oscillation parameters obtained by the updated global fit in Ref. [8].The best-fit points along with 1σ range of these parameters in NH and IH scenarios are summarized in Table 1 [8].The δ CP corresponding to NH scenario is quite consistent with π well within 1σ, which is consistent with zero CP-violation in the neutrino sector.The existing uncertainties on this parameter are also much bigger compared to the mixing angles and mass-squared differences.Adding the δ CP parameter does not restrict the parameter space any further and hence we choose to keep δ CP = 0 throughout our analysis.

Constraints from Collider experiments
While fitting the neutrino physics observables, one also needs to ensure that the fitted model particle spectra obey other experimental constraints.The bRPV term gives rise to mixing among the sneutrino and Higgs states of the MSSM whereas the charged sleptons now mix with the charged Higgs sectors.In addition to that, to fit the neutrino data we are also varying µ and tan β, which impacts the mass of the lightest CP even Higgs boson 5 and its coupling strengths with SM particles.These are quite precisely measured and as a result, restrict the choices of the parameters affecting them.Apart from the neutrino and Higgs sector, we also consider constraints arising from the branching of rare B-hadron decays such as BR(B → X s γ) and BR(B s → µ + µ − ).Last but not least, the existing constraints on SUSY particles from various direct searches are also duly taken into account.

Constraints from Higgs Sector
The measured mass of SM-like Higgs boson obtained from the combined data of the ATLAS and CMS experiments is 125.09± 0.21(stat.)± 0.11(syst.)GeV [89].Taking into account the theoretical uncertainty of the Higgs mass calculation within the SUSY framework, we consider ±3 GeV window for Higgs mass around the bestfit value [90].Apart from the mass, the signal strengths of 125 GeV Higgs are also precisely measured by both the CMS and ATLAS collaborations [61,62].The updated results of the coupling strength modifiers (κ i ), i.e., BSM over SM ratios of the coupling strengths for a particular decay mode i, along with their 1σ uncertainties obtained by CMS collaboration from LHC Run-II data with luminosity L = 137 f b −1 are summarized in Table 2.
Note that, although some of the data points have very similar best-fit and uncertainty ranges (e.g., κ z , κ τ ), there are considerable differences in some other measurements (e.g., κ b , κ µ ) obtained from the ATLAS collaboration [62].These differences  especially in the best-fit points with similar uncertainty can slightly change the favored parameter space.Hence, we cross-check our results using ATLAS data [62] as well and comment on how much change is expected.The coupling strengths in the present model framework have been computed using SPheno which does a full two-loop calculation for the Higgs sector.Note that the CMS collaboration also quotes their measurement of the Higgs-gluon-gluon effective coupling strength as κ g = 1.16 +0.12 −0.11 , which we have not included in our analysis.This is because the effective Higgs-gluon-gluon coupling strength is quite sensitive to the choice of the SUSY breaking scale.Hence one would in principle have to vary the scale as well as a parameter, which affects among other things, the 125 GeV Higgs mass itself.That in turn prevents us from fixing some Higgs sector parameters, such as A t .To reduce the number of input parameters and thereby computation time and since our main focus remains on the neutrino sector, we avoid this scenario.

Constraints from direct searches of sparticles from the LHC:
The LHC collaboration has extensively searched for the supersymmetric partners of the SM particles (sparticles) from Run-I and Run-II data in various final states and without any statistically significant deviation of data over the SM prediction, the LHC has imposed stringent lower limits on the sparticle masses.For a summary of the ATLAS and CMS SUSY searches see Ref. [100,101].Here we briefly mention the most stringent bounds which are valid for simplified scenarios with specific assumptions on branching ratios and mostly for massless or relatively light neutralino.While doing so, we summarise the limits corresponding to both R-parity conserving and violating scenarios.Unless the R-parity violating couplings are large enough, in some cases the bounds corresponding to the R-parity conserving scenario may also be applicable.Eventually, the applicable limits depend on the mass and decay branching ratios of the relevant particle.We have applied the limits accordingly.
• In the simplified RPC-SUSY framework with different choices of decay mode and branching, the ATLAS and CMS collaborations have now pushed the lower limit of gluino mass to ∼ 2.0 -2.3 TeV for m χ 0 1 upto 600 GeV [102][103][104][105].The LHC has also pushed the light squarks mass to ∼ 1.85 TeV [102,105,106], the lightest stop to ∼ 1.0 -1.3 TeV [104,[107][108][109][110] for m χ 0 1 upto 300 GeV.In RPV SUSY scenarios if LSP decays to charged leptons via LLE couplings, then the limits obtained by the LHC are relatively stronger.For example, the ATLAS collaboration has excluded gluino mass up to 2.5 TeV6 in such scenarios [70].We have kept the squarks and gluino masses at 3 TeV to evade the current LHC constraints.
• The limits on electroweak (EW) sparticles i.e., sleptons and electroweakinos are relatively weaker compared to strong sparticles.For example, in RPC scenarios, the LHC collaboration has searched for electroweakinos for different decay modes like slepton mediated, WZ and Wh mediated final states and has excluded wino like χ± 1 upto ∼ (1.0 -1.4) TeV [106,[112][113][114].For slepton pair production, slepton mass upto ∼ 700 GeV is excluded for massless neutralino [112] for the universal slepton mass scenario.However, it should be noted that these strong limits are not always applicable to the overall parameter space of the realistic SUSY scenarios, e.g., compressed SUSY scenarios [115].The ATLAS and CMS collaborations have also interpreted the limits in RPV SUSY scenarios with LLE, U DD and bRPV couplings.For LLE type coupling, slepton and chargino (wino type) masses are excluded upto 1.2 and 1.6 TeV respectively [70].Limits in models with U DD coupling get drastically reduced as compared to RPC scenarios [71].Again for pure higgsino type χ± 1 χ0 2 pair production in RPC scenarios, the lower bound on m χ± 1 = χ0 2 is reduced to ∼ 210 GeV [72].On the other hand, higgsinos in bRPV scenarios are excluded upto 440 GeV [73].It is worth mentioning that a recent article [116] has provided an updated detailed summary of the possible gaps in RPV-MSSM searches at the LHC.In the Refs.[37,116], the authors have meticulously classified the various possible trilinear RPV-MSSM signatures at the LHC.In their analysis, they have studied both direct and indirect production of various LSPs and derive limits on SUSY masses, which are comparable or an improvement on those obtained in the R-parity conserving scenarios.However, these limits are not directly applicable in our bilinear RPV scenario.For more details, refer to [116].
• The most stringent limit on M A comes from the heavy Higgs searches in the H/A → τ + τ − decay channel and typically M A < 1.5 (1.0) TeV is excluded for tan β < 21 (8).[117,118].

Survey of parameter space
In the bRPV model, we have nine RPV parameters -three bRPV couplings (ϵ 1 , ϵ 2 , ϵ 3 ), three corresponding soft coupling parameters (B 1 , B 2 , B 3 ), and three sneutrino VEV parameters (v 1 , v 2 , v 3 ).The resulting light neutrino masses are quite sensitive to the choices of all these nine parameters.Apart from that, we also have different MSSM parameters which are essential to achieve our objective.Among them, µ and tan β are the most relevant ones.Given the large number of independent parameters in low scale MSSM, we fix some parameters which are not directly affecting the neutrino sector.For example, trilinear coupling for the third generation squark (A t ) is a very important parameter to achieve a 125 GeV Higgs in the model but it does not have a direct impact on the light neutrino masses and mixing angles.As large values of A t is required to obtain m h ∼ 125 GeV, we have chosen A t = -3.5 TeV.
We have also fixed all the three generation squarks soft masses (m q) and slepton soft masses (m l) at 3 and 2 TeV respectively to ensure that none of the current exclusion bounds on sparticles masses affect our parameter space.Similarly, we have fixed the bino (M 1 ), wino (M 2 ), gluino (M 3 ) soft masses and M A at 0.3, 1.2, 3.0 and 3.0 TeV respectively.
Following a literature survey [24,55] and some preliminary computation, we decided on exhaustive ranges for the bRPV model parameters ϵ i , v i and B i to ensure that the light neutrino mass square differences and mixing angles are generated in the correct order.Since our objective is to probe both the normal and inverted hierarchy scenarios, we do not presume any hierarchy in the choices of these parameters generation-wise and keep the ranges uniform over all three generations.We keep the choices conservative for the other two input parameters µ and tan β.The ranges of these input parameters are enlisted in Table 3.
Table 3: Ranges of eleven input parameters considered in our analysis.

Analysis set-up
Given the data and set of model parameters, as described earlier, we now proceed to calculate the posterior probability distribution in order to locate the favored parameter space.This is obtained using the MCMC technique which maximizes the likelihood function (or minimizes χ 2 ) defined as where L is the negative of the log-likelihood and calculated using where Γ obs i represents the set of n obs observed data points with corresponding errors σ i on them and Γ th i is the calculated value of each observable using our theoretical model.Altogether, we have total 15 independent observables (two neutrino masssquared differences and three mixing angles, SM-like Higgs mass and seven coupling modifiers, two flavor constraints from rare b-decay) and 11 free parameters (ϵ i , B i , v i , tan β and µ), so that the degrees of freedom (DoF) for the χ 2 distribution is 4.
The bRPV SUSY spectrum is generated by SPheno [95,96] which calculates the Higgs masses considering up to two-loop correction [119] and all the other particle masses at one-loop level.The bRPV model was implemented in SPheno using SARAH [97][98][99].For the MCMC-based likelihood analysis, we use publicly available code emcee [120] which is a Python implementation of the affine-invariant ensemble sampler.The code itself ensures an efficient exploration of parameter space even if there are strong degeneracies among those.We use a flat prior on all the parameters as mentioned in Table 3.To get the desirable acceptance fraction for the proposed MCMC steps, we use a relatively higher number of random walkers (500) each with a sufficient number of steps (chain length).An auto-correlation analysis has also been carried out which ensures that the convergence criterion for each chain is well-satisfied.

Results and discussion
In this section, we present our results for both the hierarchy scenarios -NH and IH.For each scenario, we present the marginalized distributions for different input parameters.We also mention the best-fit and mean values of the parameters along with the χ 2 min .Finally, we compare the allowed parameter space in NH and IH scenarios.

Normal Hierarchy scenario
In the NH scenario, the heaviest of the three light neutrinos is dominantly τ flavored.The second heaviest state is a nearly equal admixture of all three flavors (e, µ, τ ) whereas the lightest neutrino mass eigenstate is dominantly e flavored.This hierarchy is expected to be highlighted by the choices of the neutrino sector parameters in the most probable region.
The best-fit, mean values along with 95% C.L. of input parameters are listed in the Table 4.The observable values for the best-fit point are also mentioned in Table 4.The best-fit point for NH corresponds to χ 2 min = 3.46, degrees of freedom (DoF) = 4 and χ 2 min /DoF = 0.8657 The contribution to the χ 2 min from the neutrino observables for the best-fit point is ∼ 0.097.It is evident that the model can fit the neutrino oscillation data quite nicely.The neutrino masses at the above-mentioned best-fit point are: m ν 1 = 3.58 × 10 −6 eV, m ν 2 = 8.67 × 10 −3 eV and m ν 3 = 5.05 × 10 −2 eV and their sum is m ν i = 0.059 eV.It may be noted that the 2σ upper limit on the sum of neutrino masses, in NH scenario, coming from the cosmological data [8] is m ν i < 0.12 eV.
The marginalized posterior distribution for various input parameter values are presented in Fig. 3 and Fig. 5.The dark blue and the light blue regions in Fig. 3 subjected to the specific choices of A t = −3.5 TeV and M A = 3 TeV as mentioned before.Clearly, even the 2σ allowed region for tan β is highly constrained.Note that, this stringent constraint is entirely due to neutrino oscillation data which proves to be far stricter than the flavor and Higgs sector observables.tan β affects the tree and loop contributions to the light neutrino masses in a contrasting manner as discussed in Sec. 2 and is expected to be highly constrained given the small margins of uncertainty in the neutrino oscillation data.
Figure 4 represents how the two light neutrino mass squared differences vary with tan β when all other parameters are kept fixed at their best-fit values.It clearly shows that there is only one region of tan β where both ∆m 2  21 and ∆m 2 31 can be fit simultaneously within their 2σ allowed ranges.Understandably, the allowed tan β here is just a number (12.11, i.e., best-fit value), and if one varies the other parameters, the allowed range at 95% C.L. is obtained as shown in Figure 3.The allowed regions for tan β in the µ − tan β plane at 68% and 95% C.L. are ∼ [10.7-15.3]and [9.1-18.0]respectively.The choice of µ is also mostly restricted from the neutrino   ) and the horizontal magenta (cyan) shaded region shows the corresponding 2σ allowed region obtained from global analysis of neutrino oscillation data [8].The mass-squared differences are taken in eV 2 units.sector data.Similarly, the 68% and 95% C.L. allowed ranges of µ found out to be around [1075-1460] and [1000-1510] GeV respectively for the fixed set of other parameters mentioned in Sec.3.
The dark blue and light blue regions represent contours at 68% and 95% C.L. The dashed grey lines indicate the best-fit values for both parameters.
Fig. 5 shows the 2D marginalized posterior probability distributions in (a) along with the allowed regions at 68% C.L. and 95% C.L. The heaviest of the three neutrinos gets its mass at tree level with some additional correction from the loop diagrams.This neutrino mass is, therefore, mostly driven by the ϵ parameter.The heaviest of the three in this case has to be dominantly τ -flavored.As a result, the ϵ 3 and v 3 are expected to be the largest among the respective ϵ i and v i parameters.On the other hand, first generation parameters ϵ 1 and v 1 are expected to be the smallest since the heaviest neutrino has next to zero admixture of electron neutrino (see Fig. 5(a)-(d)).The second heaviest neutrino gets most of its contribution from the BB-loop and it has a comparable admixture of all three neutrino flavors.As a result, we can observe nice correlations between two B i 's as shown in Fig. 5(e) and Fig. 5(f).These contributions are already loop-suppressed.In addition to that the value of tan β is also restricted from tree-level calculation by the smallness of neutrino mass scale (see Sec.2).Hence for these contributions to neutrino masses to be significant, the B i parameters have to be much larger compared to ϵ i parameters as evident from Fig 5 .The ϵB-loop contributions are further suppressed due to their dependence on ϵ.As a result, B 1 is expected to be relatively larger than B 2 since the lightest state is dominantly electron neutrino-like.From the marginalized posterior distribution in ϵ 1 -ϵ 2 and ϵ 1 -ϵ 3 plane, it is evident that all the ϵ i are tightly constrained.We obtain the 2σ ranges of ϵ 1 , ϵ 2 , ϵ 3 as (-1.65 to 1.08)×10 −2 , (-2.30 to 0.53)×10 −2 and (-5.76 to -0.99)×10 −2 GeV respectively.Similarly, the 2σ allowed regions for v 1 , v 2 , v 3 from Fig. 5 (c) and (d) are (1.13-5.11)×10−4 GeV, (1.85-6.07)×10−4 GeV and (6.81-13.6)×10−4 GeV.The same regions for B 1 , B 2 and B 3 are (250-860) GeV, (0-575) GeV and (1380-2550) GeV respectively.
In Fig. 8 (Appendix-A), we also present the corner plot of all the input parameters, which shows all the 1D and 2D marginalized distribution of the posterior probability and reflects the covariances between parameters.As mentioned in Sec.1, to get the light mass eigenstates as the three neutrino masses from the 7 × 7 neutralino mass matrix, the alignment between ϵ i and v i parameters is required (i.e., v i ∝ ϵ i ).This is also reflected in the corner plot in v i -ϵ i plane (for example the first column-third row plot in Fig. 8 from the top represents the correlation in v 1 -ϵ 1 plane).We have also compared our results with the previous analysis in the bRPV scenarios.In most of the previous works [53,54,[56][57][58]60], the authors had analyzed the mSUGRA scenarios considering the then available neutrino data and these works were published before the Higgs discovery era.Among these analyses, relevant parameter space scanning was presented in Ref. [53,60] with contemporary neutrino data and our best-fit point along with 2σ allowed regions have significant overlap with the previous results.While deriving these results, all other SUSY parameters are kept fixed at values mentioned at the beginning of Sec.We perform a similar analysis assuming that the light neutrino masses obey an inverted hierarchy.The choices of input parameters and ranges are the same as the previous analysis and details are given in Sec.3.3.In Table 5, we present the best-fit and mean values along with the 95% C.L. allowed regions of individual input parameters.The observable values for the best-fit point are also listed in the last two columns of Table 5.The best-fit point for IH corresponds to χ 2 min = 3.38, DoF = 4 and χ 2 min /DoF = 0.845.The contribution to the χ 2 min from the neutrino observables for the best-fit point is small, 0.133.The neutrino masses at the above mentioned best-fit point are m ν 1 = 4.96 × 10 −2 eV, m ν 2 = 5.03 × 10 −2 eV and m ν 3 = 1.23 × 10 −5 eV and their sum is 3 i=1 m ν i ≈ 0.1 eV.This evades the 2σ upper limit on the sum of neutrino masses coming from the cosmological data i.e., 3  i=1 m ν i < 0.15 eV [8].The resultant contour plot in the µ − tan β plane is shown in Fig. 6 and the dark (light) purple regions correspond to 68% (95%) confidence contours.The 2σ allowed regions for µ and tan β are [1270-1635] GeV and [6.2-10.2]respectively.The heaviest neutrino state in this scenario contains an almost equal admixture of all the three neutrino flavors.Tree level contribution from any single flavor cannot be very large and as a result the mixing of any single flavor neutrino with neutralino states typically would have to be smaller compared to that in a normal hierarchy scenario.That explains the comparatively smaller tan β and heavier µ values obtained in this case following Eq.2.4 and Eq.2.5.
The marginalized posterior distributions of the RPV parameters are shown in Fig. 7 with the same color convention as the previous figure.In the IH scenario, the lightest state is dominantly τ -flavored with substantial admixture from µ flavor as well.The hierarchy of the other two mass eigenstates remains the same as in normal hierarchy.However, their physical masses are now larger compared to the NH case.This leads to substantial changes in the mixing pattern of the light neutrinos.Take the heaviest state as an example.Unlike in the NH scenario, this state has an almost equal admixture of all three neutrino flavors.In addition to that, the mass of this state is larger than the corresponding state in NH.This explains the slightly different ranges of the ϵ i and the v i parameters.The hierarchy in values depends upon the relative contributions required in the admixture of the states.The ranges of B 1 and B 2 are quite different from the NH scenario (see Fig. 5(e)-(f) and Fig. 7(e)-(f)).Again, this has to do with the different amount of loop correction required from different generations depending on the changed hierarchy.Take the state ν 1 for example.In both scenarios, the mass of this state is generated at one loop, but in the IH scenario,   the physical mass of this state is larger than that in the NH scenario.This requires larger loop contributions, which explain the larger values required for B 1 and B 2 .The 2D marginalized distributions suggest that the allowed regions are more constrained in the IH scenario compared to the NH scenario.We present the corner plot for all the input parameters i.e., all the 1D and 2D marginalized distribution of the posterior probability in Fig. 9 of Appendix-B.

Addressing Anomalous Muon magnetic moment
We now proceed with our allowed parameter region to explore the possibility of explaining the existing 4.2σ excess (see Eq.1.2) in the measurement of muon (g -2) [67,68].The excess contribution arising from both RPC and RPV processes is denoted as ∆a µ .The RPC SUSY parameter space has been studied widely8 [69,[124][125][126][127][128][129][130][131][132][133][134][135][136][137][138][139].We aim to explore whether the additional contributions arising in presence of RPV couplings can make a difference.As mentioned earlier, the new loop contributions that we obtain in the present scenario are through the mixing of sneutrino-higgs, neutralino-neutrino and chargino-charged lepton states.However, given the allowed regions of parameter space obtained in our study, it is evident that these new contributions will always be subleading because they depend on the ϵ i and v i parameters which have to be quite small to address experimental observations.
One would have ideally added ∆a µ as another observable in the analysis itself to get a complete picture, but obtaining the required excess in ∆a µ also requires us to vary bino, wino as well as the slepton soft mass parameters.That would have made our analysis more complex and time-consuming.Hence we have kept our ∆a µ analysis separate and present our results in terms of some chosen benchmark points.In Table 6, we have enlisted all the input parameters and relevant observables along with the obtained value of ∆a µ .We have tried to find benchmark points that can explain ∆a µ within the 2σ allowed range of the global average mentioned in Eq. 1.2.We have considered the NH scenario here for this analysis.One can obtain similar results for the IH scenario as well.For the benchmark points (BP-I and BP-II), we have fixed the following parameters as: A t = -3.5 TeV, M A = 3 TeV, M 3 = 3 TeV, m q = 3 TeV (all 3 generations) and m ẽL/R = m τL/R = 2 TeV to evade the current LHC bounds on the sparticle masses as mentioned in section 3.3.The relevant smuon mass m μ is mentioned in the Table 6 for the two benchmark points.BP-I (BP-II) satisfies the observed value of ∆a µ within 1σ (2σ) limit.This difference happens due to the larger smuon mass in BP2.Note that, the new physics contributions obtained for Table 6: Details of benchmark points (BP) satisfying the muon (g -2) data along with other observables.
these spectra are very similar to R-parity conserving scenario.Except for the LSP, the decay of all other SUSY particles is dominantly R-parity conserving because of the smallness of the RPV parameters.Hence the applicable collider limits for these sparticles happen to be the same as in R-parity conserving MSSM9 except for the LSP.The LSP decays through RPV couplings and for its mass existing RPV limits have been taken into account [140].

Conclusion
The existence of non-zero masses and non-trivial mixing among the light neutrino states have been established beyond any doubt by the neutrino oscillation experiments.The standard model cannot address this phenomenon within its framework and that makes this observation one of the most robust indications for new physics beyond the standard model.Any completely new physics model, therefore, should be able to address this issue.Theoretically, Supersymmetry remains one of the most well-motivated new physics scenarios till date and hence various supersymmetric scenarios have been extensively studied both by theoretical and experimental collaborations.Supersymmetric extension of the standard model with conserved R-parity requires an additional mechanism to explain the neutrino oscillation data while with bilinear R-parity violating scenario, one can easily address this long-standing issue without extending the model any further.In the absence of any robust indication of new physics at the LHC, the need of the hour is to study the existing models under the light of the plethora of experimental data at our disposal and try to constrain the new physics parameter regions as much as possible.Neutrino oscillation data is more precise and capable of constraining the new physics parameter space compared to direct search data, as shown in this paper.In addition to that, taking into account the precision Higgs data and flavor data, we have computed the most restricted bRPV SUSY parameter space that can be obtained (with χ 2 min / DoF ∼ 1 for both the neutrino hierarchical scenarios) using MCMC analysis.Moreover, our analysis also considers the existing direct search limits.
Our results show that the ϵ i and the v i parameters in particular have to be quite small which means that except for the LSP, the decays of the neutralinos and charginos are expected to be mostly similar to the R-parity conserving scenario.The B i parameters can be comparatively much larger since they only contribute to the neutrino masses and mixing angles at the loop level.Since these parameters are responsible for slepton-charged Higgs and sneutrino-neutral Higgs mixing, the phenomenology of these particles is expected to differ accordingly from the R-parity conserving scenario.Thus, the results derived here will be extremely helpful in making any future studies more focused and results more predictive.Understandably, a slight change in one of the experimental data points can change the best-fit point quoted in this paper, but unless the experimental results change too much, the 2σ allowed regions are expected to remain similar.We further proceed to address the existing anomalous muon (g -2) result with the constrained parameter space obtained in this work.We find out that unless we are in the compressed region where the slepton and gaugino masses are lying close to each other, it is very difficult to abide by the collider bounds and still explain the muon (g -2) excess.In this regard, our results are quite similar to what is expected in the R-parity conserving scenario since the smallness of the RPV couplings ensures that there are no new significant contributions arising because of R-parity violation.

Figure 1 :
Figure 1: Depiction of generation of light neutrino masses at the tree level.The cross represents mass insertion.

Figure 2 :
Figure 2: Contribution to neutrino masses from one loop diagrams via -(a) BB loop and (b) ϵB loop.For the ϵB loop diagram there will be another diagram with i ↔ j.The cross represents mass insertion.

Figure 3 :
Figure 3: Marginalized posterior distribution with 68% (dark blue) and 95% (light blue) C.L. contours in the µ − tan β plane for NH scenario.The dashed grey lines indicate the best-fit values.

Figure 4 :
Figure 4: The mass square differences (∆m 2 21 and ∆m 2 31 ) vs tan β plot where all other parameters are kept fixed at the best-fit point.The line with magenta (cyan) color corresponds to ∆m 2 21 (∆m 2 31) and the horizontal magenta (cyan) shaded region shows the corresponding 2σ allowed region obtained from global analysis of neutrino oscillation data[8].The mass-squared differences are taken in eV 2 units.
which looks as shown in equation 3.1.12c 23 − c 12 s 13 s 23 c 12 c 23 − s 12 s 13 s 23 c 13 s 23 s 12 s 23 − c 12 s 13 c 23 −c 12 s 23 − s 12 s 13 c 23 c 13 c 23Here, c ij and s ij represent cos θ ij and sin θ ij (i, j are generation indices) respectively.For simplicity, we have kept the CP-violating phase as zero and work with the five aforementioned observables, namely, ∆m 2 21 , |∆m 2 31 |, θ 12 , θ 13 and θ 23 .Several groups have done global fits with neutrino oscillation data obtained from different experiments

Table 1 :
[8]e have considered both normal and inverted hierarchy List of neutrino sector observables provided by the updated global fit analysis[8]of neutrino oscillation data obtained from different neutrino experiments.Here NH and IH refers to Normal Hierarchy and Inverted Hierarchy scenario.
scenarios separately in our analysis for comparative study.The choice of hierarchy is expected to be reflected in the resulting parameter space.Note that, the neutrino oscillation experiments also measure the CP-violating phase (δ CP ).The updated global

Table 2 :
Higgs boson coupling strength modifiers obtained by the CMS collaboration using LHC Run-II 137 f b −1 data [61].

Table 4 :
, represent the 68% and 95% confidence contours in µ − tan β plane.This result is Best-fit, mean values along with the 95%C.L. of all the free parameters and best-fit values for the observables in NH scenario are shown here.The last row represents χ 2 min and χ 2 min /DoF.

Table 5 :
3.3.Varying the SUSY parameters further can alter the results.For illustration, we have presented the results with different choices of M A and A t in Appendix.C. Best-fit, mean values along with the 95%C.L. of all the free parameters and best-fit values for the observables in the NH scenario are shown here.The last row represents χ 2 min and χ 2 min /DoF .
Marginalized posterior distribution with 68% (dark purple) and 95% (light purple) C.L. contours in the µ − tan β plane for IH scenario.The dashed grey lines indicate the best-fit values.