Vacuum (meta-)stability in the $\mu\nu$SSM

We perform an analysis of the vacuum stability of the neutral scalar potential of the $\mu$-from-$\nu$ Supersymmetric Standard Model ($\mu\nu$SSM). As an example scenario, we discuss the alignment-without-decoupling limit of the $\mu\nu$SSM, for which the required conditions on the Higgs sector are derived. We demonstrate that in this limit large parts of the parameter space feature unphysical minima that are deeper than the electroweak minimum. In order to estimate the lifetime of the electroweak vacuum, we calculate the rates for the tunneling process into each unphysical minimum. We find that in many cases the resulting lifetime is longer than the age of the universe, such that the considered parameter region is not excluded. On the other hand, we also find parameter regions in which the electroweak vacuum is short-lived. We investigate how the different regions of stability are related to the presence of light right-handed sneutrinos. Accordingly, a vacuum stability analysis that accurately takes into account the possibility of long-lived metastable vacua is crucial for a reliable assessment of the phenomenological viability of the parameter space of the $\mu\nu$SSM and its resulting phenomenology at the (HL)-LHC.


Introduction
The Standard Model (SM) of particle physics successfully describes the fundamental interactions of the known matter fields at the present level of experimental accuracy. In particular, the discovery of a scalar particle with a mass of ≈ 125 GeV at the Large Hadron Collider (LHC) [1,2], that within the current experimental uncertainties [3][4][5] behaves as predicted by the realisation of the Brout-Englert-Higgs mechanism within the SM [6,7], has confirmed the existence of a non trivial vacuum structure of the universe. Despite the great success of the SM, there are also phenomena that cannot be accommodated, and furthermore the SM has conceptual shortcomings. In particular, the SM does not provide a sizable contribution to the relic abundance of cold dark matter, the neutrino oscillation data cannot be explained with strictly massless SM neutrinos, there is no mechanism to account for the observed matter-antimatter asymmetry of the universe, and there is no consistent formulation of the SM that includes gravitational interactions.
Another conceptual problem of the SM is related to the fact that its scalar sector, incorporating a single SU(2) doublet scalar field, just provides a minimal parameterization of electroweak symmetry breaking. The renormalizable scalar potential of the SM can be expressed in terms of only two parameters, e.g., the vacuum expectation value (vev), v ≈ 246 GeV, and the mass of the predicted Higgs boson, which within the SM is identified with the mass of the detected particle, m h 125 ≈ 125 GeV. This minimal construction is sufficient to generate the masses of the gauge bosons, the Higgs boson, the quarks and the charged leptons via spontaneous symmetry breaking. However, within the SM there is no dynamical origin of the Higgs potential, and the mass of the Higgs boson of the SM is not protected by symmetries against large corrections from physics at high scales. This gives rise to the so-called hierarchy problem [8,9] regarding the observed value of m h 125 and more generally the stability of the electroweak (EW) scale.
One of the most thoroughly studied frameworks of BSM physics is Supersymmetry (SUSY). Models based on SUSY provide a solution for the hierarchy problem, since the symmetry between bosonic and fermionic degrees of freedom leads to a cancellation between contributions with a quadratic dependence on new physics at higher energies [10][11][12], thus protecting the masses of the Higgs bosons of the model. Furthermore, for the observed gauge bosons and the states in the Higgs sector SUSY predicts new fermionic superpartners, as well as new bosonic superpartners for the three generations of quarks and leptons. As a result, SUSY extensions of the SM have a substantially richer matter sector, and in addition to the hierarchy problem many other shortcomings of the SM can be addressed as well. For instance, there are several possible candidates for cold dark matter.
A particularly well motivated SUSY model is the µ-from-ν Supersymmetric Standard Model (µνSSM) [13,14]. 1 Beyond the well-known appealing features of commonly studied SUSY models, in the µνSSM the tiny neutrino masses and their mixings can be accommodated via an EW seesaw mechanism, where it is required that the matter content is enlarged w.r.t. the SM by right-handed neutrinos [16][17][18][19][20]. Their superpartners, the "right-handed" 2 scalar neutrinos (sneutrinos), are gauge singlet scalar fields. If the right-handed sneutri-nos acquire vevs, the so-called µ-term of the Minimal Supersymmetric Standard Model (MSSM) can be generated effectively, in complete analogy to the Z 3 symmetric Next-to-MSSM (NMSSM). Consequently, the µνSSM also provides a solution to the so-called µproblem [21,22]. By construction, the µνSSM does not allow for a consistent assignment of conserved R parity charges, giving rise to the fact that there is no stable lightest SUSY particle. Compared to the (N)MSSM, the collider constraints from the LHC are therefore substantially weaker [23][24][25][26][27][28]. Regarding dark matter, it is possible to accommodate the measured relic abundance by means of decaying, but long-lived, gravitinos or axinos [29][30][31].
One of the key features of the µνSSM is the mixing of the two Higgs doublet fields H d and H u with the three left-handed and right-handed sneutrino fields ν iL and ν iR , leading to a total of eight neutral CP-even particles and seven neutral CP-odd scalar particles. 3 While the mixing of the left-handed sneutrinos with the Higgs doublet fields is suppressed by the smallness of lepton number breaking, the mixing of the right-handed sneutrinos with the Higgs fields can be sizable, while being in agreement with the current experimental limits [27,32]. The phenomenological consequences of the latter possibility have been thoroughly studied, paying special attention to the prediction for the mass of the SM-like Higgs boson [32][33][34]. In addition to the more complicated scalar particle spectrum, the enlarged Higgs sector compared to the SM and the (N)MSSM also results in a substantially more complicated scalar potential. While it is always possible to ensure the existence of a phenomenologically viable EW vacuum by conveniently choosing the vevs of the neutral scalar fields as input parameters (as will be discussed in Sect. 2), the scalar potential of the µνSSM in general contains further unphysical local minima. If one or more of the unphysical minima are deeper than the EW minimum, the question arises whether the EW vacuum can become short-lived in comparison to the age of the universe, and if this is the case the corresponding parameter point should be rejected.
Up to now no detailed analysis of constraints on the parameter space of the µνSSM arising from possible instabilities of the EW vacuum has been carried out in the literature. To the best of our knowledge, only Ref. [16] took into account such constraints, demanding, however, that the EW minimum is the global minimum of the potential, i.e. the possibility of a long-lived metastable EW vacuum was not taken into account. Apart from that, heuristic constraints that have been obtained for the presence of color-breaking unphysical minima in the MSSM could be considered [35]. Besides their limited reliability, constraints of this kind do not capture the possibility that charge-and color-conserving minima can nevertheless be unphysical. Accordingly, we study in this paper the impact of genuine effects of the µνSSM on the stability of the EW vacuum.
In this paper, we perform a detailed analysis of the neutral scalar potential of the µνSSM regarding EW vacuum instabilities. In a first step, we determine the unphysical minima of the tree-level neutral scalar potential by making use of an homotopy continuation method. In a second step, we estimate the rates for the decays of the EW vacuum into all unphysical minima below the EW minimum using a semi-classical approximation. By subsequently comparing these rates to the age of the universe, we categorize parameter points into ones with a short-lived (unstable), long-lived (metastable) and an absolutely stable EW vacuum, and we discuss the impact of the resulting constraints on the parameter space of the model.
As an example scenario, in our numerical analysis we apply the investigation of vacuum stability constraints to the alignment-without-decoupling limit of the µνSSM. It turns out that in this limit the EW minimum is not the global minimum of the potential. However, in large parts of the analyzed parameter regions the decay rates into the unphysical minima are so small that the EW vacuum is long-lived compared to the age of the universe, and those parameter regions with a metastable vacuum are therefore phenomenologically viable.
The paper is organized as follows: In Sect. 2 we introduce the Higgs sector of the µνSSM, followed by a discussion of the alignment without decoupling limit in Sect. 3. In Sect. 4 we describe the determination of the minima of the neutral scalar potential and the procedure that is applied for assessing the stability of the EW vacuum. In Sect. 5 we apply this approach to example scenarios in the alignment limit. Finally, we summarize our results and give a brief outlook on possible avenues for future research in Sect. 6. In the appendix we briefly describe the public code munuSSM that can be used for performing the vacuum stability test in the µνSSM parameter space. We furthermore provide explicit expressions for the coefficients of the scalar potential relative to the EW vacuum for some special cases.

The neutral scalar sector of the µνSSM
The neutral scalar potential of the µνSSM consists of contributions from the superpotential W and the soft SUSY-breaking Lagrangian L soft . The part of W that contains the neutral superfields is given by [13,14] are the Higgs doublet superfields,L i are the leftchiral lepton superfields with the left-chiral neutrino superfieldsν i in the upper component, andν c are the right-chiral neutrino superfields. i, j = 1, 2, 3 are the family indices, and a, b = 1, 2 are the weak isospin indices ( 12 = 1). The portal couplings λ i give rise to the mixing between the right-handed sneutrinos and the Higgs doublet fields, and also the µterm is generated after electroweak symmetry breaking (EWSB) (µ = λ i v iR / √ 2). The selfcouplings κ ijk generate masses for the right-handed (s)neutrinos after EWSB, such that they determine the scale of the seesaw mechanism. Since the seesaw scale is in this way related to the EW scale, for the neutrino Yukawa couplings one has to demand Y ν ij Y e ≈ 10 −6 , with Y e the electron Yukawa coupling, in order to obtain left-handed neutrino masses at the sub-eV level.
The part of L soft that contains the neutral scalar fields can be written as [13,14] − where H d , H u , ν iL and ν iR denote the scalar components of the superfieldsĤ d ,Ĥ u ,ν i and ν c i , respectively. As usual, soft mass parameters that are non-diagonal in field space are neglected. Thus, terms of the form ((m 2 H d L ) i H a * d L a iL + h.c.) are absent, and the matrices m 2 L and m 2 ν are assumed to be diagonal. For the trilinear scalar couplings, we will make use of the decomposition which is motivated in models of supergravity with diagonal Kähler metric [36].
The soft terms shown above, in combination with the D-and F -terms derived from the superpotential, give rise to the neutral scalar potential The required existence of a physical EW minimum can be made explicit by using the decomposition where the superscripts R and I denote CP-even and CP-odd components of each scalar field, respectively. The numerical prefactor ensures that the kinetic terms of the real fields H R d,u , H I d,u , ν R iR,L and ν I iR,L are canonically normalized, provided that the complex scalar fields are canonically normalized (see the discussion in Ref. [37]). The EW minimum is then defined by the vevs In order to simplify the notation, the CP-even and the CP-odd components of the scalar fields will be collectively denoted by , respectively. For practical reasons, we focus only on the CP-even part of the scalar potential to reduce the number of field dimensions, i.e. in our analysis we assume that σ i = 0, and the parameters v d , v u , v iR , v iL are assumed to be real. 4 The vevs defined in Eq. (7) will be used as input parameters. 5 The soft mass parameters m 2 H d , m 2 Hu , (m 2 L ) ii and (m 2 ν ) ii are then given as a function of the vevs demanding that the minimization conditions Thus, our analysis does not take into account the possibility of further CP-or charge-breaking minima. 5 Regardless of the prefactors 1/ √ 2 in Eqs. (8), we will refer to the parameters v u , v d , v iL and v iR as vevs.
are fulfilled. Moreover, for phenomenological reasons the Hessian matrix of V as a function of φ i and σ i has to be positive definite in the EW vacuum in order to avoid tachyonic states. It should be noted that in principle a parameter point with a tachyonic state at treelevel could still be physical, since sufficiently large radiative corrections to the scalar masses might give rise to a change of the curvature in a particular field direction, thus leading to the presence of only non-tachyonic physical states. In our analysis, we do not take into account radiative corrections to the scalar potential, and accordingly we discard parameter points with tachyonic states at tree-level. In general, the inclusion of loop corrections for the analysis of the vacuum stability would demand a conceptually different approach than the one pursued here, because we make use of the fact that the tree-level potential only contains polynomial contributions (see also Sect. 4). Here it should be kept in mind that the question whether the usual approach of taking into account higher-order corrections in the form of an effective Coleman-Weinberg potential actually improves the calculation of the lifetime of the EW minimum is still an open issue. This is related to the fact that the effective potential formulation does not capture momentum-dependent contributions, which were demonstrated to be relevant for the calculation of the decay rates, leading to an inconsistent perturbative expansion when truncated at a certain order ofh [38,39] (see also Ref. [37] for a discussion). In addition to the absence of tachyonic states, it is required for a phenomenologically viable parameter point that the vevs of the fields charged under the EW symmetry have to fulfill the condition v GeV, such that the observed values of the gauge boson masses are properly reproduced. It is convenient to define tan β = v u /v d in order to make a connection to the MSSM. Using the above relations, v and tan β can be used as input parameters instead of v d and v u . From Eq. (12) one can deduce, using T ν ij = A ν ij Y ν ij , that a solution to this equation The scale of v iR , on the other hand, is given by the SUSY-breaking scale via the parameters (m 2 ν ) ij , T κ ijk and T λ i in Eq. (11). The CP-even part of the scalar potential V at lowest order is a polynomial in eight (field) dimensions with degree four. Consequently, it can have numerous coexisting local minima besides the EW minimum. 6 Some of these minima are physically redundant, because the potential contains accidental discrete symmetries yielding degenerate stationary points. For instance, one can see that if φ = (v d , v u , v iR , v iL ) is a solution to Eqs. (9)-(12), then there is always a second physically identical solution at φ = (−v d , −v u , v iR , −v iL ). However, there can also be solutions of the minimization conditions that are potentially dangerous for the stability of the EW vacuum. Any field configuration φ i that solves Eqs. (9)- (12), that has a positive definite Hessian matrix, and for which the value of the potential in the minimum is smaller than the one of the EW minimum, constitutes a deeper local minimum of the potential. The existence of such a deeper local minimum gives rise to the possibility of a rapid EW vacuum decay. We will refer to such minima as unphysical minima, since they do not fulfill v ≈ 246 GeV. Furthermore, due to the fact that the values of v iL are different in the unphysical minima, the predicted neutrino masses are also modified compared to the prediction based on the EW minimum.
Consequently, for a phenomenologically viable parameter point one has to verify that the EW minimum is either the global minimum of the potential (but there are constraints that apply even in this case, see below), or that, if deeper unphysical minima exist, the lifetime of the (metastable) EW vacuum is large in comparison to the age of the universe. As discussed in Refs. [40,41], taking into account the thermal history of models with multiple scalar fields, one finds that a parameter point featuring a global EW minimum at zero temperature may still be unphysical. This happens if the transition to the EW vacuum would have occurred via a first-order phase transition, but the transition probabilities turn out to be never large enough to allow the onset of the bubble nucleation of true EW vacuum bubbles in the early universe. In this case, the universe would be trapped in one of the unphysical vacua in the limit of zero temperature despite the fact that the EW minimum is deeper than the unphysical one. In our analysis in this paper we restrict ourselves to the case of zero temperature. One of the important findings of our analysis will be that the rich vacuum structure of the µνSSM neutral scalar potential gives rise to the fact that several unphysical vacua can be present simultaneously, such that there are several possibilities for vacuum trapping of the universe in an unphysical field configuration. The analysis of the stability of the EW vacuum at T = 0, which is the focus of the present paper, is carried out under the assumption that the thermal history of the universe has been such that the universe is actually in the EW minimum at T = 0. This analysis that we will perform can be used to determine constraints on the parameter space of the µνSSM essentially independently of the question which additional constraints would arise from a dedicated analysis of the thermal history of the universe. We leave the latter kind of analysis, which in particular takes into account constraints arising from the above-mentioned effect of vacuum trapping, for future work.
Since the multi-dimensional parameter space of the µνSSM cannot be covered entirely, we focus our numerical discussion on the alignment without decoupling limit, which is theoretically and experimentally well motivated in view of the fact that the properties of the Higgs boson at 125 GeV are SM-like. The conditions that apply for the parameter space of the model in this limit and the resulting phenomenological features will be introduced in the following section.

Alignment without decoupling in the µνSSM
We start by investigating the alignment without decoupling conditions for the µνSSM (for partial results see also Ref. [42]). If these conditions are respected at least approximately, it is possible to obtain a particle state at 125 GeV with properties resembling the ones of the SM Higgs boson without relying on a decoupling of the remaining doublet-like scalars from the EW scale. In the following discussion we will make use of the relations The first two inequalities arise from the requirement to suppress lepton-number breaking. The third expression, i.e., treating κ ijk as diagonal, is a reasonable simplification since the superpotential can always be written in a basis in which the right-handed sneutrino selfcouplings κ ijk are diagonal. The diagonal structure of κ ijk is also preserved by the RGE evolution of the couplings [32]. The general strategy for determining the alignment conditions follows the approach of Ref. [43], therein applied to the NMSSM. First, one rotates the scalar squared mass matrix from the interaction basis into the so-called Higgs basis, in which only one scalar field obtains an EWSB vev. This transformation can be expressed in terms of a unitary transformation matrix U HB that will be specified below. One can then identify relations between the free parameters of the model for which the non-diagonal mass matrix elements between the field that obtains the EWSB vev and the other scalar fields vanish. The field that is aligned with the EWSB vev then couples to the gauge bosons and fermions of the SM in the same way as the Higgs boson that is predicted by the SM [44]. In this limit, the CP-even Higgs sector of the µνSSM consists of a SM-like Higgs boson, H SM , a second doublet-like heavy Higgs boson, H NSM , three singlet-like Higgs bosons, H 1,2,3R , and three decoupled left-chiral sneutrino states, H 1,2,3L . 7 We will make use of a similar unitary transformation in terms of the matrix U AB for the CP-odd part of the Higgs sector. This will allow us to define a tree-level mass parameter M A in analogy to the (N)MSSM. This mass parameter is roughly equal to the tree-level mass of the massive doublet-like CP-odd Higgs boson A of the µνSSM in the alignment limit.
The neutral scalar mass terms in the interaction basis are given by where the squared mass matrices M s and M p are determined by the curvature of the scalar potential V in the direction of the fields φ and σ in the EW minimum. The tree-level entries of these matrices can be found in Ref. [32]. In the following it is implied that in the mass matrices the soft mass parameters m 2 H d , m 2 Hu , (m 2 ν R ) ii and (m 2 ν L ) ii were replaced by the vevs v d , v u , v iR , and v iL , making use of Eqs. (9)- (12). We transform L into the Higgs basis via Here the scalar fields in the Higgs basis are expressed in terms of the field vector where G 0 is the neutral Goldstone boson and A is the doublet-like particle state. For the elements that are given with two different signs, the upper and lower signs refer to U HB and U AB , respectively. The different signs arise from the fact that the chiral superfields H d and ν iL have the opposite hypercharge as H u . Thus, a field redefinition via a global SU(2) transformation and a complex conjugation of the doublet fields H d and ν iL is applied in order to use a basis in which all fields with nonzero hypercharge have the same value of the hypercharge. For the imaginary components of the scalar fields, σ i , (see the definition as specified in Eq. (7) and the text below) the complex conjugation translates into the differences of the signs of the elements of U AB and U HB . 8 Inverting the transformations defined in Eq. (15), one can replace φ i = (U HB ) −1 ij H j and σ i = (U AB ) −1 ij A j in Eq. (14) to obtain the mass matrices in the Higgs basis, such that As already mentioned before, the alignment conditions are precisely defined by demanding that the state H SM is aligned with the vacuum expectation value and has vanishing mixing with the other states, i.e., M HB 1i = 0, where i = 1. The resulting conditions on the model parameters will be discussed in the following. In order to make a distinction between the fields in the Higgs basis H i and the (loop-corrected) mass eigenstates, we will use the notation h i for the latter, where the index i reflects the mass hierarchy of the scalars.

First alignment condition
The first alignment condition results from requiring that the mixing between the two Higgs doublet states H d and H u is such that (ignoring for now a possible singlet admixture) one mass eigenstate h i ≈ H SM couples to the SM gauge bosons with a coupling that is approximately equal to the one of the SM Higgs boson, and the other doublet-like state h j ≈ H NSM has a vanishing coupling to gauge bosons as it does not acquire an EWSB vev. The alignment condition arises from the requirement where GeV, summation over i = 1, 2, 3 is implied, and we use the short-hand notation c x ≡ cos x and s x ≡ sin x. Using also that it is convenient to write (with t x = tan x) The expression shown in Eq. (18) is the squared tree-level mass of the SM-like Higgs boson in the exact alignment limit. It is known that in SUSY models large radiative corrections from the stop sector can be present [45]. As discussed in Ref. [43], to a very good approximation these radiative corrections only enter in M HB 11 , but not in the off-diagonal entry M HB 12 . Thus, in order to account for those radiative corrections it is sufficient to substitute M HB 11 → m 2 h SM ≈ (125 GeV) 2 in Eq. (19), and the condition shown in Eq. (17) can be written as This is what we will refer to in the following as the first alignment condition, which is a straightforward generalization of the NMSSM condition [43].

Second, third and fourth alignment conditions
In the NMSSM, there is a second alignment condition that ensures that there is no mixing between H SM and the gauge singlet field. In the µνSSM, there are three gauge singlet right-handed sneutrino fields. Each of the them can potentially mix with the SM-like Higgs boson. 9 In order to avoid a singlet admixture in the state H SM , one has to fulfill three more conditions, i.e., Thus, alignment without decoupling is obtained if these three conditions are fulfilled together with the condition of Eq. (20). Inserting the mass matrix elements, one finds the conditions The corresponding NMSSM condition was derived in Ref. [46], and it was already generalized to the µνSSM in Ref. [42]. However, a phenomenological exploration of the µνSSM parameter space in combination with the first alignment condition shown in Eq. (20) has not been carried out yet. One can also obtain a relation in closer analogy to the NMSSM by realizing that, if all three conditions shown above are fulfilled, then also the sum vanishes, such that This expression allows one to replace the terms ∼ A λ i = T λ i /λ i by the tree-level mass parameter of the doublet-like CP-odd Higgs boson given by where the summation over i = 1, 2, 3 is implied and, as mentioned in the beginning, terms ∼ v iL are not written out. This leads to the condition where again the summation over i = 1, 2, 3 is implied. The first term in Eq. (25) exactly coincides with the NMSSM result [43]. The second term has become considerably more complicated due to the presence of more than one gauge singlet field. However, assuming that two of the three singlets are decoupled, e.g., λ 2,3 → 0, the NMSSM formula is recovered when substituting v 2 1R = 2µ 2 /λ 2 1 in the surviving term. It should be noted, however, that in the µνSSM the condition shown in Eq. (25) is a necessary, but not a sufficient criterion for achieving the alignment limit. In other words, Eq. (25) must be true in the alignment limit, i.e., when the conditions shown in Eq. (22) are fulfilled separately, but not all parameter configurations that respect Eq. (25) correspond to the alignment limit. Finally, by noting that in order to fulfill the first alignment condition one has to require rather large values of λ 2 1 + λ 2 2 + λ 2 3 0.6 2 , a possibility to avoid Landau poles below the GUT scale is to assume κ i λ i . In this case the second term in Eq. (25) is small, and it is sufficient to use as a condition M 2 A s 2 2β /(4µ 2 ) ≈ 1, just as in the NMSSM [43].

Vacuum stability
The main goal of this work is to investigate the stability of the EW vacuum of the µνSSM at temperature T = 0. This analysis can be divided into two different tasks. First, one has to check whether there are local minima present in the scalar potential that are deeper than the EW minimum. If there is no such unphysical minimum, the EW minimum is the global one, and the EW vacuum is stable (see, however, the discussion on vacuum trapping in Sect. 2). If one or more deeper unphysical minima are present, one has to investigate in a second step whether the EW vacuum can be regarded as sufficiently long-lived or whether it would rapidly decay into one of the unphysical vacua. In Sect. 4.1 we will briefly introduce our approach for determining the different minima of the neutral scalar potential. Subsequently, we give details on the calculation of the lifetime of metastable EW vacua in Sect. 4.2.
The computations that are carried out for testing vacuum stability as described in the following are implemented in the public code munuSSM v.1.1.0 [32][33][34]. More details on the implementation and simple user instructions are given in App. A.

Finding potentially dangerous unphysical vacua
Due to the large hierarchy between the doublet and singlet vevs, v d , v u , v iR ≈ v, and the left-handed sneutrino vevs, v iL ≈ Y ν ij v, it is very challenging to find the local minima of the potential V with the help of usual minimization algorithms. In particular, gradient based methods poorly converge due to the nearly flat directions in the potential originating from the large hierarchy of parameters. Moreover, such algorithms crucially depend on the initial conditions (such as the starting values for the fields), and it is very challenging to find all (or at least most) of the various different minima of the µνSSM scalar potential in this way.
We therefore use an approach that is based on directly solving the polynomial minimization conditions shown in Eqs. (9)- (12). The solutions yield the stationary points of the potential. In a second step we determine the stationary points corresponding to a local minimum of the potential by calculating the eigenvalues of the Hessian matrix of the potential for each solution that was found. For all local minima the value of the scalar potential is calculated. In this way the minima that are deeper than the EW minimum, and therefore potentially dangerous for the stability of the EW vacuum, are determined. For each minimum of this kind a calculation of the transition probability is carried out as described in Sect. 4.2.
The solutions to the stationary conditions were obtained using the public code HOM4PS-2.0 [47], which is an implementation of the polyhedral homotopy continuation method [48][49][50]. In general, the code is efficient in finding all existing solutions. In our analysis it happened only in rare cases that a solution was missed by HOM4PS-2.0 (see also the discussion in Sect. 5 below), which could be tested by computing the same parameter point several times and comparing the number of solutions that were found. This can have different reasons. First, it can happen that the code converges twice to the same solution from two different starting points on the unit circle in the complex plane, for instance, if at intermediate stages of the algorithm two different solutions only differ by the tiny values of the left-handed sneutrino fields ν R iL . Secondly, one has to define a numerical uncertainty for the imaginary parts of the solutions to the stationary conditions in order to identify which of the solutions correspond to a real solution, since they are the ones we are interested in here. If this uncertainty is set to too small values, i.e., below the numerical precision of the algorithm, a real solution is misidentified as a complex one and dismissed erroneously (see also Ref. [37]).
We finally remark that it is important to remove the redundant degrees of freedom related to the gauge symmetries. Otherwise, any stationary point would be transformed into flat directions along the redundant degrees of freedom, and the method introduced above would fail. In our analysis, this problem is absent since we restrict the analysis to the case in which only the CP conserving real parts of the neutral scalar fields are allowed to develop a vev, such that no degeneracies related to the gauge symmetries are present.

Lifetime of the EW vacuum in the presence of deeper minima
In case one or more deeper unphysical minima are present, the EW vacuum could decay via quantum tunneling effects into an unphysical vacuum. The probability for this vacuum decay mainly depends on the euclidean bounce action S E for the classical path in field space from the EW into the unphysical minimum. The computation of this bounce action is a complicated task, in particular for cases where the number of field dimensions is quite large. 10 Due to the fact that in the µνSSM there are (at least) eight real fields that have to be considered, and since there are large hierarchies between the field values, one has to find a numerically efficient, fast and reliable way of calculating the tunneling probabilities. We follow here the approach of Ref. [37] (see also Refs. [52][53][54]), in which S E is obtained by making use of a semi-analytical approximation, and a simple criterion is applied in order to determine whether a certain parameter point features a short-lived, long-lived, or stable EW vacuum. We briefly describe this approach in the following.
At tree level the CP-even scalar Potential V (φ i ) as given in Eq. (3) can be brought into the form [37] V where the fields were redefined via a shift to the EW vacuum configuration, such that φ i = φ EW i + ϕ i , and the field vector ϕ i is then expressed in terms of its norm ϕ = (ϕ i ϕ i ) 1/2 and the unit vectorφ i , such that ϕ i = ϕφ i . Since by construction φ EW solves the stationary conditions of the potential there is no linear term in Eq. (26). In addition, the field independent terms were discarded in Eq. (26), so that V EW (0,φ) = 0.
After having identified the unphysical minima, one can investigate V EW along the direction from the EW minimum towards each of the unphysical minima using a straight-path approximation [37]. For a fixed unit vectorφ i ≡φ i,d pointing towards a deeper unphysical has two local minima, the EW minimum at ϕ = 0 and the unphysical minimum at ϕ d ≡ ∆φ, which are separated by a potential barrier originating from the trilinear A-term (see also Fig. 3 in the numerical discussion below). The unphysical minimum at φ i,d is therefore the global minimum of V EW (ϕ), and V EW (ϕ d ) is equal to the potential difference ∆V EW between the unphysical minimum and the EW minimum.
For the one-dimensional potential V EW (ϕ), which has been obtained as described above making use of a straight-path approximation, it is possible to calculate the bounce action for a transition from the EW into the deeper unphysical minimum. This calculation was carried out in Ref. [55] for potentials of the shape as defined in Eq. (26), and it was shown that a very good numerical fit to the result for S E is given by In our numerical analysis we compared the results for S E based on this semi-classical approximation with the ones obtained with the public code cosmoTransitions [56], and found deviations between both calculations below 1% for all parameter points considered. 11 If the 10 See, for instance, Ref. [51] for more details and different numerical approaches applied in the literature to compute the bounce action. 11 For this comparison, we implemented the one-dimensional scalar potential V EW (ϕ) into cosmoTransitions. Thus, the calculation does not benefit from the path deformation method as provided by the code for the case of multiple field dimensions. Given the full CP-even scalar potential V of the µνSSM, this method could in principle be applied in order to obtain a more precise result. However, we observed that the path deformation algorithm as implemented in cosmoTransitions did not converge successfully. We attribute this to the large number of fields and the large hierarchy of the field values in the µνSSM. kinetic terms are canonically normalized, the decay rate Γ for the quantum tunneling into the unphysical minimum under consideration per space volume V is then given by [57,58] where K is a sub-leading prefactor that incorporates higher-order corrections. For practical purposes, it appears to be sufficient to estimate K only based on dimensional arguments (see the discussion in Ref. [37]), such that K ≈ M 4 . Here M represents a typical mass scale of the model. In order to estimate whether the EW vacuum of a parameter space point of the µνSSM is sufficiently long-lived, or whether it would rapidly decay into one of the unphysical vacua, we followed the criteria formulated in Ref. [37]. A parameter point is considered to feature a long-lived EW vacuum if the values for S E for all possible transitions into unphysical vacua are larger than 440. On the contrary, if for one of the possible transitions we find S E < 390, the EW vacuum is considered to be short-lived compared to the lifetime of the universe, and the corresponding parameter point should be rejected. The region in between, 390 < S E < 440, arises from the uncertainty stemming from the unknown prefactor K. In our numerical analysis we will separately display parameter points for which S E falls into this region where in our approach no clear distinction between a short-lived and a sufficiently long-lived EW vacuum is possible. Before starting the numerical discussion in Sect. 5, we point out that in the µνSSM more than one transition into different unphysical minima with S E < 440 can be present for a single parameter point. Here it should be noted that the lifetime of the EW vacuum is practically unaffected by the number of deeper unphysical minima into which the EW vacuum could decay, because the numerically dominant exponential term e −S E is determined based on each possible vacuum decay individually. Taking into account that the prefactor K in our approach is estimated only based on dimensional arguments, it is therefore sufficient to equate the total decay rate of the EW vacuum with the maximum of the individual decay rates belonging to each possible transition. Consequently, our above considerations remain valid also in the presence of several unphysical vacua.

Numerical analysis
As already mentioned in Sect. 1, in our analysis of the possible impact of EW vacuum stability constraints on the parameter space of the µνSSM we will focus on the alignment without decoupling region. In the µνSSM an important difference to the corresponding limit in the NMSSM arises from the fact that there are three gauge singlet scalar fields present in the µνSSM. Each of these singlet sneutrino fields is coupled to the Higgs doublet fields via a portal coupling λ i , whereas in the NMSSM there is only one such coupling λ. Consequently, while in the NMSSM the first alignment condition fixes the value λ as a function of tan β, in the µνSSM the corresponding condition as given in Eq. (20) can be fulfilled for different individual values of λ i , while only their squared sum is given as a function of tan β. In order to be as general as possible, we decided to choose different values for λ i in a hierarchical structure, such that one right-handed sneutrino is rather strongly coupled to the Higgs doublet fields, another one is moderately coupled, and the third one is largely decoupled.  The precise values of the free parameters are summarized in Tab. 1, and we will briefly discuss the choice of their values in the following.
We have chosen tan β = 2, such that according to Eq. (20) the λ i have to fulfill the constraint λ 2 1 + λ 2 2 + λ 2 3 ≈ 0.65. Taking into account the above, we use values of λ 1 = 0.54, λ 2 = 0.36 and λ 3 = 0.005. These values yield a relatively large enhancement of the tree level mass of the SM-like Higgs boson. Thus, no large radiative corrections are required to achieve a physical mass of ≈ 125 GeV, and the values for the stop mass parameters, m Q, t = M S = 2.5 TeV, the soft trilinear stop coupling, A t = 500 GeV, and the gluino mass parameter, M 3 = 2 TeV, are set accordingly. 12 In order to fulfill the remaining alignment conditions as given in Eq. (22), we set the values of A λ i and κ i as given in Tab. 1. The parameters v iR and A κ i were varied in a range that yield right-handed sneutrino masses in the vicinity of 125 GeV, since this is the mass scale at which the sneutrinos can be mixed sizably with the SM-like Higgs boson. The corresponding parameter space is also of phenomenological interest as it can be probed experimentally by the LHC and possible future colliders [32]. As we will discuss below, the parameters v iR and A κ i are particularly important for the presence of a potential barrier between the EW minimum and deeper unphysical minima. Therefore, also the decay rates of the EW vacuum strongly depend on these parameters. Combining the impact of the sneutrino mixing with the SM-like Higgs boson on the one hand, and the considerations regarding the EW vacuum stability as discussed below on the other hand, there is an interesting complementary between experimental and theoretical constraints on the parameter space under investigation.
While we choose the same values of the parameters v iR , A ν i and A κ i for the three gen- we set the values of κ i slightly different. These differences are commonly included in studies of the µνSSM, since they ensure that there are not two right-handed sneutrinos with almost exactly the same masses in case the values for λ i are chosen to be equal, due to an accidental symmetry in the scalar potential (see Refs. [23,32] for details). Even though the λ i are not uniform here, we still follow this approach in order to increase the mass differences between the three families of right-handed sneutrinos slightly more. This leads to the fact that also the values of A λ i are not uniform, as their values are fixed by the second, third and fourth alignment conditions shown in Eq. (22).
The remaining parameters are mostly related to the EW seesaw mechanism and not directly relevant for the discussion here. In particular, for v iL = 10 −4 GeV and Y ν ij = δ ij Y ν i = 1.7 · 10 −7 one finds left-handed neutrino masses at the sub-eV level. Note, however, that we did not include a precise fit to the neutrino oscillation data in terms of squared mass differences and mixing angles, as such details have no impact on our discussion. Finally, the gaugino mass parameters are set to M 1 = 200 GeV and M 2 = 400 GeV, such that they are of roughly the same order as µ ≈ 500 GeV and the left-handed sneutrino masses m ν iL ≈ 560 GeV, where m ν iL are controlled by the values of A ν i = −1.5 TeV. Since the left-handed sneutrinos are substantially heavier than the right-handed sneutrinos and the SM-like Higgs boson, and since they are practically not mixed with each other, we do not discuss the left-handed sneutrinos any further in the following.
As shown in Tab. 1, we vary the two parameters v R ≡ v iR and A κ ≡ A κ i . We start our numerical investigation by varying each of the two parameters individually and analyze how the EW vacuum stability depends on them. At the end of this section, we will summarize our numerical discussion by showing the results in the two-dimensional parameter plane v R -A κ arising from a grid scan over both parameters. In the first scan, we set v R = 800 GeV and vary A κ . In the second scan, we set A κ = −880 GeV while varying v R in the given range. For each parameter point, we first calculated the radiatively corrected Higgs boson spectrum. Afterwards, we checked the point against constraints from BSM Higgs boson searches and the signal rate measurements of the Higgs boson at ≈ 125 GeV. For the parameter points that passed the constraints, we determined the unphysical minima and finally calculated the bounce actions for the available transitions from the EW vacuum into all potentially dangerous unphysical vacua. The whole analysis chain can be performed with the public code munuSSM v. 1 [47]. For the investigation of the stability of the EW vacuum, we have extended the code munuSSM by incorporating the new subpackage vacuumStability. 13 In Fig. 1 we show the radiatively corrected Higgs boson masses m h 1,...4 , corresponding to a SM-like Higgs boson at ≈ 125 GeV (≡ h 125 ) and the three CP-even right-handed sneutrino states with masses not far below or above 125 GeV. The colors of the points indicate the gauge singlet component of each Higgs boson, which we define as where R H is the loop corrected mixing matrix that transforms the CP-even fields φ i from the interaction basis into the fields h i in the mass eigenstate basis. 14 One can see that, as expected from the discussion in Sect. 3.2, h 125 has only a small singlet component over almost the entire scan range, despite the fact that the right-handed sneutrinos are close in mass. The only exceptions are the points for which h 125 is practically degenerate with the lightest right-handed sneutrino and where a "level-crossing" occurs, i.e. h 125 changes its character from being the next-to-lightest state, h 2 = h 125 , to the lightest state, h 1 = h 125 . In the region where the level-crossing takes place the properties of a SM-like Higgs boson are in fact shared between the states h 1 and h 2 , and both Higgs bosons contribute to the signal rate measurements at the LHC for the particle state at ≈ 125 GeV. Using HiggsSignals we applied a χ 2 -test to each parameter point shown in Fig. 1 in order to compare the predicted signal rates of h 125 with the experimental measurements. 15 We discuss the results of this test in terms of ∆χ 2 = χ 2 µνSSM − χ 2 SM , where χ 2 µνSSM is the fit result of the parameter points and χ 2 SM = 84.4 is the fit result assuming a SM Higgs boson at 125 GeV. For the points of the scan over A κ (left plot of Fig. 1), we find a good compatibility with the measured properties of h 125 , resulting in ∆χ 2 values of ∆χ 2 ≈ −1 to 3, except for the points for which m h 1 ≈ m h 2 ≈ 125 GeV, where ∆χ 2 ≈ 4 to 6. Here it should be noted that A κ does not appear in the alignment conditions. Accordingly, the variation of A κ does not have a sizable impact on the mixing of h 125 with the right-handed sneutrino states, except for the degenerate region with m h 1 ≈ m h 2 , in which the assumptions on which the alignment condition rely are not fulfilled. Consequently, no sizable impact on the properties of h 125 is expected from the variation of A κ in the analyzed parameter space outside of the degenerate region.
The situation is different when v R is varied (right plot of Fig. 1), because it gives rise to a variation of the µ parameter that appears in Eq. (22). The constraint ∆χ 2 < 6 was used on the scan range, except for the points for which h 1 is degenerate with h 2 , where we find values of ∆χ 2 of up to ∆χ 2 ≈ 12. Since we are mainly interested in the impact of the constraints on vacuum stability, and since the mass-degenerate points do not show any particular feature regarding the constraints from vacuum stability compared to the other points, we do not discard points with larger values of ∆χ 2 in the regions with m h 1 ≈ m h 2 ≈ 125 GeV (we note that in this region the predicted signal rates, to which both h 1 and h 2 contribute, show sizable deviations from the measured values). For the upper end of the displayed v R range ∆χ 2 increases as a consequence of the decrease in the predicted value for the mass of the SM-like Higgs boson. For values of v R 830 GeV we find m h 125 122 GeV, which is outside of the allowed range that we have employed based on a theoretical uncertainty of the Higgsmass prediction of 3 GeV. The lower end of the scan range of v R arises from constraints from direct searches, as described below.
Regarding the HiggsBounds test, which confronts each parameter point with the available constraints from BSM Higgs-boson searches at colliders, all points of the scan over A κ pass the experimental constraints. This is related to the fact that the particles corresponding to the right-handed sneutrinos are very singlet-like and have strongly suppressed couplings to the SM particles. 16 Only for the smallest values of A κ one can see an increase of the doublet component of the heaviest right-handed sneutrino (h 4 ), such that the channel pp → h 4 → ZZ * becomes important, which was searched for by CMS in the two-lepton and four-lepton final states and including width effects [77]. However, the predicted cross sections remain roughly 20% below the experimental limits even in this range of A κ . The same experimental channel is also the most sensitive one in the scan over v R , yielding an exclusion for values of v R < 770 GeV. Accordingly, we use this value as the lower limit of the scan over v R .
In Fig. 2 we show the results of the analysis of vacuum stability for the parameter space under investigation. On the left-hand side the points from the A κ -scan are shown, and on the right hand side we show the points from the scan over v R . In the upper row we show the values of the potential V for the EW minimum and for all unphysical minima with potential values below the EW minimum. The points corresponding to the EW vacuum are shown in gray. For the unphysical minima the colors of the points indicate the values of the euclidean bounce action S E for the tunneling from the EW minimum into the unphysical minima. Parameter points for which we find that the EW vacuum is short-lived compared to the lifetime of the universe, corresponding to S E < 390, are furthermore highlighted in red. The regions with 390 < S E < 440, for which in our approach no clear distinction between a short-lived and a sufficiently long-lived vacuum can be drawn, are highlighted in orange (see the discussion in Sect. 4.2).
In both scans it can be observed that the EW minimum is not the global minimum of the potential for any of the points. In the A κ -scan we find a total number of 12 deeper unphysical minima, where a few of those have the same depth V due to the accidental symmetries mentioned in Sect. 2, such that only eight lines are visible in the upper left plot of Fig. 2. The situation is similar for the upper right plot, with the exception that for this scan the lines of two unphysical minima merge for small values of v R . The different unphysical minima below the EW one can be more easily distinguished from each other in the plots displayed in the lower row where the field space distance between the unphysical minimum and the EW minimum, ∆φ, is shown. Here one can see for the A κ -scan the already mentioned number of 12 unphysical minima below the EW one, while for the v R -scan up to 14 unphysical minima below the EW one are visible. The lower right plot as a function of v R shows the interesting feature that for a certain value of v R ≈ 803 GeV two of the unphysical minima bifurcate into two different ones each. Hence, in the v R -scan either 12 or 14 unphysical minima below the EW minimum are present (see also the discussion below). Even though there are numerous unphysical minima present below the EW vacuum for all the investigated parameter points, most of the unphysical minima do not give rise to a rapid decay of the EW vacuum, because the tunneling rates for all possible transitions turn out to be very small. This is in particular the case also for the global minimum of each parameter point, which is characterized by non-zero values of the right-handed sneutrino fields φ iR and large negative values of the Higgs fields φ d and φ u . As a result of the latter, the tunneling rates for transitions into the global minimum are suppressed by the large distance in field space between the EW minimum with positive values of φ d and φ u .
The displayed results clearly indicate that a detailed investigation of the tunneling rates is crucial for determining the constraints from vacuum stability on the parameter space of the µνSSM. Analyses in which just the minima of the potential are determined and parameter points where the EW vacuum is not the global one are discarded, as they have sometimes been carried out in the literature for other models, would obviously be completely misleading for the case of the µνSSM.
As another interesting feature that is visible in Fig. 2, one can observe in the scan over A κ the presence of a local minimum with a vanishing potential value V = 0 over the whole scan range. This minimum is located at the origin of field space, and it is therefore of particular interest. 17 According to our approach for investigating vacuum stability, we find that the presence of this minimum at the origin of field space does not give rise to exclusions of parameter points in this scan, because the tunneling rate from the EW vacuum into this minimum is too small. It should however be noted that the analysis carried out here relies on the assumption that the EW minimum was actually adopted by the universe in its thermal history. If a local minimum at the origin exists in the universe now, one may wonder, however, whether such a local minimum at the origin may have existed already during the thermal history of the early universe. While in our present analysis of vacuum stability we restrict ourselves to the case of temperature T = 0, a detailed analysis of the thermal history of the potential would be required in order to determine whether the EW vacuum was actually reached during the thermal evolution. Consequently, an analysis of the evolution of the field values as a function of the temperature (which is inversely proportional to the time) could potentially give rise to even stronger exclusion bounds compared to the zero-temperature analysis performed here (see also the discussion in Sect. 6). The minimum in the origin is also present for parts of the parameter space of the scan over v R . It exists in the range 770 GeV ≤ v R ≤ 802 GeV, as can be seen in the upper right plot of Fig. 2. For larger values of v R the potential values bend down from the horizontal line at V = 0 towards negative potential values, indicating that the minimum shifts away from the origin. The origin becomes unstable along the direction of φ u . Accordingly, in the range 803 GeV ≤ v R ≤ 840 GeV we find an unphysical minimum with non-zero values of 53 GeV ≤ |φ u |≤ 377 GeV and all other fields vanishing.
The above discussion leads to the following important conclusions. We have found that the µνSSM gives rise to a rich vacuum structure, featuring a metastable (i.e., long-lived in comparison to the age of the universe) or short-lived EW vacuum together with several deeper unphysical minima. Thus, the analysis of the validity of a certain parameter point should include a test of vacuum stability. Furthermore, as mentioned above, it would not be sufficient to simply require that the EW minimum is the global minimum of the potential, as such a criterion would exclude large parts of actually allowed parameter space. Instead, a calculation of the tunneling rates of the EW vacuum has to be performed in order to reliably determine whether a parameter point is phenomenologically viable. In the following we further discuss the structure of the potentially dangerous unphysical minima, including an analytical discussion of the dependence on the main parameters entering the calculation of S E , in order to determine which parts of the µνSSM parameter space are expected to be susceptible to vacuum instabilities.
We note in this context that parameter points where the EW vacuum is short-lived feature transitions into unphysical minima which are not the deepest ones. On the contrary, we can see from the plots in the upper row of Fig. 2 that the most dangerous minima are in fact the ones with the highest values of V , which therefore have the smallest potential difference compared to the EW minimum, ∆V EW . This indicates that for the considered case ∆V EW is not the most important quantity in the evaluation of the transition rate. Instead, as can be seen in the lower row of Fig. 2, the distance in field space ∆φ = |φ EW − φ| between the EW minimum at φ EW and an unphysical minimum at φ has a bigger impact. In the lower row of Fig. 2, ∆φ is shown as a function of A κ on the left and of v R on the right for all potentially dangerous (i.e., deeper than the EW vacuum) unphysical minima that occur for each parameter point. It can be seen that the transitions with the lowest values of S E are typically those into the unphysical minima with the lowest values of ∆φ (as expected, a smaller distance in field space is correlated with a higher decay rate of the EW vacuum).
Our results indicate that for the analyzed µνSSM scenarios the distance in field space is more important for the determination of the tunneling rate of the EW vacuum than the difference of the potential depths.
In the lower right plot of Fig. 2, and to a lesser extent also in the lower left plot, one can see that a few points are missing in the curves (in particular near the bifurcation point in the lower right plot). This is due to the fact that the numerical approach applied here to determine the stationary conditions missed one of the minima (see also the discussion in Sect. 4.1). Nevertheless, given the fact that this happened only in rare cases, our conclusions regarding the impact of vacuum stability constraints on the µνSSM are not affected by this issue. In practice, it was always possible to find also the minima missing in Fig. 2 by marginally changing the input parameters. However, we decided to show the results as they were obtained during a single run in order to give an impression of the performance of the applied procedure.
The fact that the field space distance ∆φ appears to be more important than the potential difference ∆V EW , and also that the dangerous minima were found at distances ∆φ ≈ v R , makes it apparent that the field values of the right-handed sneutrino fields are particularly important for the determination of the tunneling rates. In order to shed more light on the field configurations of the unphysical minima, we show in Tab. 2 the field values of the most dangerous minima for three selected parameter points together with the field values of the EW minimum (φ d and φ u in the EW minimum are determined by tan β = 2 and v ≈ 246 GeV). We only display the unphysical minima for which the transition rate from  the EW minimum corresponds to values of S E below 700. The first parameter point under consideration (A) has A κ = −880 GeV and v R = 800 GeV, and it features the unphysical minima AI and AII with S E < 700. The most striking feature of these minima is that, while the EW minimum has universal field values φ iR = v R , in AI and AII the field value of one of the right-handed sneutrino fields is almost zero, while the other two approximately maintain values of about v R . A similar observation can be made for the second parameter point (B), with A κ = −900 GeV and v R = 800 GeV, featuring the unphysical minima BI-BIV. Here, the minima BIII and BIV are the analogues to the minima AII and AI of the first parameter point, respectively, i.e. they lie on the same branch of points in the plots in Fig. 2. In addition, two more minima are shown for the parameter point B, where only one of the fields φ iR has a value of ≈ v R , while the other two have φ iR ≈ 0. For example, BI lies exactly on the φ 3R axis. Two minima of this kind are also present for the parameter point A, but they correspond to EW vacuum decay rates with S E > 700 and are therefore not shown in the table. The second point has a larger value of |A κ |. This suggests that more negative values of A κ can give rise to an unstable EW vacuum (as also shown in Fig. 2). Finally, the third point (C) with A κ = −880 GeV and v R = 800 GeV features the unphysical minima CI-CIV. While they show a similar field configuration as BI-BIV, it is interesting to compare the different parameters related to the tunneling probabilities. For this reason we show in Tab. 3 the values of S E associated to each minimum, together with the coefficients that define the form of V EW as defined in Eq. (26). Values of S E giving rise to a (potentially) short-lived EW vacuum are highlighted in bold font. In addition, the table also lists the values for the field space distance ∆φ and the potential difference −∆V EW between the EW minimum and the unphysical minima. One can see that all three points have at least one unphysical deeper minimum associated with a value of S E < 440. For AII we find S E = 401.5, such that for the parameter point A the EW vacuum is potentially short-lived.   In order to definitively answer the question whether point A features a sufficiently long-lived EW vacuum or not, one would have to improve the computation of S E , incorporating in particular a more elaborate treatment of the prefactor K introduced in Eq. (28) (see also the discussion below). For the parameter point B, the presence of both BIII and BIV with S E = 241.1 and S E = 380.4, respectively, indicates that the EW vacuum is short-lived, and the parameter point should be rejected. The same observation can be made for the parameter point C, for which we find S E = 353.89 for CII and S E = 180.7 for CIV. One can gain further analytical insight into the tunneling dynamics by comparing the minima of different parameter points with similar field configurations, i.e. minima that lie on the same branch of points in Fig. 2. As already mentioned before, for the point A with the minimum AI and S E = 660.0 the analogues minimum for the point B is the minimum BIV with S E = 380.4. Hence, the change from A κ = −880 GeV to A κ = −900 GeV gives rise to an instability of the EW vacuum. The reason for this can be understood analytically by realizing that the only field that substantially changes if the EW vacuum decays into AI and BIV is φ 3R . In the EW minium we have φ 3R = v R = 800 GeV, and the field evolves to φ 3R ≈ 0 during the transition. Thus, the potential V EW as defined in Eq. (26) is given for this kind of transitions in very good approximation by choosing the unit field vectorφ i in the direction of −φ 3R , such thatφ AI,BIV = (0, 0, 0, 0, −1, 0, 0, 0). Then the coefficients of ∆V EW are given by where in the second row the terms ∼ v iL and ∼ Y ν ij were neglected, and where in our scans v R = v 1R = v 2R = v 3R . Furthermore, given the small value of λ 3 = 0.05 in this scenario, the first term in the expression for m 2 is suppressed compared to the second term. For the determination of S E the ratio A 2 /(m 2 λ) is particularly important (see also the definition of δ in Eq. (27)), which is related to the fact that terms in V EW with odd powers of the fields give rise to the potential barrier that separates the EW minimum from the unphysical minimum: increasing the trilinear coefficient A leads to larger values of S E , while S E becomes smaller with increasing coefficients m 2 and λ of the terms with even powers of the fields. Focusing on the terms ∼ A κ 3 , we find that this ratio scales with −κ 3 |A κ 3 |/v R . 18 Thus, if |A κ 3 | becomes larger, the decrease in the potential barrier caused by A is more important than the decrease of the bilinear coefficient m 2 , such that S E , and therefore the lifetime of the EW vacuum, becomes smaller.
It should be noted that the change from an unstable to a metastable EW vacuum sensitively depends on the value of A κ , as we demonstrated here for the points A and B, where the A κ values differ by only ≈ 2% for points with a metastable and an unstable vacuum. Consequently, one can expect that the uncertainty that is related to our treatment of the prefactor K in Eq. (28), giving rise to the intermediate regime 390 < S E < 440 in which the lifetime of the EW vacuum is of comparable size as the age of the universe, affects only a relatively small fraction of the µνSSM parameter space. For most parts of the parameter space, the value of S E for possible vacuum decays should either be substantially below or above the intermediate regime with S E ≈ 400, such that its impact on vacuum stability can clearly be determined. 19 The smallest values of S E are found for the parameter point C, which has a smaller value of v R = 780 GeV compared to v R = 800 GeV for the points A and B. We already mentioned before that this can intuitively be understood due to the fact that the field space difference ∆φ becomes smaller for the relevant unphysical minima. We find, for example, a value of S E = 353.89 for the minimum CII, while for the analogous minimum AII of the first parameter point with an identical value of A κ we find S E = 401.5. Following the above reasoning, the smaller S E value for CII compared to AII can be understood by the fact that (for λ 3 1) A 2 /(m 2 λ) scales with 1/v R . Comparing to Tab. 3, one can additionally see that also the absolute value of the potential difference |∆V EW | changes with v R . However, since the rate for the vacuum decay into CII is larger compared to the decay into AII, even though |∆V EW | is larger for AII, one can conclude (as before) that the field space difference ∆φ has quantitatively a larger impact on the lifetime of the EW vacuum.
The most dangerous minima of all three parameter points cause EW vacuum decays for which φ 1R (instead of φ 3R ) evolves from φ 1R = v R to φ 1R ≈ 0. For these transitions an important difference arises from the fact that λ 1 = 0.56 λ 3 , such that the interactions between φ 1 and the Higgs doublet fields φ d and φ u cannot be neglected. As a result, also the latter fields change during the transition from the EW vacuum into the vacua AII, BIII and CIV, while φ d and φ u remain approximately constant for the transitions into AI, BIV and CII. For the cases where the doublet fields change during the transition, the expressions for 18 Note that A κ i < 0 is required in order to avoid tachyonic CP-odd states. 19 A similar observation was made in Ref. [37] regarding charge-and color-breaking minima in the MSSM and their relation to the precise values of the soft trilinear couplings. the coefficients of V EW become significantly more complicated and the analytic treatment does not provide much insight. 20 Instead we show in Fig. 3 the potential difference ∆V EW in the direction of the unphysical minimum corresponding to the smallest value of S E for different values of A κ and v R .
In the left plot of Fig. 3 we show the potential difference ∆V EW for three different values of A κ = −900, −880, −860 GeV and for fixed v R = 800 GeV, while the remaining parameters are fixed as shown in Tab. 1. Hence, the curve with A κ = −880 GeV (orange) corresponds to the parameter point A discussed in relation to Tab. 2 and Tab. 3, and ∆V EW is shown in the direction of the minimum AII. One can see that the variation of A κ induces a variation of both the barrier height and the depth of the unphysical minimum, whereas the field space difference ∆φ (indicated by the vertical dashed lines) is largely unaffected. This follows the expectation from Eq. (30): For a transition which at least partially evolves in the direction −φ iR , i = 1, 2, 3, there is a positive contribution to the coefficient A which is proportional to v iR κ 2 i , and a negative contribution proportional to −κ i |A κ i |. Thus, the smaller |A κ |, the larger is the potential barrier between both minima for fixed values of v iR and κ i , and the lifetime of the EW vacuum increases. On the other hand, for large values of |A κ i | the lifetime of the EW vacuum decreases, and constraints arising from the vacuum stability become important.
This observation is not surprising since it is known that for |A κ i | v iR κ i the CP-even right-handed sneutrinos can become tachyonic [23,32], pointing to the fact that vacuum instabilities might occur. Similar observations were also made in the NMSSM where only one gauge singlet field is present [40,78]. 21 In the µνSSM, the presence of three such fields leads to the fact that there are more ways in which dangerous unphysical minima can arise. For instance, there can be minima in which one, two or all three fields φ iR take on values of ≈ 0. In addition, the way in which the singlet states are coupled to the doublet Higgs fields yields a much richer structure of vacuum configurations, with several options that could give rise to EW vacuum instabilities. Thus, compared to the NMSSM, the constraints from vacuum instabilities can be expected to be even more important in the µνSSM.
In the right plot of Fig. 3 we show ∆V EW in the direction of the most dangerous unphysical minimum for different values of v R = 780, 800, 820 GeV and a fixed value of A κ = −880 GeV. As before, the orange curve belongs to the parameter point A of Tab. 2, with ∆V EW shown in the direction of the minimum AII. There are two main effects of the variation of v R on the stability of the EW vacuum. As already mentioned before, smaller values of v R reduce the value of ∆φ (indicated by the vertical dashed lines in Fig. 3) for minima of the type shown in Tab. 2, in which the field values of one or more right-handed sneutrinos go from v R to zero during the transition. Hence, smaller values of v R are also associated with a smaller lifetime of the (metastable) EW vacuum. This effect is further enhanced by the fact that also the potential barrier becomes smaller when the value of v R becomes smaller.
In combination of both observations discussed above, i.e., smaller values of v iR and larger values of |A κ i | can lead to an unstable EW vacuum, there is a clear phenomenological consequence. One can expect that an instability of the EW vacuum in the µνSSM can be caused by the presence of relatively light right-handed sneutrinos, with masses at or below the EW scale. Since this is the parameter region in which there are prospects to potentially probe the existence of these particle states at the LHC or other future colliders, either directly or indirectly via signal rate measurements of the Higgs boson at 125 GeV, the constraints from vacuum stability can have an important impact on scenarios that can be probed by collider experiments, and thus should be taken into account in analyses of the collider prospects of the µνSSM.
As an illustration of our results, we show in Fig. 4 the parameter points in the v R -A κ plane, in which the parameter space is divided into a red region in which the EW vacuum is short-lived, an orange region in which the EW vacuum is potentially short-lived (see the discussion in Sect. 4.2 about this intermediate region), and a blue region in which the EW vacuum is metastable and long-lived. There are no points in the lower left corner due to the appearance of tachyonic CP-even states at tree level (see also Fig. 1). For none of the points the EW minimum is the global minimum of the potential. Demanding that the EW minimum should be the global minimum would imply the exclusion of the phenomenologically viable blue region. This illustrates the importance of taking into account the possibility of a sufficiently long-lived metastable EW vacuum. Since the lifetime of the EW vacuum has a sensitive dependence on the precise values of v R and A κ , the orange region constitutes only a small narrow band of the investigated parameter space. As the parameter points in the red region feature a short-lived EW vacuum, they are not physically viable and should be excluded. 21 The heuristic NMSSM criterion A 2 κ > 9m 2 S [78], with m 2 S being the soft mass parameter of the singlet scalar, relies on the condition that the EW minimum is the global minimum of the potential, and does not take into account the possibility of a metastable EW vacuum. The lightest (loop-corrected) scalar masses that we find in the red, blue and orange regions are m h 1 = 70.77 GeV, 112.4 GeV and 116.5 GeV, respectively. As a result, we find that in this scenario all points with m h 1 < 112 GeV are excluded due to EW vacuum stability constraints. This demonstrates the importance of the interplay between the constraints from vacuum stability and possible collider phenomenology of the model at current or future colliders. As is indicated by the gray area, only two small parameter regions are excluded via cross section limits from direct searches, and both of these regions lie within the red parameter region. Hence, we observe in the scenario investigated here that the vacuum stability constraints have a larger impact and give rise to exclusion limits that by far exceed the present limits from direct searches for the BSM Higgs bosons of the µνSSM.
In Fig. 4 we also indicate the parameter points for which the prediction for the mass of the SM-like Higgs-boson does not lie within the range (125 ± 3) GeV (pink area). However, the points in this region should not be considered as strictly excluded, since m h 125 could easily be adjusted to lie within the required interval by modifying the parameters related to the stop sector, from which the radiative corrections to m h 125 mainly arise. Related to m h 125 , one can also see that large parts of the green area, which indicates a worse fit to the signal-rate measurements of h 125 , overlap with the pink area. The green area is defined as the parameter region in which the points have ∆χ 2 = χ 2 µνSSM − χ 2 SM > 5.99, i.e. the points are disfavored at the 95% confidence level in the two-dimensional parameter scan considered here. In contrast to the exclusions from direct searches, the exclusion limits related to the properties of h 125 (pink and green areas) stretch out over parts of the blue parameter region.

Conclusions and outlook
In this paper we have presented an investigation of the stability of the EW vacuum of the µνSSM. We have described in detail the approach that was used in order to determine the dangerous unphysical minima in the eight-dimensional field space of the CP-even scalar fields. Moreover, we have utilized a numerically robust semi-classical approximation for the computation of the EW vacuum decay rates, which allowed us to categorize parameter points for which the EW minimum is not the global minimum into points featuring a metastable long-lived or an unstable EW vacuum (as well as an intermediate region between the two cases).
Focusing on the alignment without decoupling limit of the model, we have demonstrated in our numerical discussion how the analysis of vacuum stability can provide important constraints on the parameter space. We have generically found that the EW minimum is not the global minimum of the potential. In comparison to the NMSSM, the presence of three gauge singlet scalar fields in the µνSSM leads to the possibility of various potentially dangerous unphysical minima existing in the neutral scalar potential for a single parameter point. We have found that the largest rates of the EW vacuum decays are related to the unphysical minima that are the closest in field space to the EW minimum, despite the fact that there might be even deeper unphysical minima present. Furthermore, our results clearly indicate that simply requiring that the EW minimum is the global minimum of the potential would exclude large parts of the µνSSM parameter space that are actually phenomenologically viable. Instead, it is crucial to take into account the possibility of a metastable EW vacuum and to perform an analysis of the EW vacuum decay rates.
We have shown analytically and numerically that an unstable EW vacuum can be expected to occur quite generically in regions of the parameter space in which relatively light right-handed sneutrinos are present. Due to the fact that these particles are gauge singlets, thus not coupled directly to the SM particles, they may not be detectable at collider experiments even if they have small masses. It is therefore an important finding that the corresponding parameter space can be constrained in other ways, such as the vacuum stability investigation as presented here. On the other hand, for the case where the sneutrinos have a significant coupling to the SM particles via a mixing with the Higgs boson at 125 GeV, constraints from collider experiments in combination with constraints from the analysis of the EW vacuum stability can be utilized in a complementary way in order to exclude parts of the parameter space of the µνSSM.
As an outlook to possible future research, we emphasize that further constraints on the parameter space of the µνSSM may be obtained by complementing the present analysis at T = 0 with an investigation of the thermal history of the universe, see recent studies in other scalar extensions of the SM [40,41]. While the EW minimum might be sufficiently long-lived at zero temperature, the universe might have adopted one of the unphysical minima during its thermal history. It can then happen that the transition rate for the phase transition into the EW vacuum would have never been large enough to allow for the onset of EW symmetry breaking, and the universe would have been trapped until zero temperature in one of the unphysical minima. This vacuum trapping effect is especially important if the formation of vacua at the origin, or vacua in which only the singlet fields obtain vacuum expectation values, are possible, and we encountered such cases in our analysis. It should be noted that minima of the latter kind usually form at temperatures much larger than the EW scale, because vanishing field values of the singlet fields are not stabilized via the SM interactions. In the event that an unphysical minimum at the origin (or any minimum that features non-zero vevs only for the right-handed sneutrino fields) can exist until zero temperature there are then two possibilities: The universe is trapped in such a vacuum, or the universe undergoes a first-order EW phase transition. While the vacuum trapping might yield important constraints on the parameter space of the µνSSM, the possibility of first-order EW phase transitions can facilitate baryogenesis or give rise to the formation of a stochastic gravitational wave background in the early universe, thus pointing towards particularly interesting regions of parameter space. An investigation of the thermal history of the µνSSM is left for future work. The definition of each item of the dictionaries can be found in Ref. [34]. The most timeconsuming call is the one of the function check stability, which depending on the input parameters usually takes roughly a minute to finish. The results of the vacuum stability check are saved in:

B Coefficients of V EW for special cases
We consider here two special cases, for which we provide explicit expressions for the coefficients of the scalar potential relative to the EW vacuum, V EW (ϕ,φ), as given in Eq. (26). If one of the right-handed sneutrino fields φ iR is mixed with the SM-like Higgs boson, then transitions along the directions −φ iR and φ u occur (see AII, BIII and CIV in Tab. 2). Choosing for instance i = 1, the unit vectorφ for such transitions can be approximately expressed asφ = N (0, D, 1/D, 0, 0, 0, 0, 0) , with N = D √ D 4 + 1 , where D parametrizes the direction in field space, and the variation of the other fields was neglected. Then the coefficients of V EW are given by (34)  . Another kind of transition that was encountered in the numerical discussion was one in which two of the right-handed sneutrino fields, e.g., φ 1R and φ 2R , evolve from v R to approximately zero (see BI and CI in Tab. 2). Neglecting the variation of the Higgs doublet fields φ d and φ u for these transitions, one can parametrize the unit vector bŷ ϕ = −N (0, 0, D, 1/D, 0, 0, 0, 0) , and one finds for the coefficients As a cross check one can use the feature that since both φ 1R and φ 2R are gauge-singlet fields, the expressions shown above are invariant under the transformation D → 1/D and at the same time exchanging the indices 1 ↔ 2.