1 Introduction

The measurements of the properties and signal rates of the discovered scalar boson at the LHC [1, 2], indicate that it is compatible with the expectations of the Standard Model (SM). Besides, no hints for new physics have been detected yet despite of numerous searches and tremendous efforts of the experimental collaborations. As a consequence, extensions of the SM such as low-energy Supersymmetry (SUSY) are being severely constrained, namely the parameter space of SUSY models is shrinking considerably. This renders the detailed analyses of Higgs properties, signal rates and couplings to SM particles as important as the search for new particles.

Concerning the latter, the search of SUSY particles has been focused mainly on signals with missing transverse energy inspired by R-parity conserving (RPC) models, such as the Minimal Supersymmetric Standard Model (MSSM) [3,4,5,6,7]. There, significant bounds on sparticle masses have been obtained [8], especially for strongly interacting sparticles whose masses must be above about 1 TeV. Less stringent bounds of about 100 GeV have been obtained for weakly interacting sparticles, with the bino-like neutralino basically not constrained due to its small pair production cross section. Qualitatively similar results have also been obtained in the analysis of simplified R-parity violating (RPV) scenarios with trilinear lepton- or baryon-number violating terms [9], assuming a single channel available for the decay of the LSP into leptons. However, this assumption is not possible in other RPV scenarios, such as the ‘\(\mu \) from \(\nu \)’ supersymmetric standard model (\(\mu \nu \mathrm{SSM}\)) [10], where the several decay branching ratios (BRs) of the lightest supersymmetric particle (LSP) significantly decrease the signal. This implies that the extrapolation of the usual bounds on sparticle masses to the \(\mu \nu \mathrm{SSM}\) is not applicable. For example, it was shown in Refs. [11, 12] that the LEP lower bound on masses of slepton LSPs of about 90 GeV obtained in trilinear RPV [13,14,15,16,17,18] is not applicable in the \(\mu \nu \mathrm{SSM}\). For the bino LSP,Footnote 1 only a small region of the parameter space of the \(\mu \nu \mathrm{SSM}\) was excluded [23] when the left sneutrino is the next-to-LSP (NLSP) and hence a suitable source of binos. In particular, the region of bino (sneutrino) masses 110–150 (110–160) GeV.

Concerning Higgs physics, the \(\mu \nu \mathrm{SSM}\) expands the singlet superfield of the NMSSM [24, 25] to three right-handed neutrino superfields [10, 26]. In the NMSSM, various works, using different methods (see e.g. [27,28,29,30,31,32,33]), have been dedicated to sample the parameter space in the light of a given set experimental data, and vast regions have been explored. In this work, we use a powerful likelihood data-driven method based on the algorithm called MultiNest [34,35,36] for sampling the Higgs sector of the \(\mu \nu \mathrm{SSM}\). Since three families of right-handed neutrino superfields are present in the model in order to solve the \(\mu \) problem and simultaneously reproduce neutrino physics, the new couplings and sneutrino vacuum expectation values (VEVs) produce a substantial mixing among the three right sneutrinos and the doublet-like Higgses. Although a detailed analysis of this sector was performed in Ref. [26], finding viable regions that avoid false minima and tachyons, as well as fulfill the Landau pole constraint, it was carried out prior the discovery of the SM-like Higgs boson, and therefore the issue of reproducing Higgs data was missing. In Ref. [22], this issue was taking into account to perform an analytical estimate of all the new two-body decays for the SM-like Higgs in the presence of light scalars, pseudoscalars and neutralinos. More recently, in Refs. [37, 38], in addition to performing the complete one-loop renormalization of the neutral scalar sector of the \(\mu \nu \mathrm{SSM}\), interesting benchmark points (BPs) with singlet-like eigenstates lighter than the SM-like Higgs boson were studied.

Given the increasing data including the properties of the SM-like Higgs and the exclusion limits on scalars from extended Higgs-sneutrino sectors provided by the combined 7-, 8- and 13-TeV searches at the LHC, and also by other results such as flavor observables, it appears relevant to re-investigate the \(\mu \nu \mathrm{SSM}\) parameter space to simultaneously accommodate this new scalar and its properties, the exclusion limits and to explore the phenomenological consequences respecting various experimental results. To carry this out, the likelihood data-driven method used in our analysis presents advantages over traditional ones such as those based on random grid scans or Chi-square methods, since it is much more efficient in the computational effort required to explore a parameter space. Also, since it uses a Bayesian approach, it allows to take easily into account all relevant sources of uncertainties in the likelihood. In addition, given the accumulation of data from various experimental collaborations, this method provides a convenient approach to qualitatively explore beyond standard models compared to simplified methods.

The paper is organized as follows. In Sect. 2, we will briefly review the \(\mu \nu \mathrm{SSM}\) and its relevant parameters for our analysis of the Higgs sector. This sector will be studied in detail in Sect. 3, where the mixing among doublet-like Higgses and right and left sneutrinos will be explained. We will pay special attention to accommodate the correct mass of the SM-like Higgs, depending on the values of the couplings \(\lambda \) among right sneutrinos and doublet-like Higgses, and the masses of the singlet-like eigenestates. Subsequently, in Sect. 4 we will discuss the strategy that we will employ to perform scans searching for points of the parameter space of our scenario compatible with current experimental data on Higgs physics, as well as flavor observables. The results of these scans will be presented in Sect. 5, and applied to show that there are large viable regions of the parameter space of the \(\mu \nu \mathrm{SSM}\). Our conclusions and prospects for future work are left for Sect. 6. Finally, useful formulae, figures, and BPs are given in the Appendices. In Appendix A, the Higgs-right sneutrino mass submatrices are written. In Appendix B, results from the \(\lambda -\kappa \) plane are shown for different values of the other parameters, using several figures for each scan performed. In Appendix C, several BPs showing interesting characteristics of the model are given.

2 The \(\mu \nu \mathrm{SSM}\)

The \(\mu \nu \mathrm{SSM}\) [10, 26] is a natural extension of the MSSM where the \(\mu \) problem is solved and, simultaneously, the neutrino data can be reproduced [10, 19, 20, 26, 39, 40]. This is obtained through the presence of trilinear terms in the superpotential involving right-handed neutrino superfields \({\hat{\nu }}^c_i\), which relate the origin of the \(\mu \)-term to the origin of neutrino masses and mixing angles. The simplest superpotential of the \(\mu \nu \mathrm{SSM}\) [10, 26, 41] with three right-handed neutrinos is the following:

$$\begin{aligned} W= & {} \sum _{a,b}\sum _{i,j} \epsilon _{ab} \left( Y_{e_{ij}} {\hat{H}}_d^a {\hat{L}}^b_i {\hat{e}}_j^c + Y_{d_{ij}} {\hat{H}}_d^a {\hat{Q}}^{b}_{i} {\hat{d}}_{j}^{c} + Y_{u_{ij}} {\hat{H}}_u^b {\hat{Q}}^{a}_i {\hat{u}}_{j}^{c} \right) \nonumber \\&\quad + \sum _{a,b}\sum _{i,j}\epsilon _{ab} \left( Y_{{\nu }_{ij}} {\hat{H}}_u^b {\hat{L}}^a_i {\hat{\nu }}^c_j - \lambda _i {\hat{\nu }}^c_i {\hat{H}}_u^b {\hat{H}}_d^a \right) \nonumber \\&\quad + \sum _{i,j,k} \frac{1}{3} \kappa _{ijk} {\hat{\nu }}^c_i{\hat{\nu }}^c_j{\hat{\nu }}^c_k, \end{aligned}$$
(1)

with \(a,b=1,2\) \(SU(2)_L\) indices with \(\epsilon _{ab}\) the totally antisymmetric tensor \(\epsilon _{12}=1\), and \(i,j,k=1,2,3\) the usual family indices of the SM.

The simultaneous presence of the last three terms in Eq. (1) makes it impossible to assign R-parity charges consistently to the right-handed neutrinos (\(\nu _{iR}\)), thus producing explicit RPV (harmless for proton decay). Note nevertheless, that in the limit of neutrino Yukawa couplings \(Y_{{\nu }_{ij}} \rightarrow 0\), \({\hat{\nu }}^c_i\) can be identified in the superpotential as pure singlet superfields without lepton number, similar to the singlet of the NMSSM, and therefore R parity is restored. Thus, \(Y_{\nu }\) are the parameters which control the amount of RPV in the \(\mu \nu \mathrm{SSM}\), and as a consequence this violation is small since the size of \(Y_{\nu }\lesssim 10^{-6}\) is determined by the electroweak-scale seesaw of the \(\mu \nu \mathrm{SSM}\) [10, 26].

The tree-level neutral scalar potential \(V=V_{\text {soft}} + V_F + V_D\), receives in addition to the usual D- and F-term contributions that can be found e.g. in Refs. [26, 41], the following contribution from the soft SUSY-breaking Lagrangian:

$$\begin{aligned} V_{\text {soft}}= & {} \sum _{i,j,k}\left( T_{{\nu }_{ij}} H_u^0 {\widetilde{\nu }}_{iL} {\widetilde{\nu }}_{jR}^* - T_{{\lambda }_{i} } {\widetilde{\nu }}_{iR}^* H_d^0 H_u^0 \right. \nonumber \\&\quad \left. + \frac{1}{3} T_{{\kappa }_{ijk}} {\widetilde{\nu }}_{iR}^* {\widetilde{\nu }}_{jR}^* {\widetilde{\nu }}_{kR}^*\ + \text {h.c.} \right) \nonumber \\&\quad + \sum _{i,j} \left( m_{\widetilde{L}_{ij}}^2\widetilde{\nu }_{iL}^* {\widetilde{\nu }}_{jL} + m_{\widetilde{\nu }_{ij}}^2 \widetilde{\nu }_{iR}^* {\widetilde{\nu }}_{jR} \right) \nonumber \\&\quad + m_{H_d}^2 {H^0_d}^* H^0_d + m_{H_u}^2 {H^0_u}^* H^0_u. \end{aligned}$$
(2)

If we follow the assumption based on the breaking of supergravity that all the trilinear parameters are proportional to their corresponding couplings in the superpotential [42], we can write

$$\begin{aligned} T_{{\nu }_{ij}} = A_{{\nu }_{ij}} Y_{{\nu }_{ij}}, \;\; T_{{\lambda }_i}= A_{{\lambda }_i} \lambda _i, \;\; T_{{\kappa }_{ijk}}= A_{{\kappa }_{ijk}} \kappa _{ikj}, \end{aligned}$$
(3)

and the parameters A substitute the T as the most representative. We will use both type of parameters in our discussions. It is worth noticing here that we do not use the summation convention on repeated indices throughout this work.

The soft terms of Eq. (2) induce the electroweak symmetry breaking (EWSB) in the \(\mu \nu \mathrm{SSM}\). The minimization equations, with the choice of CP conservation,Footnote 2 can also be found in Refs. [26, 41]. For neutral Higgses (\(H^0_{u,d}\)) and right (\({\widetilde{\nu }}_{iR}\)) and left (\({\widetilde{\nu }}_{iL}\)) sneutrinos defined as

$$\begin{aligned} H_d^0= & {} \frac{1}{\sqrt{2}} \left( H_{d}^{\mathcal {R}} + v_d + i\ H_{d}^{\mathcal {I}}\right) ,\nonumber \\ H^0_u= & {} \frac{1}{\sqrt{2}} \left( H_{u}^{\mathcal {R}} + v_u +i\ H_{u}^{\mathcal {I}}\right) , \end{aligned}$$
(4)
$$\begin{aligned} \widetilde{\nu }_{iR}= & {} \frac{1}{\sqrt{2}} \left( \widetilde{\nu }^{{\mathcal {R}}}_{iR}+ v_{iR} + i\ \widetilde{\nu }^{{\mathcal {I}}}_{iR}\right) , \nonumber \\ \widetilde{\nu }_{iL}= & {} \frac{1}{\sqrt{2}} \left( \widetilde{\nu }_{iL}^{\mathcal {R}} + v_{iL} +i\ \widetilde{\nu }_{iL}^{\mathcal {I}}\right) , \end{aligned}$$
(5)

the following vacuum expectation values (VEVs) are developed:

$$\begin{aligned} \langle H_{d}^0\rangle= & {} \frac{v_{d}}{\sqrt{2}},\quad \langle H_{u}^0\rangle = \frac{v_{u}}{\sqrt{2}},\nonumber \\ \langle \widetilde{\nu }_{iR}\rangle= & {} \frac{v_{iR}}{\sqrt{2}},\quad \langle \widetilde{\nu }_{iL}\rangle = \frac{v_{iL}}{\sqrt{2}}, \end{aligned}$$
(6)

with \(v_{iR}\sim \) TeV whereas \(v_{iL}\sim 10^{-4}\) GeV. The latter result is because of the contributions proportional to \(Y_{\nu }\) to the \(v_{iL}\) minimization equations. They enter through \(V_F\) and \(V_{\text {soft}}\) (assuming \(T_\nu \) as in Eq. 3), and are small due to the electroweak-scale seesaw mentioned before that determines \(Y_{\nu }\lesssim 10^{-6}\). Note in this respect that the last term in the superpotential (1) generates dynamically Majorana masses for the right-handed neutrinos \(\sim \) TeV:

$$\begin{aligned} {{\mathcal {M}}}_{ij} =\sum _k {2}\kappa _{ijk} \frac{v_{kR}}{\sqrt{2}}. \end{aligned}$$
(7)

On the other hand, the fifth term generates an effective \(\mu \)-term \(\sim \) TeV:

$$\begin{aligned} \mu =\sum _i \lambda _i \frac{v_{iR}}{\sqrt{2}}. \end{aligned}$$
(8)

Given the structure of the scalar potential, the free parameters in the neutral scalar sector of the \(\mu \nu \mathrm{SSM}\) at the low scale \(M_{EWSB}= \sqrt{m_{{\tilde{t}}_l} m_{{\tilde{t}}_h}}\), are therefore: \(\lambda _i\), \(\kappa _{ijk}\), \(Y_{{\nu }_{ij}}\), \(m_{H_{d}}^{2}\), \(m_{H_{u}}^{2}\), \(m_{\widetilde{\nu }_{ij}}^2\), \(m_{\widetilde{L}_{ij}}^2\), \(T_{{\lambda }_i}\), \(T_{{\kappa }_{ijk}}\) and \(T_{{\nu }_{ij}}\). Using diagonal sfermion mass matrices, in order to avoid the strong upper bounds upon the intergenerational scalar mixing (see e.g. Ref. [43]), from the eight minimization conditions with respect to \(v_d\), \(v_u\), \(v_{iR}\) and \(v_{iL}\) one can eliminate the above soft masses in favor of the VEVs. In addition, using \(\tan \beta \equiv {v_u}/{v_d}\) and the SM Higgs VEV, \(v^2 = v_d^2 + v_u^2 + \sum _i v^2_{iL}={4 m_Z^2}/{(g^2 + g'^2)}\approx \) (246 GeV)\(^2\) with the electroweak gauge couplings estimated at the \(m_Z\) scale by \(e=g\sin \theta _W=g'\cos \theta _W\), one can determine the SUSY Higgs VEVs, \(v_d\) and \(v_u\). Since \(v_{iL} \ll v_d, v_u\), one has \(v_d\approx v/\sqrt{\tan ^2\beta +1}\). Besides, we can use diagonal neutrino Yukawa couplings, since data on neutrino physics can easily be reproduced at tree level in the \(\mu \nu \mathrm{SSM}\) with such structure, as we will discuss below. Finally, assuming for simplicity that the off-diagonal elements of \(\kappa _{ijk}\) and soft trilinear parameters T vanish, we are left with the following set of variables as independent parameters in the neutral scalar sector:

$$\begin{aligned} \lambda _i, \kappa _{i}, Y_{{\nu }_{i}}, \tan \beta , v_{iL}, v_{iR}, T_{{\lambda }_i}, T_{{\kappa }_{i}}, T_{{\nu }_i}, \end{aligned}$$
(9)

where \(\kappa _i\equiv \kappa _{iii}\), \(Y_{{\nu }_i}\equiv Y_{{\nu }_{ii}}\), \(T_{{\nu }_i}\equiv T_{{\nu }_{ii}}\) and \(T_{{\kappa }_i}\equiv T_{{\kappa }_{iii}}\). Note that now the Majorana mass matrix is diagonal, with the non-vanishing entries given by

$$\begin{aligned} {{\mathcal {M}}}_{i} ={2}\kappa _{i} \frac{v_{iR}}{\sqrt{2}}. \end{aligned}$$
(10)

The rest of (soft) parameters of the model, namely the following gaugino masses, scalar masses, and trilinear parameters:

$$\begin{aligned} M_1, M_2, M_3, m_{{\tilde{Q}}_{iL}}, m_{{\tilde{u}}_{iR}}, m_{{\tilde{d}}_{iR}}, m_{{\tilde{e}}_{iR}}, T_{u_{i}}, T_{d_{i}}, T_{e_{i}}, \end{aligned}$$
(11)

are also taken as free parameters and specified at low scale.

A further sensible simplification that we will also use in the next sections when necessary, is to assume universality of the parameters in Eq. (9) with the exception of those connected directly with neutrino physics such as \(Y_{{\nu }_i}\) and \(v_{iL}\), that must be non-universal to generate correct neutrino masses and mixing angles. Neither we will impose universality for \(T_{{\nu }_i}\), since they are connected with sneutrino physics as we will discuss in the next section, and a hierarchy of masses in that sector can be phenomenologically interesting [12]. We are then left with the following set of low-energy free parameters:

$$\begin{aligned} \lambda , \kappa , Y_{{\nu }_{i}}, \tan \beta , v_{iL}, v_{R}, T_{\lambda }, T_{\kappa }, T_{\nu _i}, \end{aligned}$$
(12)

where \(\lambda \equiv \lambda _i\), \(\kappa \equiv \kappa _i\), \(v_R\equiv v_{iR}\), \(T_{\lambda }\equiv T_{{\lambda }_i}\) and \(T_\kappa \equiv T_{{\kappa }_i}\). In this case, the three non-vanishing Majorana masses are equal \({{\mathcal {M}}}_i={{\mathcal {M}}}\), with

$$\begin{aligned} {{\mathcal {M}}} ={2}\kappa \frac{v_{R}}{\sqrt{2}}, \end{aligned}$$
(13)

and the \(\mu \)-term is given by

$$\begin{aligned} \mu =3\lambda \frac{v_{R}}{\sqrt{2}}. \end{aligned}$$
(14)

The new couplings and sneutrino VEVs in the \(\mu \nu \mathrm{SSM}\) induce new mixing of states. The associated mass matrices were studied in detail in Refs. [20, 26, 41]. Summarizing, there are eight neutral scalars and pseudoscalars (Higgses–sneutrinos), where after rotating away the pseudoscalar would be Goldstone boson we are left with seven pseudoscalar states. There are also eight charged scalars (charged Higgses–sleptons), five charged fermions (charged leptons–charginos), and ten neutral fermions (neutrinos–neutralinos).

Since reproducing neutrino data is an important asset of the \(\mu \nu \mathrm{SSM}\), in the following we will briefly review this issue. The neutral fermions have the flavor composition \((\nu _{iL},\widetilde{B}^0,{\widetilde{W}}^0,{\widetilde{H}}_{d}^0,{\widetilde{H}}_{u}^0,\nu _{iR})\). Thus, with the low-energy bino and wino soft masses, \(M_1\) and \(M_2\), of the order of the TeV, and similar values for \(\mu \) and \({\mathcal {M}}\) as discussed above, this generalized seesaw produces three light neutral fermions dominated by the left-handed neutrino (\(\nu _{iL}\)) flavor composition. In fact, as mentioned before, data on neutrino physics [44,45,46,47] can easily be reproduced at tree level [10, 19, 20, 26, 39, 40], even with diagonal Yukawa couplings [19, 39], i.e. \(Y_{{\nu }_{ii}}=Y_{{\nu }_{i}}\) and vanishing otherwise as used already in Eq. (9). A simplified formula for the effective mixing mass matrix of the light neutrinos is [39]:

$$\begin{aligned} (m_{\nu })_{ij}= & {} \frac{m_{{{\mathcal {D}}_i}} m_{{{\mathcal {D}}_j}} }{3{{\mathcal {M}}}} (1-3 \delta _{ij}) -\frac{(v_{iL}/\sqrt{2})(v_{jL}/\sqrt{2})}{2M^{\text {eff}}} \nonumber \\&\quad - \frac{ m_{{{\mathcal {D}}_i}}m_{{{\mathcal {D}}_j}} }{2M^{\text {eff}}} \frac{1}{3\lambda \tan \beta } \nonumber \\&\quad \times \left( \frac{v_{iL}/\sqrt{2}}{m_{{{\mathcal {D}}_i}}} + \frac{v_{jL}/\sqrt{2}}{m_{{{\mathcal {D}}_j}}} + \frac{1}{3\lambda \tan \beta } \right) , \end{aligned}$$
(15)

where we have defined the Dirac mass for neutrinos as

$$\begin{aligned} m_{{{\mathcal {D}}_i}}\equiv Y_{{\nu }_{i}} \frac{v_u}{\sqrt{2}}, \end{aligned}$$
(16)

and

$$\begin{aligned} M^{\text {eff}}\equiv & {} M - \frac{\left( v/\sqrt{2}\right) ^2 }{2\mu \left( {{\mathcal {M}}} \frac{v_R}{\sqrt{2}}+ 2\lambda \left( \frac{v}{\sqrt{2}}\right) ^2 \frac{\tan \beta }{1+\tan ^2\beta }\right) } \nonumber \\&\times \left[ 2{{\mathcal {M}}} \frac{v_R}{\sqrt{2}} \frac{\tan \beta }{1+\tan ^2\beta } +\lambda \left( \frac{{v}}{\sqrt{2}}\right) ^2\right] , \end{aligned}$$
(17)

with

$$\begin{aligned} \frac{1}{M} = \frac{g'^2}{M_1} + \frac{g^2}{M_2}. \end{aligned}$$
(18)

Here we have assumed universal \(\lambda _i = \lambda \), \(v_{iR}= v_{R}\), and \(\kappa _{i}=\kappa \) as in Eq. (12). The first term of Eq. (20) is generated through the mixing of \(\nu _{iL}\) with \(\nu _{iR}\)-Higgsinos, and the other two also include the mixing with gauginos. These are the so-called \(\nu _{R}\)-Higgsino seesaw and gaugino seesaw, respectively [39].

We are then left in general with the following subset of variables of Eqs. (9) and (11) as independent parameters in the neutrino sector:

$$\begin{aligned} \lambda _i, \kappa _i, Y_{\nu _i}, \tan \beta , v_{iL}, v_{iR}, M. \end{aligned}$$
(19)

In the numerical analyses of the next sections, it will be enough for our purposes to consider the sign convention where all these parameters are positive.

Under several assumptions, the formula for \((m_{\nu })_{ij}\) can be further simplified. Notice first that the third term is inversely proportional to \(\tan \beta \), and therefore negligible in the limit of large or even moderate \(\tan \beta \) provided that \(\lambda \) is not too small. Besides, the first piece inside the brackets in the second term of Eq. (17) is also negligible in this limit, and for typical values of the parameters involved in the seesaw also the second piece, thus \(M^{\text {eff}}\sim M\). Under these assumptions, the second term for \((m_{\nu })_{ij}\) is generated only through the mixing of left-handed neutrinos with gauginos. Therefore, we arrive to a very simple formula where only the first two terms survive with \(M^{\text {eff}}= M\) in Eq. (15):

$$\begin{aligned} (m_{\nu })_{ij} = \frac{m_{{{\mathcal {D}}_i}} m_{{{\mathcal {D}}_j}} }{3{{\mathcal {M}}}} (1-3 \delta _{ij}) -\frac{(v_{iL}/\sqrt{2})(v_{jL}/\sqrt{2})}{2M}. \end{aligned}$$
(20)

This expression can be used to understand easily the seesaw mechanism in the \(\mu \nu \mathrm{SSM}\) in a qualitative way. From this discussion, it is clear that \(Y_{\nu _i}\), \(v_{iL}\) and M are crucial parameters to determine the neutrino physics.

Let us finally point out that the accommodation of the SM-like Higgs boson discovered at the LHC is mandatory for any SUSY model. In the next section, we will review this subject in some detail in the context of the \(\mu \nu \mathrm{SSM}\). In order to carry it out, we will study the Higgs sector of the model and its enhanced particle content. Although the parameter space of the model given by Eqs. (9) and (11) is large, we will see in the next sections that some of these parameters are not relevant for the study of the Higgs sector, as also occurs for the neutrino sector, and therefore its analysis is simplified.

3 The Higgs sector of the \(\mu \nu \mathrm{SSM}\)

In the \(\mu \nu \mathrm{SSM}\), doublet-like Higgses are mixed with left and right sneutrinos, giving rise to \(8\times 8\) (‘Higgs’) mass matrices for scalar and pseudoscalar states. However, the \(5\times 5\) Higgs-right sneutrino submatrix is almost decoupled from the \(3\times 3\) left sneutrino submatrix due to the small values of \(Y_{{\nu }_{ij}}\) and \(v_{iL}\) in the off-diagonal entries [19, 26]. Thus, to accommodate the SM-like Higgs in the \(\mu \nu \mathrm{SSM}\), we can focus on the analysis of the Higgs-right sneutrino mass submatrix. The tree-leel entries of the scalar and pseudoscalar mass submatrices are shown in Appendix A using the parameters of Eq. (9). Upong diagonalization of the scalar submatrix, one obtains the SM-like Higgs, the heavy doublet-like neutral Higgs, and three singlet-like states. Similarly, upon diagonalization of the pseudoscalar submatrix, and after rotating away the pseudoscalar would be Goldstone boson, we are left with the doublet-like neutral pseudoscalar, and three singlet-like pseudoscalar states.

In what follows, we will concentrate on the study of the properties of the SM-like Higgs, given the amount of experimental data, and on the ones of the new states with respect to the MSSM, i.e. the right sneutrino-like states. Although not relevant for accommodating the SM-like Higgs, for completeness we will also review the left sneutrino states and the charged Higgs sector.

3.1 The SM-like Higgs

As explained before, we can focus on the analysis of the Higgs-right sneutrino mass submatrix of Appendix A.1 to accommodate the SM-like Higgs in the \(\mu \nu \mathrm{SSM}\). Through the mixing with the right sneutrinos, which appears through \(\lambda _i\), the tree-level mass of the lightest doublet-like Higgs receives an extra contribution with respect to the MSSM. We want to emphasize that this analysis has a notable similarity with that of the NMSSM (although in the NMSSM one has only one singlet), however RPV and an enhanced particle content offer a novel and unconventional phenomenology for the Higgs-right sneutrino sector of the \(\mu \nu \mathrm{SSM}\) [19,20,21,22, 26, 37, 38, 48,49,50]. Taking into account all the contributions, the mass of the SM-like Higgs in the \(\mu \nu \mathrm{SSM}\) can be schematically written as [22, 26]:

$$\begin{aligned} m^{2}_{h} = m^2_{0h} + \Delta _{\text {mixing}} + \Delta _{\text {loop}}, \end{aligned}$$
(21)

where

$$\begin{aligned} m^{2}_{0h} = m^2_Z \cos ^2 2\beta + ({v}/{\sqrt{2}})^2\ {\varvec{\lambda }}^2 \sin ^2 2\beta , \end{aligned}$$
(22)

corresponds to neglect the mixing of the SM-like Higgs with the other states in the mass squared matrix, \(\Delta _{\text {mixing}}\) encodes those mixing effects lowering (raising) the mass if it mixes with heavier (lighter) states, and \(\Delta _{\text {loop}}\) refers to the radiative corrections. Note that \(m^{2}_{0h}\) contains two terms, where the first is characteristic of the MSSM and the second of the \(\mu \nu \mathrm{SSM}\) with

$$\begin{aligned} {\varvec{\lambda }}\equiv \left( {\sum _i \lambda ^2_i}\right) ^{1/2}=\sqrt{3}\ \lambda , \end{aligned}$$
(23)

where the last equality is obtained if one assumes universality of the parameters \(\lambda _i=\lambda \).

We can write \(m^{2}_{0h}\) in a more elucidate form for our discussion below as

$$\begin{aligned} m^{2}_{0h} = {m^2_Z} \left\{ \left( \frac{1 {-} \mathrm{tan}^2\beta }{1 {+} \mathrm{tan}^2\beta }\right) ^2 {+} \left( \frac{v/\sqrt{2}}{m_Z}\right) ^2\ {\varvec{\lambda }}^2 \left( \frac{\mathrm{2\ tan\beta }}{1 {+} \mathrm{tan}^2\beta }\right) ^2 \right\} , \nonumber \\ \end{aligned}$$
(24)

where the factor \(({v/\sqrt{2} m_Z})^2\approx 3.63\), and we see straightforwardly that the second term grows with small tan\(\beta \) and large \({\varvec{\lambda }}\). In the case of the MSSM this term is absent, hence the maximum possible tree-level mass is about \(m_Z\) for \(\mathrm{tan}\beta \gg 1\) and, consequently, a contribution from loops is essential to reach the target of a SM-like Higgs in the mass region around 125 GeV. This contribution is basically determined by the soft parameters \(T_{u_3}, m_{{\widetilde{u}}_{3R}}\) and \(m_{\widetilde{Q}_{3L}}\). On the contrary, in the \(\mu \nu \)SSM one can reach this mass solely with the tree-level contribution for large values of \({\varvec{\lambda }}\) [26]. Following the work of Ref. [22], we choose for this analysis three regions in \({\varvec{\lambda }}\) values. In particular, for convenience of the discussion of Sect. 5, where the last equality of Eq. (23) is used, our regions are:

(a) Small to moderate (\(0.01 \le {\varvec{\lambda }}/\sqrt{3} < 0.2\))

In this range, the maximum value of \(m_{0h}\) using Eq. (24) with \(\varvec{\lambda }/\sqrt{3}=0.2\) goes as \(\approx 78.9\) GeV for \(\mathrm{tan}\beta =2\), which is \(\approx 18\) GeV more compared to a similar situation in the MSSM. It is thus essential to have additional contributions to raise \(m^{2}_{0h}\) up to around 125 GeV. A possible source of extra tree-level mass can arise when the right sneutrinos are lighter compared to the lightest doublet-like Higgs. In this situation, the later feels a push away effect from the former states characterized by \(\Delta _{\text {mixing}}>0\), pushing \(m_h\) a bit further towards 125 GeV. Unfortunately, for most of this range of \(\varvec{\lambda }\) values, away from the upper end, the push-up effect is normally small owing to the small singlet-doublet mixing which is driven by \(\lambda _i\) (see Eq. A.1.4 of Appendix) [26, 49]. The additional contribution to accommodate the 125 GeV doublet-like Higgs is coming then from loop effects. The situation is practically similar to that of the MSSM, where large masses for the third-generation squarks and/or a large trilinear A-term are essential [51,52,53]. A small trilinear \(A_{u_3}\)-term is possible only by decoupling the scalars to at least 5 TeV [53]. A light third generation squark, especially a stop, is natural in the so-called maximal mixing scenario [54], where \(|X_{u_3}/m_{\widetilde{Q}_{3L}}| \approx \sqrt{6}\) with \(X_{u_3}\equiv A_{u_3} - \mu /\tan \beta \).

These issues indicate that the novel signatures from SUSY particles (e.g., from a light stop or sbottom) are less generic in this region of \(\varvec{\lambda }\). Nevertheless, novel differences are feasible for Higgs decay phenomenology, especially in the presence of singlet-like lighter states [20,21,22, 37, 38, 48,49,50].

(b) Moderate to large (\(0.2 \le \varvec{\lambda }/\sqrt{3} < 0.5\))

For this range of \(\varvec{\lambda }\) values, \(m_{0h}\) can go beyond \(m_Z\), especially for \(\mathrm{tan}\beta \lesssim 5\) and \(\varvec{\lambda }/\sqrt{3} \gtrsim 0.29\). For example with \(\varvec{\lambda }/\sqrt{3}\approx 0.4\), \(\mathrm{tan}\beta =2~(5)\) gives \(m_{0h}\sim 112~(96)\) using Eq. (24). This is \(\approx 100\%~(14\%)\) enhancement compared to the MSSM scenario with the same \(\mathrm{tan}\beta \). Note that this value \(\varvec{\lambda }\approx 0.7\) (\(\lambda \approx \) 0.4) is the maximum possible value of \({\varvec{\lambda }}\) maintaining its perturbative nature up to the scale of a grand unified theory (GUT), \(M_{\text {GUT}}\sim 10^{16}\) GeV. As discussed in Ref. [26], using the renormalization group equations (RGEs) for \(\lambda \) and \(\kappa \) between \(M_{\text {GUT}}\) and the low scale \(\sim \) 1 TeV, neglecting the contributions from the top and gauge couplings, one can arrive straightforwardly to the simple formula: \(2.35\ \kappa ^2 + 1.54\ \varvec{\lambda }^2\lesssim 1\). This gives the bound \(\varvec{\lambda }\lesssim 0.8\), similar but slightly larger than the one of 0.7 mentioned before. Nevertheless, one should expect a final bound slightly stronger when all contributions to the RGEs are taken into account. The numerical analysis indicates that a better approximate formula is

$$\begin{aligned} 2.77\ \kappa ^2 + 2\ \varvec{\lambda }^2\lesssim 1, \end{aligned}$$
(25)

which produces the bounds \(\varvec{\lambda }\lesssim 0.7\) and \(\kappa \lesssim 0.6\).

For this region of \(\varvec{\lambda }\) the singlet-doublet mixing is no longer negligible as we will see in the next subsection, particularly as \(\varvec{\lambda }/\sqrt{3}\rightarrow 0.4\). Thus, a state lighter than 125 GeV with the leading singlet composition appears difficult without a certain degree of tuning of the other parameters, e.g. \(\kappa _i\), \(v_{iR}\), \(T_{\kappa _i}\), \(T_{\lambda _{i}}\), etc. In this situation, the extra contribution to the tree-level value of \(m_h\) is favourable through a push-up action from the singlet states compared to small to moderate \(\varvec{\lambda }\) scenario. Once again a contribution from the loops is needed to reach the 125 GeV target. However, depending on the values of \(\varvec{\lambda }\) and \(\mathrm{tan}\beta \) the requirement sometime is much softer compared to small to moderate \(\varvec{\lambda }\) scenario. Thus the necessity of very heavy third-generation squarks and/or large trilinear soft-SUSY breaking term may not be so essential for this region [51]. It is also worth noticing that the naturalness is therefore improved with respect to the MSSM or smaller values of \(\varvec{\lambda }\).

When \(\varvec{\lambda }/\sqrt{3}\rightarrow 0.5\), \(m_{0h}\) as evaluated from Eq. (24) can be larger than 125 GeV. For example, for the upper bound of this range \(\varvec{\lambda }/\sqrt{3}= 0.5\), \(\mathrm{tan}\beta = 2\) gives \(m_{0h}\sim 132\) GeV. In this case we have to relax the idea of perturbativity up to the GUT scale, as we will discuss below.

(c) Large \({\varvec{\lambda }}\) (\(0.5 \le \varvec{\lambda }/\sqrt{3} < 1.2\))

Assuming e.g. a scale of new physics around \(10^{11}\) GeV, and following similar analytical arguments as above using the RGEs, the perturbative limit gives approximately \(1.48\ \kappa ^2 + 0.96\ \varvec{\lambda }^2\lesssim 1\) producing the bounds \(\varvec{\lambda }\lesssim 1\) and \(\kappa \lesssim 0.82\). For \(\varvec{\lambda }\), taking into account the contributions from the top and gauge couplings to the RGEs as above, one can find numerically [26] the final bound \(\varvec{\lambda }\lesssim 0.88\), i.e. \(\varvec{\lambda }/\sqrt{3} \lesssim 0.5\).

Pushing the scale of new physics further below to 10 TeV, the approximate analytical analysis gives the perturbative limit

$$\begin{aligned} 0.25\ \kappa ^2 + 0.14\ \varvec{\lambda }^2\lesssim 1, \end{aligned}$$
(26)

producing now the bounds \(\varvec{\lambda }\lesssim 2.6\) (i.e. \(\varvec{\lambda }/\sqrt{3}\lesssim 1.5\)) and \(\kappa \lesssim 2\). Given that the full numerical analysis produces typically stronger bounds, we will use \(\varvec{\lambda }/\sqrt{3}\lesssim 1.2\) in our scan of Sect. 5. A similar scenario in the context of the NMSSM has been popularised as \(\lambda \)-SUSY [55]. The constraint in this case [56] is slightly different than ours because of the presence of only one singlet.

In this region of \(\varvec{\lambda }\) values, \(m_{0h}\) as evaluated from Eq. (24) can remain well above 125 GeV even up to \(\mathrm{tan}\beta \sim 8\) for \(\varvec{\lambda }/\sqrt{3}\sim 1.2\). For \(\varvec{\lambda }/\sqrt{3} =0.58\), \(m_{0h}\) for \(\mathrm{tan}\beta =2,5\) and 10 is estimated as \(\sim 150\) GeV, 108 GeV and \(\sim 96\) GeV, respectively. With \(\varvec{\lambda }/\sqrt{3}=1.2\) these numbers increase further, for example, \(\sim 113\) GeV when \(\mathrm{tan}\beta =10\). The requirement of an extra contribution to reach the target of 125 GeV is thus rather small and even negative in this corner of the parameter space unless \(\mathrm{tan}\beta \) goes beyond 10 or 15 depending on the values of \(\varvec{\lambda }\). A singlet-like state lighter than 125 GeV is difficult in this corner of the parameter space due to the large singlet-doublet mixing. In fact even if one manages to get a scalar lighter than 125 GeV with parameter tuning, a push-up action can produce a sizable effect to push the mass of the lightest doublet-like state beyond 125 GeV, especially for \(\mathrm{tan}\beta \lesssim 10\) taking \(\varvec{\lambda }/\sqrt{3}=1.2\). Moreover, a huge doublet component makes these light states hardly experimentally acceptable. In this region of the parameter space a heavy singlet-like sector is more favourable which can push \(m_h\) down towards 125 GeV, due to \(\Delta _{\text {mixing}}<0\). In addition, for such a large \(\varvec{\lambda }\) value, new loop effects from the right sneutrinos proportional to \(\varvec{\lambda }^2\) can also give a sizeable negative contribution [28, 57]. A set of very heavy singlet-like states, even with non-negligible doublet composition is also experimentally less constrained.

It is needless to mention that the amount of the loop correction is much smaller in this region compared to the two previous scenarios. Following the above discussion (b) for large values of \(\varvec{\lambda }\), this region of the parameter space also favours third-generation squarks lighter than 1 TeV, which can be produced with enhanced cross sections and can lead to novel signatures of this model with RPV at the LHC, even when the singlet-like states remain heavier, as stated earlier.

In Sect. 5, we will analyze these \(\varvec{\lambda }\) regions using three scans, and we will check how much room is left for new physics in the light of the current precise measurements of the SM-like Higgs properties. Let us now study the right sneutrino-like sector which, as pointed out before, is crucial to determine the properties of the SM-like Higgs in the \(\mu \nu \mathrm{SSM}\).

3.2 The right sneutrino-like states

From the scalar and pseudoscalar mass submatrices in Appendix A, it is clear that \(\kappa _i\) and \(T_{\kappa _i}\) are crucial parameters to determine the masses of the singlet-like states, originating from the self-interactions. The remaining parameters \(\lambda _i\) and \(T_{\lambda _i}\) (\(A_{\lambda _i}\) assuming the supergravity relation \(T_{\lambda _i}= \lambda _i A_{\lambda _i}\) of Eq. 3) not only appear in the said interactions, but also control the mixing between the singlet and the doublet states and hence, contribute in determining the mass scale. Note that the contributions of the parameters \(T_{\nu _i}\) are negligible assuming \(T_{\nu _i}=Y_{\nu _i} A_{\nu _i}\), given the small values of neutrino Yukawas. We conclude, taking also into account the discussion below Eq. (24), that the relevant independent low-energy parameters in the Higgs-right sneutrino sector are the following subset of parameters of Eqs. (9) and (11):

$$\begin{aligned} \lambda _i, \kappa _i, \tan \beta , v_{iR}, T_{\kappa _i}, T_{\lambda _i}, T_{u_3}, m_{{\widetilde{u}}_{3R}}, m_{\widetilde{Q}_{3L}}. \end{aligned}$$
(27)

In the limit of vanishingly small \(\lambda _i\) (considering simultaneously very large \(v_{iR}\) in the case that we require the lighter chargino mass bound of RPC SUSY \(\mu \gtrsim 100\) GeV), not only the off-diagonal entries of the right sneutrino submatrices (A.1.6) and (A.1.12) are negligible, but also the off-diagonal entries (A.1.4), (A.1.5) and (A.1.10), (A.1.11) of the Higgs-right sneutrino matrices. As a consequence of the latter, the singlet states are decoupled from the doublets. It is thus apparent, that \(\lambda _i\) are undoubtedly the most relevant parameters for the analysis of these states. Another aspect of the parameters \(\lambda _i\), namely to yield additional contributions to the tree-level SM-like Higgs mass has already been discussed in the previous subsection. Thus, one can write the right sneutrino masses as

$$\begin{aligned} m_{\widetilde{\nu }^{{\mathcal {R}}}_{iR} }^2= & {} \left( \frac{T_{{\kappa }_i}}{\kappa _i} +2{\mathcal M}_{i}\right) \frac{{{\mathcal {M}}}_{i}}{2}, \end{aligned}$$
(28)
$$\begin{aligned} m_{\widetilde{\nu }^{{\mathcal {I}}}_{iR} }^2= & {} -\frac{3}{2} \frac{T_{{\kappa }_i}}{\kappa _i} {{\mathcal {M}}}_{i}, \end{aligned}$$
(29)

where in the case of supergravity we can use the relation \(T_{\kappa _i}/\kappa _i= A_{\kappa _i}\). In addition, also in this limit \({{\mathcal {M}}}_{i}\) coincide approximately with the masses of the right-handed neutrinos, since they are decoupled from the other entries of the neutralino mass matrix:

$$\begin{aligned} m_{\nu _{iR}} = {{\mathcal {M}}}_{i}. \end{aligned}$$
(30)

With the sign convention adopted in Sect. 2, \({{\mathcal {M}}}_{i}>0\), and from Eq. (29) we deduce that negative values for \(T_{\kappa _i}\) (or \(A_{\kappa _i}\)) are necessary in order to avoid tachyonic pseudoscalars. Using also that equation, we can write (28) as

$$\begin{aligned} m_{\widetilde{\nu }^{{\mathcal {R}}}_{iR}}^2= & {} {{\mathcal {M}}}^2_{i} - \frac{1}{3} m_{\widetilde{\nu }^{{\mathcal {I}}}_{iR} }^2. \end{aligned}$$
(31)

Thus, the simultaneous presence of non-tachyonic scalars and pseudoscalars implies that [22]

$$\begin{aligned} m_{\widetilde{\nu }^{{\mathcal {I}}}_{iR}}< & {} \sqrt{3} {{\mathcal {M}}}_{i}, \end{aligned}$$
(32)
$$\begin{aligned} m_{\widetilde{\nu }^{{\mathcal {R}}}_{iR}}< & {} {{\mathcal {M}}}_{i}. \end{aligned}$$
(33)

Hence, light scalars/pseudoscalar states are assured when light neutralinos (i.e. basically the product \(\kappa _i v_{iR}\)) are present. From Eq. (28) we also see that the absence of tachyons implies the condition

$$\begin{aligned} \frac{-T_{{\kappa }_i}}{\kappa _i} < 2 {{\mathcal {M}}}_{i}, \end{aligned}$$
(34)

and therefore the value of \(T_{{\kappa }_i}/\kappa _i\) and the product \(\kappa _i v_{iR}\) have to be chosen appropriately to fulfill it. Also from that equation we see that singlet scalars lighter than the SM-like Higgs can be obtained when

$$\begin{aligned} 2 {{\mathcal {M}}}_{i} - \frac{ 2m^2_{\text {Higgs}} }{{\mathcal M}_{i}}< \frac{-T_{{\kappa }_i}}{\kappa _i}. \end{aligned}$$
(35)

Thus, for a given value of \({{\mathcal {M}}}_{i}\), only a narrow range of values of \({-T_{{\kappa }_i}}/{\kappa _i}\) is able to fulfill simultaneously both conditions (34) and (35). We will come back to this issue in Sect. 5.

On the other hand, even in a region of small to moderate \(\lambda _i\), to obtain approximate analytical formulas for tree-level scalar and pseudoscalar masses turn out to be rather complicated due to the index structure of the parameters involved. As discussed in detail in Ref. [22], the expressions for their masses can be simplified in the limit of complete degeneracy in all relevant parameters as in Eq. (12), i.e. when \(\lambda _i=\lambda \) (\(\lambda = \varvec{\lambda }/ \sqrt{3}\) as defined in Eq. 23), \(\kappa _i=\kappa , v_{iR}=v_R, T_{\kappa _i}=T_\kappa , T_{\lambda _i}=T_\lambda \). In this case, the \(3\times 3\) scalar and pseudoscalar mass submatrices in Eqs. (A.1.6) and (A.1.12) have the form

$$\begin{aligned} \begin{pmatrix} a&{}\quad b&{}\quad b\\ b&{}\quad a&{}\quad b\\ b&{}\quad b&{}\quad a\\ \end{pmatrix}, \end{aligned}$$
(36)

with the three eigenvalues given by \(a-b\), \(a-b\) and \(a+2b\). Then, it was shown that for both, scalars and pseudoscalars, the two mass eigenstates corresponding to the first two eigenvalues \(a-b\) get decoupled and remain as pure singlet-like states without any doublet contamination. Using the values of a and b from Appendix A, and neglecting \(T_\nu \) under the supergravity assumption of being proportional to the small \(Y_\nu \), one obtains the following degenerate masses:

$$\begin{aligned} m^2_{\widetilde{{\nu }}^{{\mathcal {R}}}_{1,2R}}= & {} \left( \frac{T_{{\kappa }}}{\kappa } +2{\mathcal M}\right) \frac{{{\mathcal {M}}}}{2} + 3\lambda ^2 \left( \frac{v}{\sqrt{2}}\right) ^2 \nonumber \\&\quad \left( \frac{1}{\mu }\frac{T_\lambda }{\lambda } \frac{\tan \beta }{1+\tan ^2\beta }-1\right) , \end{aligned}$$
(37)
$$\begin{aligned} m^2_{\widetilde{{\nu }}^{{\mathcal {I}}}_{1,2R}}= & {} -\frac{3}{2} \frac{T_{{\kappa }}}{\kappa } {{\mathcal {M}}} + 3\lambda ^2 \left( \frac{v}{\sqrt{2}}\right) ^2 \nonumber \\&\quad \times \left[ \left( \frac{1}{\mu }\frac{T_\lambda }{\lambda } + \frac{4}{ 3}\frac{\kappa }{\lambda } \right) \frac{\tan \beta }{1+\tan ^2\beta } -1\right] , \end{aligned}$$
(38)

where now \(\mu \) is defined in Eq. (14), and the Majorana mass is given in Eq. (13) corresponding approximately also to two pure right-handed neutrino states decouple from the rest of the neutralinos:

$$\begin{aligned} m_{\nu _{1,2R}} = {{\mathcal {M}}}. \end{aligned}$$
(39)

The degeneracy of these states can be broken by introducing mild splittings in \(\kappa _i\) values [22, 49], as it is obvious e.g. from Eqs. (28), (29) and (30). One thing to highlight is that the first term of Eq. (37) is like the one of Eq. (28) in the case of universality, and therefore a condition similar to (34) is welcome to avoid tachyons in the scalar spectrum. Nevertheless, depending on the values chosen for the input parameters, the second term in (37) can be positive, relaxing this condition. The latter is especially true for large values of \(\lambda \).

The mass eigenstate corresponding to the third eigenvalue, namely the one which goes as \(a+2b\), however mixes with the doublet-like states, and eventually its mass appears with a complicated form. In the case of the pseudoscalar this is given by

$$\begin{aligned} m^2_{\widetilde{{\nu }}^{{\mathcal {I}}}_{3R}} = - \frac{3}{2}\frac{T_{\kappa }}{\kappa }{\mathcal {M}} + 9\lambda \kappa \ \left( \frac{v}{\sqrt{2}}\right) ^2 \frac{T_{\lambda }/{\lambda }}{\frac{T_\lambda }{\lambda } + \frac{{{\mathcal {M}}}}{2}}\ \frac{\tan \beta }{1+\tan ^2\beta }. \nonumber \\ \end{aligned}$$
(40)

In the case of the scalar, \(m^2_{\widetilde{{\nu }}^{{\mathcal {R}}}_{3R}}\) appears with a much complicated form that can be found in Ref. [22]. Similarly, the third right-handed neutrino-like state mixes with the other MSSM-like neutralinos, and its mass is given by

$$\begin{aligned}&m_{{\nu }_{3R}} = {{\mathcal {M}}} -\frac{\lambda ^2}{2\mu } \left( \frac{v}{\sqrt{2}}\right) ^2\nonumber \\&\quad \times \left[ \frac{1+\tan ^2\beta }{\tan \beta } {-} \frac{4 M \mu }{\left( {v}/{\sqrt{2}}\right) ^{2}} \right] \left[ \frac{1+\tan ^2\beta }{\tan \beta } \frac{M \mu }{ \left( {v}/{\sqrt{2}}\right) ^{2}}-1 \right] ^{-1}, \nonumber \\ \end{aligned}$$
(41)

where M is defined in Eq. (18).

Depending on the input parameters chosen, the two degenerate states can be heavier or lighter than the third state. For example, for the pseudoscalars we see that the second term in Eq. (40) is always positive whereas the second term in Eq. (38) can be positive (larger or smaller than the previous one) or even negative.

It is also worthy to note that for further small \(\lambda \) values (i.e. \(\lesssim 0.01\)) or in the limit of a vanishingly small \(\lambda \), these formulas take simpler forms, and the three states are degenerate:

$$\begin{aligned} m_{\widetilde{\nu }^{{\mathcal {R}}}_{iR} }^2= & {} \left( \frac{T_{{\kappa }}}{\kappa } +2{\mathcal M}\right) \frac{{\mathcal {M}}}{2}, \end{aligned}$$
(42)
$$\begin{aligned} m_{\widetilde{\nu }^{{\mathcal {I}}}_{iR} }^2= & {} -\frac{3}{2} \frac{T_{{\kappa }}}{\kappa } {{\mathcal {M}}}, \end{aligned}$$
(43)
$$\begin{aligned} m_{\nu _{iR}}= & {} {{\mathcal {M}}}. \end{aligned}$$
(44)

As expected, they coincide with Eqs. (28), (29), and (30), respectively, when written for universal parameters.

Other simple formulas can be obtained in the limit of \(\tan \beta \rightarrow \infty \), where

$$\begin{aligned} m^2_{\widetilde{{\nu }}^{{\mathcal {R}}}_{1,2R}}= & {} \left( \frac{T_{{\kappa }}}{\kappa } +2{\mathcal M}\right) \frac{{{\mathcal {M}}}}{2} -3\lambda ^2 \left( \frac{v}{\sqrt{2}}\right) ^2, \end{aligned}$$
(45)
$$\begin{aligned} m^2_{\widetilde{{\nu }}^{{\mathcal {I}}}_{1,2R}}= & {} - \frac{3}{2}\frac{T_{\kappa }}{\kappa }{\mathcal {M}} -3\lambda ^2 \left( \frac{v}{\sqrt{2}}\right) ^2, \end{aligned}$$
(46)

and

$$\begin{aligned} m^2_{\widetilde{{\nu }}^{{\mathcal {R}}}_{3R}}= & {} \left( \frac{T_{{\kappa }}}{\kappa } +2{\mathcal M}\right) \frac{{{\mathcal {M}}}}{2} -3\lambda ^2 \left( \frac{v}{\sqrt{2}}\right) ^2 \left( \frac{2\mu }{m_Z}\right) ^2, \end{aligned}$$
(47)
$$\begin{aligned} m^2_{\widetilde{{\nu }}^{{\mathcal {I}}}_{3R}}= & {} - \frac{3}{2}\frac{T_{\kappa }}{\kappa }{\mathcal {M}}, \end{aligned}$$
(48)

whereas for the right-handed neutrinos one obtains

$$\begin{aligned} m_{\nu _{1,2R}}= & {} {\mathcal {M}}, \end{aligned}$$
(49)
$$\begin{aligned} m_{{\nu }_{3R}}= & {} {{\mathcal {M}}} -\frac{ \lambda ^2 }{2\mu ^2M} \left( \frac{v}{\sqrt{2}}\right) ^4. \end{aligned}$$
(50)

It is evident from this result that unless \(\lambda \) is small to moderate, it is in general hard to accommodate a complete non-tachyonic light spectrum (i.e. \(\lesssim m_{\text {Higgs}}/2\)) for both the scalars and pseudoscalars in the limit of large \(\tan \beta \) without a parameter tuning. In addition, this limit is severely constrained from diverse experimental results. This is because the BRs for some low-energy processes (e.g. \(B^0_s\rightarrow \mu ^+\mu ^-\)), depending on the other relevant parameters are sensitive to the high powers of \(\tan \beta \) and thus, can produce large BRs for these processes in an experimentally unacceptable way. The other limit, i.e. small \(\tan \beta \), on the contrary, is useful from the viewpoint of raising the mass of the lightest doublet-like scalar towards 125 GeV, specially for moderate to large \(\lambda \) values . However, as shown in the discussion of Eqs. (37), (38) and (40), not all the mass formulas for the light states are simple structured in this region and a numerical analysis is convenient.

3.3 The left sneutrinos

The behaviour of the left sneutrinos is very different from the one of the right sneutrinos, since the former are tightly associated to neutrino physics.

As discussed before, the \(3\times 3\) scalar and pseudoscalar left sneutrino submatrices are decoupled from the \(5\times 5\) Higgs-right sneutrino submatrices. Besides, their off-diagonal entries are negligible compared to the diagonal ones, since they are suppressed by terms proportional to \(Y^2_{\nu _{ij}}\) and \(v^2_{iL}\). As a consequence, the mass squared eigenvalues correspond to the diagonal entries, and in this approximation both states also have degenerate masses. Using the minimization equations for \(v_{iL}\), we can write their tree-level values as [19, 26, 38, 41]

$$\begin{aligned} m_{\widetilde{\nu }^{{\mathcal {R}}}_{iL}}^2= & {} m_{\widetilde{\nu }^{{\mathcal {I}}}_{iL}}^2 = \frac{ m_{{{\mathcal {D}}_i}}}{v_{iL}/\sqrt{2}} \frac{v_{iR}}{\sqrt{2}}\nonumber \\&\times \left[ \frac{-T_{{\nu }_i}}{Y_{{\nu }_i}}- \frac{{{\mathcal {M}}}_{i}}{2} + \frac{\mu }{\tan \beta }+\lambda _i \frac{\left( v/\sqrt{2}\right) ^2}{v_{iR}/\sqrt{2}} \frac{\tan \beta }{1+\tan ^2\beta }\right] . \nonumber \\ \end{aligned}$$
(51)

Therefore, in addition to the parameters of Eq. (19) relevant for neutrino physics, the

$$\begin{aligned} T_{\nu _i} \end{aligned}$$
(52)

are relevant parameters for the study of left sneutrino masses. The fourth term in Eq. (51) can usually be neglected as long as \(v_{iR} \gg v\) and/or \(\lambda _i\) is small, and in the limit of moderate/large \(\tan \beta \) one can also neglect the third term. Under these approximations, the condition for non-tachyonic left sneutrinos can be written as an upper bound on the Majorana masses

$$\begin{aligned} \frac{{{\mathcal {M}}}_{i}}{2} \lesssim \frac{-T_{\nu _i}}{Y_{\nu _i}}. \end{aligned}$$
(53)

Given our sign convention of positive Majorana mass, we will use negative values for \(T_{{\nu }_i}\) in the numerical analyses of Sect. 5.

Going back to Eq. (51), we see clearly why left sneutrinos are special in the \(\mu \nu \mathrm{SSM}\) with respect to other SUSY models. Given that their masses are determined by the minimization equations with respect to \(v_i\), they depend not only on left sneutrino VEVs but also on neutrino Yukawas, unlike right sneutrinos, and as a consequence neutrino physics is very relevant for them.

Considering the normal ordering for the neutrino mass spectrum, which is nowadays favored by the analyses of neutrino data [44,45,46,47], and taking advantage of the dominance of the gaugino seesaw for some of the three neutrino families, representative solutions for neutrino/sneutrino physics using diagonal neutrino Yukawas were summarized in Ref. [12]. Different hierarchies among the generations of left sneutrinos are possible, using different hierarchies among \(Y_{\nu _i}\) (and also \(v_{iL}\)).

There is enough freedom in the parameter space of the \(\mu \nu \mathrm{SSM}\) in order to get heavy as well as light left sneutrinos from Eq. (51), and the latter scenario with the left sneutrino as the lightest supersymmetry particle (LSP) was considered in Refs. [11, 12, 41]. Due to the doublet nature of the left sneutrino, masses smaller than half of the mass of the SM-like Higgs were found to be forbidden [12] to avoid dominant decay of the latter into sneutrino pairs, leading to an inconsistency with Higgs data. Let us finally remark that in those works negative values for \(T_{u_3}\) were used, in order to avoid too light left sneutrinos due to loop corrections. Although we are not specially interested in light sneutrinos in this work, we will maintain the same sign convention in what follows. To use positive values for \(T_{u_3}\) would have not modify our results, since their effect on the SM-like Higgs mass is similar.

3.4 The charged scalars

The charged scalars have a \(8\times 8\) (‘charged Higgs’) mass matrix. Similar to the neutral Higgs mass matrices where some sectors are decoupled, the \(2\times 2\) charged Higgs submatrix is decoupled from the \(6\times 6\) slepton submatrix. Thus, as in the MSSM, the mass of the charged Higgs is similar to the one of the doublet-like neutral pseudoscalar, specially when the latter is not very mixed with the right sneutrinos. In this case, both masses are also similar to the one of the heavy doublet-like neutral Higgs.

Concerning the \(6\times 6\) submatrix, the right sleptons are decoupled from the left ones, since the mixing terms are suppressed by the electron-type Yukawa couplings or \(v_{iL}\). Then, the masses of right and left sleptons are basically determined by their corresponding soft terms, \(m_{\widetilde{e}_{iR}}^2\) and \(m_{\widetilde{L}_{i}}^2\), respectively. Although the left sleptons are in the same SU(2) doublet as the left sneutrinos, they are a little heavier than the latter mainly due to the mass splitting produced by the D-term contribution, \(-m_W^2 \cos 2\beta \).

4 Strategy for the scanning

In this section we describe the methodology that we employed to search for points of our parameter space that are compatible with the latest experimental data on Higgs physics In addition, we demanded the compatibility with some flavor observables. To this end, we performed scans on the parameter space of the model, with the input parameters optimally chosen.

4.1 Sampling the \(\mu \nu \mathrm{SSM}\)

For the sampling of the \(\mu \nu \mathrm{SSM}\), we used a likelihood data-driven method employing the MultiNest [34,35,36] algorithm as optimizer. The goal is to find regions of the parameter space of the \(\mu \nu \mathrm{SSM}\) that are compatible with a given experimental data. For it we have constructed the joint likelihood function:

$$\begin{aligned} {\mathcal {L}}_{\text {tot}} ={\mathcal {L}}_{\text {Higgs}} \times {\mathcal {L}}_{\text {B physics}} \times {\mathcal {L}}_{\mu \text { decay}} \times {\mathcal {L}}_{m_{{\widetilde{\chi }}^\pm }}, \end{aligned}$$
(54)

where \({\mathcal {L}}_{\text {Higgs}}\) represents Higgs observables, \({\mathcal {L}}_{\text {B physics}}\) B-physics constraints, \({\mathcal {L}}_{\mu \text { decay}}\) \(\mu \) decay constraints, and \({\mathcal {L}}_{m_{{\widetilde{\chi }}^\pm }}\) constraints on the chargino mass.

To compute the spectrum and the observables we used a suitably modified version of SARAH v4.5.9 code [58] to generate a SPheno v3.3.6 [59, 60] version for the model. We condition that each point is required not to have tachyonic eigenstates. For the points that pass this constraint, we compute the likelihood associated to each experimental data set and for each sample all the likelihoods are collected in the joint likelihood \({\mathcal {L}}_{\text {tot}}\) (see Eq. (54) above).

4.2 Likelihoods

The likelihood functions used have already been discussed in Ref. [12]. Summarizing, we used three types of likelihood functions in our analysis. For observables for which a measurement is available we use a Gaussian likelihood function defined as follows:

$$\begin{aligned} {\mathcal {L}}(x) = \exp \left[ -\frac{(x-x_0)^2}{2\sigma _T^2} \right] , \end{aligned}$$
(55)

where \(x_0 \) is the experimental best fit set on the parameter x, \(\sigma _T^2= \sigma ^2 +\tau ^2\) with \(\sigma \) and \(\tau \) being respectively the experimental and theoretical uncertainties on the observable x.

On the other hand, for any observable for which the constraint is set as lower or upper limit, an example is the chargino lower mass bound, the likelihood function is defined as

$$\begin{aligned} {\mathcal {L}}(x)= & {} \frac{\sigma }{\sigma _T} [ 1-K(D(x)) ] \exp \left[ -\frac{(x-x_0)^2p}{2\sigma _T^2} \right] \nonumber \\&\quad + \frac{1}{ \tau } K((x-x_0)p), \end{aligned}$$
(56)

where

$$\begin{aligned} D(x) = \frac{\sigma }{\tau } \left( \frac{(x_0 - x)p}{ \sigma _T} \right) , \; K(a) = \frac{1}{2} \text {erfc}\bigg (\frac{a}{\sqrt{2}}\bigg ). \end{aligned}$$
(57)

The variable p takes \(+1\) when \(x_0\) represents the lower limit and \(-1\) in the case of upper limit, while erfc is the complementary error function.

The last class of likelihood function we used is a step function in such a way that the likelihood is one/zero if the constraint is satisfied/non-satisfied.

It is important to mention that in this work, unless explicitly mentioned, the theoretical uncertainties \(\tau \) are unknown and therefore are taken to be zero. Subsequently, we present each constraint used in this work together with the corresponding type of likelihood function.

Higgs observables

Before the discovery of the SM-like Higgs boson, the negative searches of Higgs signals at the Tevatron, LEP and LHC, were transformed into exclusions limits that must be used to constrain any model. Its discovery at the LHC added crucial constraints that must be taking into account in those exclusion limits. We have considered all these constraints in the analysis of the \(\mu \nu \mathrm{SSM}\), where the Higgs sector is extended with respect to the MSSM as discussed in Sect. 2. For constraining the predictions in that sector of the model, we interfaced HiggsBounds v5.3.2 [61,62,63,64,65] with MultiNest. First, several theoretical predictions in the Higgs sector (using a \(\pm 3\) GeV theoretical uncertainty on the SM-like Higgs boson) are provided to determine which process has the highest exclusion power, according to the list of expected limits from LEP, Tevatron and LHC. Once the process with the highest statistical sensitivity is identified, the predicted production cross sections of scalars and pseudoscalars are computed in the ‘effective coupling approximation’ of HiggsBounds with the inputs given by SPheno. The cross sections are then multiplied by the BRs and compared with the limits set by these experiments. Finally, whether the corresponding point of the parameter under consideration is allowed or not at 95% confidence level is indicated.

In constructing the likelihood from HiggsBounds constraints, the likelihood function is taken to be a step function. Namely, it is set to one for points for which Higgs physics is realized, and zero otherwise. In order to address whether a given Higgs scalar of the \(\mu \nu \mathrm{SSM}\) is in agreement with the signals observed by ATLAS and CMS, we interfaced HiggsSignals v2.2.3 [66,67,68] with MultiNest. A \(\chi ^2\) measure is used to quantitatively determine the compatibility of the \(\mu \nu \mathrm{SSM}\) prediction with the measured signal strengths and mass. The experimental data used are those of the LHC with some complements from Tevatron. The details of the likelihood evaluation can be found in Refs. [65, 66].

B decays

\(b \rightarrow s \gamma \) is a flavour changing neutral current (FCNC) process, and hence it is forbidden at tree level in the SM. However, its occurs at leading order through loop diagrams. Thus, the effects of new physics (in the loops) on the rate of this process can be constrained by precision measurements. In the combined likelihood, we used the average value of \((3.55 \pm 0.24) \times 10^{-4}\) provided in Ref. [69]. Note that the likelihood function is also a Gaussian (see Eq. 55). Similarly to the previous process, \(B_s \rightarrow \mu ^+\mu ^-\) and \(B_d \rightarrow \mu ^+\mu ^-\) are also forbidden at tree level in the SM but occur radiatively. In the likelihood for these observables Eq. (55), we used the combined results of LHCb and CMS [70]Footnote 3, \(\text {BR} (B_s \rightarrow \mu ^+ \mu ^-) = (2.9 \pm 0.7) \times 10^{-9}\) and \( \text {BR} (B_d \rightarrow \mu ^+ \mu ^-) = (3.6 \pm 1.6) \times 10^{-10}\). Concerning the theoretical uncertainties for each of these observables, we use the guesstimate \(\tau = 10 \%\) of the corresponding best fit value. We denote by \({\mathcal {L}}_{\text {B physics}}\) the likelihood from \(b \rightarrow s \gamma \), \(B_s \rightarrow \mu ^+\mu ^-\) and \(B_d \rightarrow \mu ^+\mu ^-\).

\(\mu \) decays

We also included in the joint likelihood the constraint from BR\((\mu \rightarrow e\gamma ) < 5.7\times 10^{-13}\) [72]Footnote 4 and BR\((\mu \rightarrow eee) < 1.0 \times 10^{-12}\) [74]. For each of these observables we defined the likelihood as a step function. As explained before, if a point is in agreement with the data the likelihood \({\mathcal {L}}_{\mu \text { decay}}\) is set to 1, and otherwise to 0.

Let us point out here that we did not try to explain the interesting but not conclusive 3.5\(\sigma \) discrepancy between the measurement of the anomalous magnetic moment of the muon and the SM prediction, \(\Delta a_{\mu }= a_{\mu }^{\text {exp}}-a_{\mu }^{\text {SM}}=(26.8 \pm 6.3\pm 4.3) \times 10^{-10}\) [8]. Nevertheless, we will check the level of compatibility with this value of the points fulfilling all constraints discussed in Sect. 5.

Chargino mass bound

In RPC SUSY, the lower bound on the lightest chargino mass depends on the spectrum of the model [8, 75]. Although in the \(\mu \nu \mathrm{SSM}\) there is RPV and therefore this constraint does not apply, to compute \({\mathcal {L}}_{m_{{\widetilde{\chi }}^\pm }}\) we have chosen a conservative limit of \(m_{{\widetilde{\chi }}^\pm _1} > 92\) GeV. For the theoretical uncertainty we use the guesstimate \(\tau = 5\%\) of the chargino mass.

4.3 Input parameters

In order to efficiently scan for Higgs physics in the \(\mu \nu \mathrm{SSM}\), it is important to identify first the parameters to be used, and optimize their number and their ranges of values. In Sect. 3, we found that the relevant parameters are those in Eq. (27). However, to perform scans over 19 parameters we would have to run MultiNest a extremely long time making the task very computer resources demanding. The analysis can be nevertheless much simplified assuming universality of the parameters as we did in the discussion below Eq. (35), without significantly modifying the conclusions. In addition, we will also assume in the scans for the sake of simplicity \(m_{{\tilde{Q}}_{3L}} = m_{{\tilde{u}}_{3R}}\). Thus, we will perform scans over the 8 parameters

$$\begin{aligned} \lambda , \kappa , \tan \beta , v_R, T_\kappa , T_\lambda , T_{u_{3}}, m_{{\tilde{Q}}_{3L}}=m_{{\tilde{u}}_{3}}, \end{aligned}$$
(58)
Table 1 Range of low-energy values of the input parameters in Eq. (58) that are varied in the three scans, where \(\tan \beta \) is a flat prior whereas the others are \(\log \) priors. The VEVs \(v_{R}\), and the soft parameters \(T_{\lambda }\), \(T_{\kappa }\), \(T_{u_{3}}\), \( m_{{\tilde{Q}}_{3L}}=m_{{\tilde{u}}_{3R}}\) are given in GeV

as shown in Table 1, using the sign conventions discussed in Sect. 2 and 3 . We will use log priors (in logarithmic scale) for all of the parameters, except for \(\tan \beta \) which is taken to be a flat prior (in linear scale). Let us point out, nevertheless, that we do not assume exact universality of \(\kappa _i\), to avoid an artificial degeneracy in the masses of the two scalars/pseudoscalars (and two neutralinos) which appear in the spectrum without doublet contamination (see the discussion in Sect. 3.2). Thus we take

$$\begin{aligned} \kappa _3 = 1.04\kappa _1, \kappa _2 = 1.02 \kappa _1, \kappa _1=\kappa , \end{aligned}$$
(59)

and scan over \(\kappa \).

For the choice of the scans, we will choose the ranges of \(\varvec{\lambda }\) (\(\equiv \sqrt{3} \lambda \)) discussed in Sect. 3 for convenience of the discussion. In particular, the sample denoted by \(S_1\) corresponds to small/moderate values with \(0.01 \le \lambda < 0.2\), \(S_2\) to moderate/large values with \(0.2 \le \lambda < 0.5\), and finally \(S_3\) to large values \(0.5 \le \lambda < 1.2\). For each scan, the same ranges for the other parameters are considered. In particular, the upper bound of \(\kappa \) has been motivated in the discussion of Sect. 3.1 by relaxing the idea of perturbativity up to the GUT scale, pushing the scale of new physics further below to 10 TeV (see Eq. 26). Concerning the range of \(v_{R}\), the lower and upper bounds allow to have reasonable Majorana masses for right-handed neutrinos, \({{\mathcal {M}}_{i}}=2k_iv_{R}/\sqrt{2}\) (see Eq. 10), even when \(\kappa _i\) are very large or very small, respectively. The ranges of \(T_{\lambda }\) and \(T_{\kappa }\) are also natural following the supergravity framework of Eq. (3). The lower bound on \(m_{{\widetilde{Q}}_{3L}}\) of 200 GeV is chosen to avoid too light stops/sbottoms, and the upper bound of 2 TeV is enough not to introduce too large soft masses and therefore too heavy squarks. With this range of \(m_{{\widetilde{Q}}_{3L}}\), we take the upper bound of \(-T_{u_{3}}\) at 5 TeV to be able to reproduce in the small \(\lambda \) limit the usual maximal mixing scenario when \(m_{{\widetilde{Q}}_{3L}} \sim 2\) TeV.

Table 2 Low-energy values of the input parameters that are fixed in the three scans. The VEVs \(v_{iL}\) and the soft parameters \(T_{u,d,e}\), \(m_{{{\widetilde{Q}}},{{\widetilde{u}}},{\widetilde{d}},{{\widetilde{e}}}}\), \(M_{1,2,3}\) are given in GeV

Let us remark that the MSSM limits on squark masses cannot be applied to the \(\mu \nu \mathrm{SSM}\). For example, if the stop is the LSP it can decay only via RPV channels into top plus neutrino and bottom plus lepton, and these decays can be prompt or displaced depending on the region of the parameter space of the model. Thus, dedicated analyses are necessary for recasting ATLAS and CMS results to the many possible cases of the model, and we leave them for a forthcoming publication [76]. In this sense, in the present work we choose to be conservative enough using the above lower bound of 200 GeV. We consider that with this value we cover all the potentially interesting range of the model. It is also worth noticing that the regions with small stop masses correspond to (light-red and light-blue) points not fulfilling the perturbative condition up to GUT scale, as shown in Figs. 59 and 13 below (this is also evident in Figs. 2028 and 36).

The rest of the parameters of the model, which are less relevant for the analysis, are fixed as shown in Table 2. For squarks, and right sleptons we choose a value of 1000 GeV. Note that the rest of soft masses for Higgses, right sneutrinos and left sleptons, are fixed by the minimization conditions, as discussed in Sect. 2. The relations among gaugino masses \(M_{1,2,3}\) are inspired by GUTs, and in particular we choose gluinos masses of 2.7 TeV. As for the other trilinear parameters, the values of \(T_{d_3}\) and \(T_{e_3}\) have been chosen taking into account the supergravity relations and the corresponding Yukawa couplings. Finally, the parameters \(Y_{{\nu }_i}\), \(v_{iL}\), and \(T_{\nu _i}\) are mainly determined by neutrino and sneutrino physics (see Eqs. 15 and 51).

Since reproducing neutrino data is an important asset of the \(\mu \nu \mathrm{SSM}\), a few words on the subject are worth it. As explained in Sect. 2, how the model reproduces the correct neutrino masses and mixing angles has been intensively addressed in the literature [12, 19, 39, 40, 77]. Although the parameters in Eq. (19), \(\lambda _i\), \(\kappa _i\), \(v_{iR}\), \(\tan \beta \), \(Y_{\nu _i}\), \(v_{iL}\) and M, are important for neutrino physics, the most crucial of them are \(Y_{\nu _i}\), \(v_{iL}\) and M, and they are essentially decoupled from the parameters in Eq. (27) controlling Higgs physics. Thus, for a suitable choice of \(\lambda _i\), \(\kappa _i\), \(v_{iR}\) and \(\tan \beta \) reproducing Higgs physics, there is still enough freedom to reproduce in addition neutrino data by playing with \(Y_{\nu _i}\), \(v_{iL}\) and M, as shown in Ref. [12]. As a consequence, we will not scan over the parameters \(Y_{\nu _i}\), \(v_{iL}\), \(M_1\), \(M_2\) in order to relax our already demanding computing task, and since it is not going to affect our results. For our purposes, it will be sufficient to choose these parameters mimicking the type of solutions of neutrino physics with normal ordering found in Ref. [12], imposing only the cosmological upper bound on the sum of the masses of the light active neutrinos given by \(\sum m_{\nu _i} < 0.12\) [78].

The same comment applies to the parameters \(T_{\nu _i}\) in Eq. (52), which are only relevant to determine the left sneutrino masses, and therefore we fix them to mimic also the left sneutrino physics of Ref. [12]. In that work, it was easy for \(M>0\) to find solutions with the gaugino seesaw as the dominant one for the third family. In this case, \(v_{3L}\) determines the corresponding neutrino mass and \(Y_{\nu _3}\) can be small. On the other hand, the normal ordering for neutrinos determines that the first family dominates the lightest mass eigenstate implying that \(Y_{\nu _{1}}< Y_{\nu _{2}}\) and \(v_1 < v_2,v_3\), with both \(\nu _{R}\)-Higgsino and gaugino seesaws contributing significantly to the masses of the first and second family. Taking also into account that the composition of these two families in the second mass eigenstate is similar, we expect \(v_2 \sim v_3\). Concerning left sneutrino physics, a light tau left sneutrino was required in Ref. [12] implying \(-T_{\nu _{3}} < -T_{\nu _{2}}=-T_{\nu _{1}}\). This pattern of hierarchies for \(Y_{\nu _{i}}\), \(v_{iL}\), and \(T_{\nu _{i}}\) is used in Table 2.

5 Results

By using the methods described in the previous section, we evaluate now the constraints on the parameter space of the \(\mu \nu \mathrm{SSM}\).

To find regions consistent with experimental observations we have performed about 160 million of spectrum evaluations in total, and the total amount of computer required for this was approximately 1110 CPU years.

Differently from other studies, in this work the likelihood is used to drive MultiNest. Then, the selection of the viable points is not based on the value of the likelihood but rather on the series of cuts that are applied on the samples, as will be discussed subsequently.

To carry out the analysis, we first demand Higgs physics to be fulfilled. As already mentioned in Sect. 4, we use HiggsBounds and HiggsSignals to take into account the constraints from 7-, 8- and 13-TeV LHC data, as well as those from LEP and Tevatron. In particular, we require that the p-value derived by HiggsSignals be larger than 5%. It is worth noticing here that, with the help of Vevacious [79], we have also checked that the EWSB vacua corresponding to the previous allowed points are viable. Then, we select points that lie within \(3\sigma \) from \(b \rightarrow s \gamma \), \(B_s \rightarrow \mu ^+\mu ^-\), and \(B_d \rightarrow \mu ^+\mu ^-\). In the third step, the points that pass these cuts are required to also satisfy the upper limits of \(\mu \rightarrow e \gamma \) and \(\mu \rightarrow eee\), and the lower bound on the chargino mass inspired in RPC SUSY following Sect. 4.2. After all these cuts are applied, about 91% of the points survive. At last we require all the points that passed the above set of cuts to satisfy the cosmological upper bound on the sum of the masses of the light active neutrinos.

As we will explain below, after imposing the relevant constraints from Higgs physics, only \(b \rightarrow s \gamma \) and (less importantly) the bound on neutrino masses put further constraints on the parameter space of the \(\mu \nu \mathrm{SSM}\). As already mentioned in Sect. 4.2, in our computation we have not tried to explain the discrepancy between the measurement of the muon anomalous magnetic moment and the SM prediction \(\Delta a_{\mu }\). Nevertheless, for completeness, we will discuss the level of compatibility of the SUSY contributions with this value, and possible improvements in this direction.

b \(\rightarrow \) s \(\gamma \)

The BR(\(b\rightarrow s\gamma \)) puts some constraints on the parameters space of the \(\mu \nu \mathrm{SSM}\), as shown in Fig. 1.Footnote 5 There we show the constraints from \(b\rightarrow s\gamma \) for all points of the parameter space fulfilling Higgs physics. For instance, in our setup this BR can be too small in certain regions of the parameter space. Forbidden points occur for small to moderate values of \(\lambda \), such as in \(S_1\) and \(S_2\), when \(\tan \beta \) can be large while \(m_{{\widetilde{Q}}_{3L}}\) can be small. As is well known, the most important contributions to the BR(\(b\rightarrow s\gamma \)) come from chargino/stop and charged Higgs/top mediated processes [81]. On the one hand, the charged Higgs contribution always tends to increase the SM value while that of the charginos depends on the sign of \(M_2\), \(T_{u_{3}}\) and \(\mu \), where in our case \(\mu =3\lambda v_R/\sqrt{2}\). Since we are working with \(M_2, \mu > 0\) and \(T_{u_{3}}<0\), the contribution from charginos in the loops acts destructively. Also for light sparticles (here charginos, charged Higgs and stops) and/or large \(\tan \beta \) the effects can be large. This is actually what happens in our cases. For small/moderate \(\lambda \), large \(\tan \beta \) favors increasing this effect. In the regime of destructive contribution involving light stops (when \(m_{{\widetilde{Q}}_{3L}}\) becomes small or in the maximal mixing scenario) and light Higgsinos (winos are moderately heavy since we fix \(M_2\) to 1800 GeV), this effect is large and suppresses the BR(\(b\rightarrow s\gamma \)). Note that for \(S_3\) this does not occur. The reason is that large values of \(\tan \beta \) are not needed, as we will see in detail in the next subsection, and in addition moderate values come together with relatively large values of \(m_{{\widetilde{Q}}_{3L}}\).

Fig. 1
figure 1

Constraints from \(b \rightarrow s \gamma \) in the \(\tan \beta -\mu \) plane, for scans \(S_{1,2,3}\). The grey (light-grey) color corresponds to points of the parameter space fulfilling Higgs physics that are (are not) compatible with the BR(\(b \rightarrow s \gamma \))

Fig. 2
figure 2

Constraints from \(\sum m_{\nu _i} < 0.12\) eV in the \({\mathcal {M}}-\tan \beta \) plane, for scans \(S_{1,2,3}\). The purple (light-purple) color corresponds to points of the parameter space fulfilling Higgs physics that are (are not) compatible with the cosmological upper bound on the sum of the masses of the light active neutrinos

Sum of neutrino masses

In Fig. 2, we show the constraints on the parameter space fulfilling Higgs physics imposed by the requirement \(\sum m_{\nu _i} < 0.12\) eV in the \({\mathcal {M}}-\tan \beta \) plane, with \({\mathcal {M}}=2\kappa \ v_R/\sqrt{2}\). We find that the sum of the masses of the three light neutrinos can exceed this upper bound when the Majoranna masses are small. This can be qualitatively explained using Eq. (15) with the approximations discussed below Eq. (19). Then, the gaugino seesaw contributions to neutrino masses given by the second term in Eq. (15), with \(M^{\text {eff}}=M\), is fixed in our scans. In particular, using the values of Table 2 for \(v_{iL}\) and \(M = 2640.45\) GeV from the values of \(M_{1,2}\), we can compute these contributions to the diagonal entries of the mass matrix \((m_{\nu })_{ii}\), which turn out to be in absolute value 0.002, 0.015, and 0.0286 eV for \(i=1,2,3\), respectively. This indicates that for sizable \(\nu _{R}\)-Higgsino seesaw, i.e. the first term in Eq. (15), the mass of the heaviest neutrino can easily be made too large. This occurs when \({\mathcal {M}}\) is small. For example, for \(\tan \beta = 10\) and \({\mathcal {M}}=30\) GeV the \(\nu _{R}\)-Higgsino seesaw contribution to the diagonal entries is in absolute value around 0.027, 0.108, and 0.0017 eV, respectively, and added to the gaugino seesaw at least one neutrino mass would be larger than 0.12 eV. Actually, in our scenarios the effect of \(\tan \beta \) is not very relevant, and the size of \({\mathcal {M}}\) is the most important one. In particular, as shown in Fig. 2, in scans \(S_1\), \(S_2\), and \(S_3\), for \({\mathcal {M}}\) below 123, 52, and 51 GeV, respectively, we find points excluded by the cosmological upper bound on neutrino masses.

Fig. 3
figure 3

Analysis of \(a^{\text {SUSY}}_\mu \) in the \(\tan \beta -\mu \) plane for scans \(S_{1,2,3}\). The light-brown color corresponds to points of the parameter space fulfilling Higgs physics compatible at \(2 \sigma \) with \(\Delta a_{\mu }\). Brown color corresponds to compatibility between 2 and \(3 \sigma \), whereas points with black color are compatible within 3 and \(3.5 \sigma \)

Muon anomalous magnetic moment

The difference between the experimental measurement and the SM prediction \(\Delta a_{\mu }= a_{\mu }^{\text {exp}}-a_{\mu }^{\text {SM}} = (26.8 \pm 6.3\pm 4.3) \times 10^{-10}\) [8], where the errors are from experiment and theory prediction (with all errors combined in quadrature), respectively, represents an interesting but not conclusive discrepancy of 3.5 times the combined 1\(\sigma \) error. SUSY contributions \(a^{\text {SUSY}}_\mu \) can be large in the presence of light muon sneutrino and charginos or light neutralino and smuons. We found in our scans \(S_1\), \(S_2\) and \(S_3\) that \(a^{\text {SUSY}}_\mu \) is smaller than \(16.96\times 10^{-10}\), \(16.83\times 10^{-10}\), and \(3.7 \times 10^{-10}\), respectively. Thus, although none of the points of the parameter space is compatible at \(1 \sigma \) with \(\Delta a_{\mu }\) in some regions \(a^{\text {SUSY}}_\mu \) is compatible at \(2 \sigma \). Note that we are neglecting the uncertainties in the SUSY computation. The result is shown in Fig. 3. The largest contributions to \(a^{\text {SUSY}}_{\mu }\) are found for small \(\mu \) and large \(\tan \beta \). In our scenarios, since bino- and wino-like (neutralino or chargino) eigenstates are heavy (in our scans \(M_2=2M_1=1800\) GeV) the contributions involving them are suppressed. Besides, although the Higgsino-like eigenstates can be light when \(\mu \) is relatively small, their contributions can be diluted by the small Yukawa coupling of the muon. Nevertheless, when \(\tan \beta \) is very large this effect can be more important. A way of explaining the discrepancy \(\Delta a_{\mu }\) with \(a^{\text {SUSY}}_\mu \) is to try to lower the muon left sneutrino mass, which in these scans is generically large given the input parameters chosen for neutrino physics. Changing the latter we could obtain smaller masses, and we leave the analysis of this interesting possibility for a forthcoming publication [82].

5.1 Viable regions of the parameter space

Once \(b\rightarrow s\gamma \), and mainly Higgs physics, have determined the parameter space that is viable in the \(\mu \nu \mathrm{SSM}\), we will discuss it in detail. In order to carry it out we will follow the division in the three different scans presented in Sect. 4.3.

5.1.1 Scan 1 (\(0.01 \le \lambda < 0.2\))

Let us concentrate first on the analysis of the results for Scan 1 (\(S_1\)). We show in Fig. 4 the viable points of the parameter space in the \(\kappa - \lambda \) plane. The red points represent cases where the SM-like Higgs boson is the lightest scalar. All of them fulfill perturbativity up to the GUT scale, and therefore \(\kappa \lesssim 0.6\). For the light-red points the SM-like Higgs boson is still the lightest scalar, but we have relaxed the perturbativity condition up to 10 TeV and therefore \(0.6 \lesssim \kappa \lesssim 2\). On the contrary, the blue points represent cases where the SM-like Higgs boson is not the lightest scalar. This figure can be considered as the summary of results for this scan, which we now discuss in detail.

Fig. 4
figure 4

Viable points of the parameter space for \(S_1\) in the \(\kappa -\lambda \) plane. The red and light-red (blue) colours represent cases where the SM-like Higgs is (is not) the lightest scalar. All red and blue points below the lower black dashed line fulfill the perturbativity condition up to GUT scale of Eq. (25). Light-red points below the upper black dashed line fulfill the perturbativity condition up to 10 TeV of Eq. (26)

As shown in the figure, we find viable solutions in almost the entire \(\kappa - \lambda \) plane analyzed in \(S_1\). The only small (white) region that becomes forbidden corresponds to very small values of \(\lambda \) and very large (non-perturbative up to the GUT scale) values of \(\kappa \). This can be understood taking into account that we are asking to all the points to fulfill the chargino mass lower bound of RPC SUSY, which corresponds to condition \(\mu =3\lambda v_{R}/\sqrt{2} \gtrsim 100\) GeV. Thus for a small \(\lambda \), a large \(v_R\) is needed (see also Fig. 16 in Appendix BFootnote 6). However, this gives rise to a large value of \({{\mathcal {M}}} ={2}\kappa {v_{R}/\sqrt{2}}\) and, as a consequence, the condition in Eq. (53) to avoid tachyonic left sneutrinos cannot be fulfilled for any value of \(\kappa \). In particular, combining both conditions we can write \(\frac{100\ \text {GeV}}{3\lambda } \lesssim \frac{v_R}{\sqrt{2}} \lesssim \frac{-T_{\nu _i}/Y_{\nu _i}}{\kappa }\), which cannot always be fulfilled.

This is the case for the muon left sneutrino whose ratio \(-T_{\nu _2}/{Y_{\nu _2}}=2500\) GeV is the smallest of the three families, as can be deduced from Table 2. For example, \(\lambda =\) 0.01 implies \(v_{R}/\sqrt{2} \gtrsim \) 3300 GeV, and then it is straightforward to see that \(\kappa \lesssim 0.75\) to avoid tachyons. Let us point out nevertheless, that this forbidden tachyonic region in Fig. 4 turns out to be an artifact of our simplified assumption about the neutrino (sneutrino) pattern in order to relax the demanding computing task, as discussed in Sect. 4.3. Simply breaking the degeneracy between \(T_{\nu _1}\) and \(T_{\nu _2}\), taking a larger value for \(T_{\nu _2}\), we would recover this region as viable. It is apparent that there is a less point-dense area with \(\kappa \lesssim 0.6\) and \(\lambda \gtrsim 0.15\) (the same occurs for the area with \(\lambda \gtrsim 0.45\) in Fig. 8 to be discussed below). This is just an artifact of the sampling, and it would have been filled out with more computing-time resources.

Going back to the values of \(v_R\) in Fig. 16, it is worth noticing that for large \(\lambda \) and/or large \(\kappa \) they are bounded, \(v_R/\sqrt{2} \lesssim 2000\) GeV. The reason is that for those points, to increase the value of \(v_R\) would increase the mixing term \(m_{H^{{\mathcal {R}}}_{u} H^{{\mathcal {R}}}_{d}}^2\) in Eq. (A.1.3) of Appendix A, decreasing therefore the SM-like Higgs mass, and eventually leading to the appearance of a negative eigenvalue. Note in this sense that the diagonal term \(m_{H^{{\mathcal {R}}}_{u} H^{{\mathcal {R}}}_{u}}^2\) (\(m_{H^{{\mathcal {R}}}_{d} H^{{\mathcal {R}}}_{d}}^2\)) in Eq. (A.1.2) (Eq. A.1.1) is small (large) for the large values of \(\tan \beta \) present in this scan, as we will discuss below (see Fig. 17). The mixing terms with right sneutrinos also increase with the value of \(v_R\), as can be seen in Eqs. (A.1.4) and (A.1.5), but much less than the above between Higgses, since the former go like \(v_R\) whereas the latter as \(v_R^2\). As we can see in those equations, the value of \(T_\lambda \) is also important to determine the mixing among states. In Fig. 18, we see that in most of the regions \(T_\lambda \) has an upper bound of around 200 GeV, and only for the lower right region with large \(\lambda \), but small (perturbative up to the GUT scale) \(\kappa \), it can reach up to 500 GeV. In the region to the left of the latter, although the values of \(\kappa \) are also small, \(v_R\) is large as discussed above, and smaller \(T_\lambda \) is favoured. On the other hand, assuming the supergravity relation \(A_\lambda =T_\lambda / \lambda \), one can check that in most of the regions \(A_\lambda \) has the upper bound of around 2 TeV, as shown in Fig. 19.

Concerning the values of \(\tan \beta \), we find in \(S_1\) that \(\tan \beta > 4\). Such a lower bound is expected in order to maximize the tree-level SM-like Higgs mass for small/moderate values of \(\lambda \), as discussed in Sect. 3.1. We can see in Fig. 17 that large values of \(\tan \beta \) are welcome for this task, similarly to the MSSM. Given the small singlet-doublet mixing, significant loop contributions are the main source to increase the tree-level mass of the SM-like Higgs. The values of the masses of the third-generation squarks and trilinear soft term necessary to generate the large loop corrections are shown in Fig. 5. The white region in the upper left side is excluded by the mass of the SM-like Higgs or by the existence of tachyons when \(-T_{u_3}\) is much larger than \(m_{\widetilde{Q}_{3L}}\). For the allowed regions, we can see first that for \(\kappa \) perturbative up to 10 TeV (light-red points) the values of \(-T_{u_3} \lesssim 2000\) GeV and \(m_{{\widetilde{Q}}_{3L}}\lesssim 1000\) GeV are highly correlated. Given the large value of \(\kappa \), the push-down effect in these light-red points makes necessary the maximal mixing scenario to cancel it, bringing the mass of the SM-like Higgs to the correct value. For the red points, where \(\kappa \) is smaller, the push-down effect is not so large and the maximal mixing scenario can be relaxed. We can see that the lower right side of Fig. 5 becomes populated. The same argument applies to the blue points (note that most of them are on top of red points), but now for the push-up effect which is also small. In Figs. 20 and 21 of Appendix B, we show in the \(\kappa -\lambda \) plane the values of \(m_{\widetilde{Q}_{3L}}\) and \(-T_{u_3}\), respectively. As discussed, smaller values of these parameters are needed in the perturbative region up to 10 TeV.

Fig. 5
figure 5

Viable points of the parameter space for \(S_1\) in the \(-T_{u_3}\) versus \(m_{{\widetilde{Q}}_{3L} }\) plane. The color code is the same as in Fig. 4

In Table 3 of Appendix C, we show the BP S1-R1 corresponding to the red region of Fig. 4, where the SM-like Higgs \(h_1\) is the lightest scalar, and \(h_{4,5,6}\) are the singlet-like states with masses larger than 900 GeV. Note nevertheless that the singlet-like pseudoscalars can be lighter than the SM-like Higgs, as shown in particular in this BP where they have masses around 40 GeV. As we can check from the fifth box of the table, the right sneutrinos are not very mixed among themselves because \(\lambda \) is small and therefore the off-diagonal terms in Eq. (A.1.6) are negligible. However, the singlet-like scalar \(h_6\), with a mass similar to \(h_7\) which has a dominant composition of \(H^{{\mathcal {R}}}_d\), has a significant composition of the latter (22.5%), whereas \(h_4\) and \(h_5\) are very pure singlets with dominant compositions \(\widetilde{{\nu }}^{{\mathcal {R}}}_{eR}\) and \(\widetilde{{\nu }}^{{\mathcal {R}}}_{\mu R}\), respectively. This BP corresponds to one of those shown in the right-hand side of Fig. 6 with a large doublet composition. In that figure, the singlet component of the singlet-like states is shown versus their masses, and we can see that most (but not all) of them are almost pure singlets. The fact that only one of the three singlet-like states of S1-R1 has a large doublet composition is because of our assumption of almost degenerate \(\kappa \)’s implying that there are always two almost pure singlets. In Table 4, we show a different BP corresponding to the red region of Fig. 4, S1-R2, where the three singlet-like scalars are very pure singlets.

Fig. 6
figure 6

The singlet component \(\sum _i|Z^H_{h \widetilde{{\nu }}^{{\mathcal {R}}}_{iR}}|^2\) of the singlet-like scalars h versus their masses, for \(S_1\). We only show here viable points with scalar masses smaller than 1000 GeV. In the lower part we zoom in the low-mass region

It is worth remarking here that the masses of the left sleptons are determined by the parameters related to neutrino/left sneutrino physics, as discussed in Sect. 4.3, and therefore can be modified choosing different values for \(Y_{\nu _i}\), \(v_{iL}\), \(M_{1,2}\) and \(Y_{\nu _i}\) in Table 2. This is true in general for all scans and applies therefore to all BPs studied in this work.

Let us now discuss in more detail the (narrow) region of Fig. 4 with blue points. For small values of \(\lambda \), the first term of Eq. (37) is a good approximation for right sneutrino masses. Clearly, unless one makes a tuning between the two pieces in that term, \(T_\kappa /\kappa \) and \(2{{\mathcal {M}}}=4\kappa v_R/\sqrt{2}\), one needs these two quantities to be small in order to obtain right sneutrinos lighter than the SM-like Higgs. Now, since \(v_R\) is typically large in this scan compared to the SM-like Higgs mass, small values of \(\kappa \) are necessary for this task. This is what we observe in the blue region of Fig. 4, where \(\kappa \lesssim 0.2\). There we also see that for larger values of \(\lambda \), larger values of \(\kappa \) are allowed, because the values of \(v_R\) decrease with \(\lambda \) as shown in Fig. 16. The correlation between the above two pieces for the blue points is also obvious from Eqs. (34) and (35). We show explicitly this effect in Fig. 7, where basically the line \(-T_{\kappa } = 2\kappa {{\mathcal {M}}}\) separates the tachyonic (white) region from the non-tachyonic one with blue and red points, i.e. \(-T_{\kappa }< 2\kappa {{\mathcal {M}}}\). Blue points have to be close to the line since they have to fulfill in addition the approximate condition (35). In Fig. 22 of the Appendix, we show the different values of \(-T_\kappa \) in the \(\kappa -\lambda \) plane. As we can see, for the small values of \(\kappa \) corresponding to the blue region of Fig. 4, the values of \(-T_\kappa \) are typically small. For larger values of \(\kappa \) corresponding to the regions with red and light-red points of Fig. 4, i.e. with masses of the singlet-like states larger than the SM-like Higgs mass, the tachyonic region can be avoided even with large values of \(-T_\kappa \) (up to the upper bound of 500 GeV imposed in the scan), as shown in the figure. We show for completeness in Fig. 23 the different values of the supergravity parameter \(A_\kappa = T_\kappa /\kappa \) in the \(\kappa -\lambda \) plane. Due to this relation, values of \(-A_\kappa \) as large as around 2.9 TeV can be obtained in regions with small \(\kappa \). Larger values of \(-A_\kappa \) are not possible because the condition \(-A_{\kappa }< 2 {{\mathcal {M}}}\) cannot be fulfilled since \(v_R\) is bounded, and therefore tachyons would appear.

Fig. 7
figure 7

Viable points of the parameter space for \(S_1\) in the \(-T_{\kappa }\) versus \(2\kappa {{\mathcal {M}}}\) plane. The color code is the same as in Fig. 4. In the upper right we zoom in the region with blue points

Finally, it is worth noticing that about 40% of the blue points correspond to cases where the singlet-like scalars have masses \(\lesssim m_{\text {Higgs}}/2\). As can be seen in Fig. 6, most of these states are almost pure singlets and therefore do not affect the Higgs decays, surviving as viable points. We show the BP S1-B1 with these characteristics in Table 5. There we see that the three singlet-like states \(h_{1,2,3}\) with masses around 50 GeV are lighter than the SM-like Higgs \(h_4\). Besides, these light states are significantly mixed among themselves because of the moderate value of \(\lambda = 0.1\), which makes the off-diagonal terms in Eq. (A.1.6) significant. It is also worth noticing that for S1-B1 the Majorana mass is small, \({\mathcal {M}}=52.92\) GeV, giving rise to two almost degenerate light right-handed neutrinos of masses 49.75 and 53.28 GeV, and one heavier of mass 64.49 GeV. As a consequence, although the singlet composition of \(h_4\) is small, this mixing is already sufficient to produce a significant decay of \(h_4\) to two neutralinos (with dominant right-handed neutrino composition) with BR=0.15.

The presence of light scalars (\(h_i\)), pseudoscalars (\(A_i\)) and neutralinos (\({\widetilde{\chi }}^0_i\)) such that \(m_{h_i}+m_{A_j} < M_Z\), \(m_{{\widetilde{\chi }}^0_i}+m_{{\widetilde{\chi }}^0_j} < M_Z\) or \(m_{{\widetilde{\chi }}^0_i}+ m_{{\widetilde{\chi }}^\pm _j} < M_W\) (here \({\widetilde{\chi }}^\pm _j =e, \mu , \tau \)) opens up new on-shell decay modes for the Z and W bosons. The possible signs of new physics from these new decay modes in the \(\mu \nu \mathrm{SSM}\) have been studied in Ref. [50].

On the other hand, when the masses of the singlet-like states are close to 125 GeV, it is possible to find solutions with a large doublet composition. Actually, as mentioned before, for each point of the parameter space only one of the three states has this property, given our assumption of almost degenerate \(\kappa \)’s implying that there are always two almost pure singlets. For these solutions, if the SM-like Higgs and the singlet-like state with significant doublet composition have masses within the mass resolution of the experiment, they will have their signal rates superimposed, and both will contribute to the resonance observed at 125 GeV [66]. In this scan \(S_1\), about 0.4% of the phenomenologically viable points found have singlet-like states with masses close to 125 GeV. We show in Table 6 the BP S1-2h1 with these properties. There we see that the right sneutrino \(h_4\) has a large composition of \(H^{{\mathcal {R}}}_u\) (27.58%), whereas \(h_1\) and \(h_2\) are very pure singlets with dominant compositions \(\widetilde{{\nu }}^{{\mathcal {R}}}_{eR}\) and \(\widetilde{{\nu }}^{{\mathcal {R}}}_{\mu R}\), respectively, and not contributing therefore to the superposition of Higgs-like states. As already discussed for the BP S1-R1, the singlet-like states are very little mixed among themselves given that \(\lambda \) is small.

Fig. 8
figure 8

Viable points of the parameter space for \(S_2\) in the \(\kappa -\lambda \) plane. The red and light-red (blue and light-blue) colours represent cases where the SM-like Higgs is (is not) the lightest scalar. All red and blue points below the lower black dashed line fulfill the perturbativity condition up to GUT scale of Eq. (25). Light-red and light-blue points below the upper black dashed line fulfill the perturbativity condition up to 10 TeV of Eq. (26)

5.1.2 Scan 2 (\(0.2 \le \lambda < 0.5\))

Figure 8 summarizes our results for Scan 2 (\(S_2\)). In this case, we find viable solutions in the entire \(\kappa - \lambda \) plane, since now \(\lambda \ge 0.2\) and therefore the chargino mass lower bound can be fulfilled with low values of \(v_R\), being safe from tachyonic left sneutrinos. In fact, we see in Fig. 24 that in most of the regions \(v_R/\sqrt{2} \lesssim 1000\) GeV. This bound, in order to avoid a too large mixing term \(m_{H^{{\mathcal {R}}}_{u} H^{{\mathcal {R}}}_{d}}^2\), is smaller than for \(S_1\) because now we are working with moderate/large values of \(\lambda \). Concerning the value of \(T_\lambda \), in Fig. 26 we see that in regions with large \(\kappa \), and given the values of the other relevant parameters, small values are preferred (for example, for \(\kappa >1\) about 87% of points have \(T_\lambda < 200\) GeV), whereas for lower values of \(\kappa \) the mixing term is smaller and larger values of \(T_\lambda \) are favoured (up to the upper bound of 500 GeV imposed in the scan). Assuming the supergravity relation \(A_\lambda = T_\lambda /\lambda \), we see in Fig. 27 that given the moderate/large values of \(\lambda \) for this scan, the upper bound for \(A_\lambda \) is typically smaller than for \(S_1\) in all regions, with a maximum value of around 2.5 TeV.

Concerning \(\tan \beta \), since \(\lambda \) is larger than in \(S_1\) we find that smaller values are favoured to maximize the tree-level SM-like Higgs mass, as shown in Fig. 25 of Appendix. In addition, as a consequence of the moderate/large \(\lambda \), the singlet–doublet mixing is larger and therefore the push-up effect for the blue points helps to increase the tree-level mass. All these effects together produce that the loop contributions to increase the tree-level mass of the SM-like Higgs can be relaxed. One can observed this comparing Figs. 9, 28 and 29 with the corresponding ones of \(S_1\).

Related to the above discussion, is the fact that in the perturbative region up to the GUT scale is more easy to find blue than red points. The push-down effect of the latter makes for them more difficult to reach the correct mass of the SM-like Higgs. Note also that blue points with all values of \(\kappa \) are present, since \(v_R\) is now smaller than for \(S_1\).

In Tables 7 and 8 of Appendix C, we show the two BPs S2-R1 and S2-R2, respectively, corresponding to the red region of Fig. 8. They have different singlet-like scalar masses, around 230 and 600 GeV, mainly due to the different values of \(v_R\). For S2-R1, as we can see in the fifth box, the singlet-like states are significantly mixed among themselves because of the moderate/large value of \(\lambda \), and the eigenstate \(h_5\) is the one having a significant composition of \(H^{{\mathcal {R}}}_u\) (10.29%). The SM-like Higgs with a composition of \(H^{{\mathcal {R}}}_d\) of 19.36% is phenomenologically viable because \(\tan \beta \) is as small as 2.31. The same occurs for S2-R2, where now the SM-like Higgs composition of \(H^{{\mathcal {R}}}_d\) is larger, 46.66%, but \(\tan \beta \) is smaller, 1.08. For this BP the mixing among right sneutrinos is larger, but no eigenstate has a significant composition of \(H^{{\mathcal {R}}}_u\) given their larger masses. In Fig. 10, we show the singlet component of the singlet-like scalars. As for \(S_1\), for large masses we can find scalars with a very large composition of \(H^{{\mathcal {R}}}_d\).

The correlation discussed for \(S_1\) in Fig. 7 is relaxed in this new scan, again because of the larger values of \(\lambda \), as discussed below Eq. (39). We show this in Fig. 11. In Figs. 30 and 31 of the Appendix, we can see the different values of \(-T_\kappa \) and \(A_\kappa = T_\kappa /\kappa \), respectively, in the \(\kappa -\lambda \) plane. In the perturbative region up to the GUT scale, except for areas with \(\kappa \) close to its upper bound, \(T_\kappa \) is typically small to avoid tachyonic right sneutrinos because \(v_R\) is small. As a consequence, in the case of supergravity \(A_\kappa \) is also typically small in this region. In this scan, about 11% of the blue points correspond to cases where the singlet-like scalars have masses \(\lesssim m_{\text {Higgs}}/2\). In Table 9, we show the BP S2-B1 corresponding to the blue region of Fig. 8, with singlet-like scalar masses \(> m_{\text {Higgs}}/2\). Apart for that, its characteristics are similar to the BP S1-B1 of scan \(S_1\).

Fig. 9
figure 9

Viable points of the parameter space for \(S_2\) in the \(-T_{u_3}\) versus \(m_{{\widetilde{Q}}_{3L} }\) plane. The color code is the same as in Fig. 8

Fig. 10
figure 10

The singlet component \(\sum _i|Z^H_{h \widetilde{{\nu }}^{{\mathcal {R}}}_{iR}}|^2\) of the singlet-like scalars h versus their masses, for \(S_2\). We only show here viable points with scalar masses smaller than 1000 GeV. In the lower part we zoom in the low-mass region

Concerning solutions with singlet-like states with masses close to 125 GeV, about 5% of the phenomenologically viable points found in this scan are of this type. However, not all of them have a significant doublet composition as to have their signals superimposed with that of the SM-like Higgs. We show in Table 10 the BP S2-2h1 as an example of this situation. As we can see, the right sneutrino \(h_4\) has the largest doublet composition of \(H^{{\mathcal {R}}}_{u}\) (3.77%) and \(H^{{\mathcal {R}}}_{d}\) (2.65%), but insufficient as to contribute significantly to the Higgs signals. Unlike the BP S1-2h1 of \(S_1\), now the three sneutrinos are very mixed because of the larger value of \(\lambda \). Similar to S1-B1, for this BP also the Majorana mass is small, \({\mathcal {M}}=55.8\) GeV, giving rise to two almost degenerate light right-handed neutrinos of masses 55.9 and 57.2 GeV, and one heavier of mass 76.7 GeV. As a consequence, the decay channel right sneutrino to two neutralinos (with dominant right-handed neutrino composition) opens, giving the most important contribution to the BRs. In Table 11, we show the BP S2-2h2 with the singlet-like state \(h_3\) having the largest doublet composition of \(H^{{\mathcal {R}}}_{u}\) (8.79%) and contributing significantly to the superposition of signals with the SM-like Higgs \(h_4\), unlike the previous case S2-2h1.

Fig. 11
figure 11

Viable points of the parameter space for \(S_2\) in the \(-T_{\kappa }\) versus \(2\kappa {{\mathcal {M}}}\) plane. The color code is the same as in Fig. 8. In the upper right we zoom in the region with blue points

5.1.3 Scan 3 (\(0.5 \le \lambda < 1.2\))

The results for Scan 3 (\(S_3\)) are summarized in Fig. 12. In this case with so large values of \(\lambda \), the white region in the lower right is forbidden because of the too large mixing term \(m_{H^{{\mathcal {R}}}_{u} H^{{\mathcal {R}}}_{d}}^2\) producing tachyons. To avoid that situation, in most of the allowed regions the right sneutrino VEVs take small values, \(v_R/\sqrt{2} \lesssim 500\) GeV, as shown in Fig. 32. These small values imply that \(\kappa \gtrsim 0.2\) to avoid tachyonic right sneutrinos. We show in Fig. 34 the value of the corresponding \(T_\lambda \), whereas in Fig. 35 the supergravity parameter \(A_\lambda \) is shown.

In Fig. 33 we show \(\tan \beta \) which, given the large value of \(\lambda \), can take smaller values than in \(S_2\). This region of the parameter space also favours light third generation squarks, as shown in Fig. 13 (see also Figs. 36 and 37).

Fig. 12
figure 12

Viable points of the parameter space for \(S_3\) in the \(\kappa -\lambda \) plane. The light-red (light-blue) colour represents cases where the SM-like Higgs is (is not) the lightest scalar. All light-red and light-blue points below the black dashed line fulfill the perturbativity condition up to 10 TeV of Eq. (26)

Fig. 13
figure 13

Viable points of the parameter space for \(S_3\) in the \(-T_{u_3}\) versus \(m_{{\widetilde{Q}}_{3L}}\) plane. The color code is the same as in Fig. 12

Fig. 14
figure 14

The singlet component \(\sum _i|Z^H_{h \widetilde{{\nu }}^{{\mathcal {R}}}_{iR}}|^2\) of the singlet-like scalars h versus their masses, for \(S_3\). We only show here viable points with scalar masses smaller than 1000 GeV. In the lower part we zoom in the low-mass region

As discussed in Sect. 3.1, the push-down effect (together with negative loop corrections) of a heavy singlet-like sector is more favourable to reproduce the SM-like Higgs mass. In Table 12 of Appendix C, we show the BP S3-R1 corresponding to the light-red region of Fig. 12. In this case, \(h_2\) is the SM-like Higgs, and the right sneutrinos are very mixed as expected from the large value of \(\lambda \). In Fig. 14, we show the singlet component of the singlet-like scalars. As for the other scans, scalars with large masses and with a very large doublet composition can also be present

Fig. 15
figure 15

Viable points of the parameter space for \(S_3\) in the \(-T_{\kappa }\) versus \(2\kappa {{\mathcal {M}}}\) plane. The color code is the same as in Fig. 12. In the upper right we zoom in the region with light-blue points

Although more difficult than in previous scans, we are also able to find in \(S_3\) solutions with light singlet-like scalars (light-blue region). However, as we can see in Fig. 14, no solutions with masses \(\lesssim m_{\text {Higgs}}/2\) are present. In Table 13, we show the BP S3-B1 corresponding to the blue region of Fig. 12. In Fig. 15, we show the correlation that is necessary to find these points, and in Figs. 38 and 39 we show the different values of \(-T_\kappa \) and \(A_\kappa = T_\kappa /\kappa \), respectively, en the \(\kappa -\lambda \) plane. The upper bound for \(A_\kappa \) is around 500 GeV in this case, but this is an artifact of our scan. If we had allowed in \(S_3\) values of \(T_\kappa \) up to 1 TeV, then the upper bound for \(A_\kappa \) would have been around 1 TeV.

In this scan, solutions with singlet-like states with masses close to 125 GeV are more rare. Only about 0.2% of the phenomenologically viable points found are of this type. We show in Table 14 the BP S3-2h1 as an example. As we can see, the right sneutrino \(h_1\) has the largest doublet composition of \(H^{{\mathcal {R}}}_{u}\) (1.42%) and \(H^{{\mathcal {R}}}_{d}\) (4.78%), but its mass if far away from 125 GeV. For this BP also the Majorana mass is small as for S2-2h1, \({\mathcal {M}}=50.3\) GeV, and there are three neutralinos dominantly right-handed neutrinos with masses of that order, 50.7, 51.8 and 64.8 GeV. As a consequence, the decay channel right sneutrino to two right-handed neutrinos opens for \(h_3\) and \(h_4\) (also for the SM-like Higgs \(h_5\)), but is not possible for \(h_1\). The latter can decay to right-handed neutrino plus light neutrino, but with a very small BR. In Table 15, we give another BP of this kind, S3-2h2, where now one of the singlet-like states, \(h_1\), contributes to the superposition of signals with the SM-like Higgs \(h_2\).

6 Conclusions

We performed a dedicated analysis of the parameter space of the \(\mu \nu \mathrm{SSM}\), in the light of the increasing data about the properties of the SM-like Higgs boson. For sampling the Higgs sector, we used a powerful likelihood data-driven method based on the algorithm MultiNest. The states of the Higgs sector crucial for our analysis are the two Higgs doublets and the three right sneutrinos, which are mixed among themselves. After determining the relevant parameters related to this sector (see Eq. 58), we performed scans to search for points compatible with the latest experimental data on Higgs physics. For constraining the predictions of our extended Higgs sector, we interfaced HiggsBounds with MultiNest, and to address whether a given Higgs scalar of the \(\mu \nu \mathrm{SSM}\) is in agreement with the signals observed by ATLAS and CMS we also interfaced HiggsSignals with MultiNest. In addition, we demanded the compatibility of the points with observables such as B and \(\mu \) decays, and discussed the values of muon \(g-2\) in different regions of the parameter space.

In this framework, we performed the three scans described in Table 1, which are determined by the range of \(\lambda \) couplings in the superpotential mixing Higgses and right sneutrinos, \(\sum _i \lambda _i {\hat{\nu }}^c_i {\hat{H}}_u {\hat{H}}_d\) (with \(\lambda _i=\lambda \)). In particular, we considered the three ranges \(\lambda \in \) [0.01, 0.2), [0.2, 0.5), and [0.5, 1.2). Perturbativity up to the GUT scale for \(\lambda \) is not imposed, and that is why we allow values of \(\lambda \) larger than 0.4. Neither we imposed perturbativity up to the GUT scale for \(\kappa \) couplings in the superpotential among right sneutrinos, \(\sum _{i,j,k} \kappa _{ijk} {\hat{\nu }}^c_i{\hat{\nu }}^c_j{\hat{\nu }}^c_k\) (with \(\kappa _{iii}=\kappa \) and vanishing otherwise), considering therefore the range \(\kappa \in \) [0.01, 2]. The results are summarized in Figs. 48, and 12 for the three scans. Clearly, we find viable solutions in almost the entire \(\kappa -\lambda \) plane with the exception of the scan \(S_3\) in Fig. 12, which is more constrained. This is due to the large values of \(\lambda \in \) [0.5, 1.2) that can give rise to tachyons originated in the mixing between the two Higgs doublets.

We have obtained therefore that the parameter space of the \(\mu \nu \mathrm{SSM}\) contains many viable solutions, including also many different phenomenological possibilities. For example, there are solutions where the SM-like Higgs is the lightest scalar (red and light-red points in the figures), but also solutions where right sneutrino-like states are lighter (blue and light-blue points). In the latter case, it is even possible to have these (singlet-like) scalars with masses \(\lesssim m_{\text {Higgs}}/2\). In addition, we also find solutions where several scalars are degenerated with masses close to 125 GeV, and can have their signals rates superimposed contributing to the resonance observed at 125 GeV.

Given these results, it is then important to study in detail the collider phenomenology of the solutions found. In particular, the impact of the new states, not only the right but also the left sneutrinos, and the neutralinos containing right-handed neutrinos. Novel signals associated to them might help to probe the \(\mu \nu \mathrm{SSM}\) at the LHC. These analyses will be carried out in a forthcoming publication [83].