Impact of Vacuum Stability Constraints on the Phenomenology of Supersymmetric Models

We present a fast and efficient method for studying vacuum stability constraints in multi-scalar theories beyond the Standard Model. This method is designed for a reliable use in large scale parameter scans. The minimization of the scalar potential is done with the well-known polynomial homotopy continuation, and the decay rate of a false vacuum in a multi-scalar theory is estimated by an exact solution of the bounce action in the one-field case. We compare to more precise calculations of the tunnelling path at the tree- and one-loop level and find good agreement for the resulting constraints on the parameter space. Numerical stability, runtime and reliability are significantly improved compared to approaches existing in the literature. This procedure is applied to several phenomenologically interesting benchmark scenarios defined in the Minimal Supersymmetric Standard Model. We utilize our efficient approach to study the impact of simultaneously varying multiple fields and illustrate the importance of correctly identifying the most dangerous minimum among the minima that are deeper than the electroweak vacuum.


Introduction
In the Standard Model (SM) the vacuum state is characterised by a non-zero vacuum expectation value (vev) of the Higgs field arising from the postulated form of the Higgs potential. This vacuum breaks the electroweak (EW) symmetry and gives masses to the particles via the Brout-Englert-Higgs mechanism. In the standard model of big bang cosmology, the phase transition into this EW vacuum happens within the first instants of the existence of the universe. As the EW vacuum is formed at that time, its stability is therefore required on time scales of the lifetime of the universe.
In the SM the EW vacuum is stable at the EW scale by construction. However, when extrapolating the model to high scales this behaviour can change as a consequence of the running of the quartic Higgs coupling λ [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Based on the present level of the theoretical predictions and the experimental inputs on the top-quark and the Higgs-boson mass, an instability occurs at scales 10 10 GeV with a lifetime that is significantly larger than the age of the universe. Such a configuration where the lifetime of the vacuum state is larger than the age of the universe is called "metastable" and is considered to be theoretically valid. The determination of the lifetime of the metastable vacuum follows the procedure developed in [20,21]. The presence of additional scalar degrees of freedom, as predicted in many models of physics beyond the SM (BSM), is in principle expected to improve the stability of the potential at high scales through its impact on the beta function of the quartic Higgs coupling. On the other hand, additional scalar degrees of freedom have the potential to destabilise the vacuum already at the EW scale by introducing additional minima of the scalar potential.
In BSM theories the assessment of vacuum stability at the EW scale places important constraints on the parameter space of the model under consideration. A sufficient condition to ensure vacuum stability at the EW scale is to require the EW vacuum to be the global minimum (i. e. true vacuum) of the scalar potential. In this case, the EW vacuum as the state with lowest potential energy is favoured and thus absolutely stable. In case the EW vacuum is a local minimum (false vacuum), the corresponding parameter region may still be considered as allowed if it is metastable. On the contrary, any configuration that predicts a shorter lifetime than the age of the universe is considered to be unstable and thus excluded as it is inconsistent with the observed lifetime of the universe.
In rather simple models, such as the Two-Higgs-Doublet Model, analytic conditions for absolute stability have been derived [22][23][24][25][26][27][28]. In theories with more scalars analytic approaches can still be applied to a simplifying subset of fields. However, conclusions about vacuum stability can change severely when additional degrees of freedom are considered (see e.g. [29]). Thus, numerical approaches that can account for a variety of fields simultaneously are of interest. Supersymmetry (SUSY) requires an extended Higgs sector as compared to the SM case and adds a scalar degree of freedom for every fermionic one in the SM. Therefore even the Minimal Supersymmetric Model (MSSM) corresponds to a multi-scalar theory with a much richer scalar sector than the SM. Many different approaches have been employed in order to obtain constraints from vacuum stability for the MSSM .
Public tools that allow one to efficiently obtain constraints from vacuum stability in general BSM models would clearly be very useful in this context. To our knowledge Vevacious [50,57], which is designed to check the stability of the EW vacuum including one-loop and finite temperature effects, is the only public tool that is applicable to a variety of BSM models. In this paper we present an approach that provides a highly efficient and reliable evaluation of the constraints from vacuum stability such that they can be incorporated into BSM parameter scans, which typically run over a large number of points in a multi-dimensional parameter space. Our approach is applicable to any model with a renormalisable scalar potential and has been validated on SUSY and non-SUSY models. In this work we outline our method and, as an example, apply it to the MSSM benchmark scenarios for Higgs searches at the LHC that were recently defined in [58]. These benchmark scenarios were designed to provide interesting Higgs phenomenology at the LHC. Limits from searches for SUSY particles as well as other constraints affecting the MSSM Higgs sector were taken into account in the definitions of the benchmark scenarios, while no detailed investigation of the parameter planes of the scenarios with respect to constraints from vac-uum stability has been carried out. We perform such an analysis of the various benchmark planes and indicate the parameter regions that are incompatible with the requirement of a sufficiently stable EW vacuum. A public tool based on the implementation of our method with which the numerical results in this paper have been obtained is in preparation.
This paper is organised as follows: we begin with the description of our method in Section 2 by discussing the most general renormalisable scalar potential in a form that is suitable for vacuum stability calculations and compare different methods for estimating the vacuum lifetime. In Section 3, we fix our notation for the MSSM and the relevant field space. Finally, we apply our method to the MSSM in Section 4 and discuss in detail the results for three of the benchmark scenarios defined in [58]. We conclude in Section 5. Furthermore, Appendix A contains the full MSSM scalar potential and Appendix B illustrates our results for the remaining CP-conserving benchmark scenarios proposed in [58].

Vacuum Stability in Multidimensional Field Spaces
The vacuum state of a (quantum) field theory is determined by the state of lowest potential energy. In field theory, this state is the minimum of the (effective) potential V (φ), which describes the potential energy density of a field φ. In general, the Lagrangian for such a real scalar field φ is given by where V can be an arbitrary function of the field φ that is bounded from below. In a renormalisable quantum field theory at tree-level, it may contain all interactions up to quartic terms. Formally, the effective potential is defined for classical field values φ cl that minimise the effective action. For our purpose, the field theoretical potential and the effective potential are the same when replacing field operators φ by classical commuting field values φ cl and defining the effective potential V (φ = φ cl ) as function of φ ≡ φ cl [59]. Thus, we treat all scalar fields as commuting variables.
We consider now the general case of n real scalar fields φ a with a ∈ {1, . . . , n} in a renormalisable quantum field theory at tree-level where the sum over repeated indices is implied. The totally symmetric coefficient tensors λ abcd , A abc , m 2 ab and t a as well as the constant c contain all possible real coefficients with non-negative mass dimension. This potential includes in general up to 3 n stationary points 1 out of which an initial vacuum at φ = φ v is selected as minimum with ∂V ∂φ a φ= φv = 0 , (2.3) and the mass matrix where t a vanishes due to Eq. (2.3), and we have normalised the potential energy at ϕ = 0 to zero. For particle physics applications, this normalisation plays no role. Note, however, that a constant term yields a non-vanishing cosmological constant [60]. We rewrite the fieldspace vector as ϕ → ϕφ with a unit vectorφ and its absolute value ϕ = ϕ 2 1 + . . . + ϕ 2 n and obtain where all the dependence on the normalised direction in field spaceφ has been absorbed into the coefficients λ, A and m 2 . The potential has to be bounded from below, so λ > 0 for all directionsφ. Furthermore, in order to have a minimum at ϕ = 0, the condition m 2 > 0 has to be satisfied. There is a freedom of sign choice in either ϕ or A. We always choose A > 0 without loss of generality and have defined the minus sign in Eq. (2.6) such that a possible minimum will be located in ϕ > 0. Figure 1 shows the resulting possible shapes of the potential in Eq. (2.6). Since it is a quartic polynomial in one variable it can have at most two minima, one of which we have chosen to lie at the origin. The second minimum exists as soon as and is deeper than the minimum at the origin if This discussion implies that large cubic terms A compared to the mass parameters and self-couplings are potentially dangerous for the stability of the initial vacuum at the origin. We call the directionsφ fulfilling Eq. (2.8) deep directions. This simple form is very useful for the calculation of vacuum decay in Section 2.1. However, many disjoint regions of deep directions may exist which makes the numerical search for such directions on the unit (n − 1)-sphere of directionsφ infeasible, see e. g. Ref. [56]. We instead use the numerical method of polynomial homotopy continuation (PHC) (see e. g. [61] or [62]) to find all stationary points of Eq. (2.2). From these stationary points we select the deep directions by comparing their depth to the initial vacuum.
PHC efficiently finds all solutions of systems of polynomial equations. We use it to solve and find all real solutions, i. e. the stationary points of the scalar potential. While PHC in theory never fails to find all solutions of the system, solutions may be missed due to numerical uncertainties in judging whether a solution is real or complex. This can be avoided by a careful preconditioning of the system of equations [61]. Another subtlety is that PHC only finds point-like, isolated solutions. This is especially important in the physically interesting cases of gauge theories where any vacuum is only unique up to gauge transformations. If any gauge freedom is left in the model, this turns all isolated solutions into continuous curves which cannot be found by the algorithm. For this reason it is essential to implement models with all gauge redundancies removed. For the case of at least one Higgs doublet this can be achieved by setting the charged and imaginary components of one Higgs doublet to zero without loss of generality.

Calculation of the Bounce Action
We briefly review the definition of the so-called bounce action, which describes the decay of a false vacuum. Consider a single real field Lagrangian as in Eq. (2.1). The semi-classical tunnelling and first quantum corrections were calculated in [20,21]. It was found that the decay rate Γ of a metastable vacuum state per (spatial) volume V S is given by the exponential decay law where K is a dimensionful parameter that will be specified below, and B denotes the bounce action which gives the dominant contribution to Γ. The bounce φ B (ρ) is the solution of the euclidean equation of motion with the boundary conditions U is the euclidean scalar potential, ρ is a spacetime variable and φ v is the location of the metastable minimum. The bounce action B is the stationary point of the euclidean action given by the integral In the one field case, Eqs. (2.11) and (2.13) can be solved numerically by the undershoot/overshoot method (see e.g. [63]). While all of the above equations generalise trivially to the multi-field case φ → φ, the strategy for obtaining the decay rate becomes considerably more involved. In order to judge the stability of the EW vacuum we need to obtain the minimal bounce action for tunnelling into a deeper point in the scalar potential. There exist methods for solving Eq. (2.11) numerically in multiple field dimensions [38,[64][65][66][67][68][69][70][71][72] using optimization, discretisation, path-deformation or multiple shooting. As we will show in the next section a fast evaluation of the bounce action is more important for our purposes than an extremely precise result. For this reason, we approximate the path of the bounce by the straight line in a given deep direction. The potential along this straight path is a simple quartic polynomial as given by Eq. (2.6). For this form of the potential there exists a semi-analytic result for the bounce action [ This dependence on the field normalisation arises from the equation of motion, where Eq. (2.11) only applies to fields with canonically normalised kinetic terms. A consistent expansion of the form of Eq. (2.6) therefore requires all real field components to have canonically normalised kinetic terms. It is crucial to ensure that the implementation of the scalar potential fulfils this requirement.

Lifetime of the Metastable Vacuum
The vacuum lifetime in Eq. (2.10) also depends on the quantity K. The value of K is both challenging to calculate and a subdominant effect towards the tunnelling rate as it does not enter in the exponent. Since it is a dimensionful parameter, [K] = GeV 4 , it can be estimated from a typical scale M of the theory as Comparing the vacuum decay time τ decay with the age of the universe t uni [74] yields [17] Figure 2 shows the relative lifetime τ decay /t uni as a function of B and M. As expected, the threshold of instability where τ decay ∼ t uni is highly sensitive to B and only mildly sensitive to M. In Fig. 2 we also show the contours corresponding to a 5 σ expected decay or a 5 σ expected survival of the vacuum during the evolution of the universe. The survival probability is given by where the (space-time) volume of the past light-cone isṼ light-cone ∼ 0.15/H 4 0 [10], and H 0 is the current value of the Hubble parameter [74]. The points in the green region of Fig. 2 are definitely long-lived with respect to the age of the universe, while the red points are definitely short-lived. We see that varying the scale M over a generous range from 10 GeV to 100 TeV shifts the border between metastability and instability by less than 10% in B. We therefore consider any point where B > 440 as long-lived and any point where B < 390 as short-lived. We treat the intermediate range 390 < B < 440 as an uncertainty on the stability threshold from the unknown M.

Parameter Scans and the Reliability of Vacuum Stability Calculations
Our main objective is to enable the use of vacuum stability constraints in fits and parameter scans of BSM theories. Typical BSM parameter scans consider millions of different parameter points which places strong requirements on evaluation time and reliability.
The first trade-off between speed and precision is in the calculation of the bounce action described in Section 2.1. Instead of using one of the available sophisticated solvers [71,72]which may take a lot of runtime and could encounter numerical problems-we approximate the tunnelling path with a straight line in field space and use the semi-analytic solution Eq. (2.14) [73]. A comparison between such simple approximations and the multipleshooting method of [72] has been performed in [75] where agreement within O(10%) for polynomial potentials has been found. The approximation according to Eq. (2.14) is evaluated instantaneously while the available solvers typically take between a few seconds and several minutes per tunnelling calculation.
This approximation, as well as the discussion of Section 2, relies on the potential being a quartic polynomial which is only true at tree-level. In contrast, vacuum stability in the SM has been studied up to full NNLO [1,2,7,10,11,14,17] precision involving the two-loop effective Higgs potential and NNLO running and threshold effects. Even for BSM models the public tool Vevacious [50,57] can perform vacuum stability calculations using the oneloop Coleman-Weinberg potential. However, it has recently been shown [76,77] that the use of the loop-corrected effective potential for stability calculations is not in general a consistent perturbative expansion. This happens because the effective action is a perturbative expansion in both the usual powers of and the momentum transfer, where the effective potential corresponds to the zeroth order term of the momentum expansion. Truncating this second expansion does not in general provide a good approximation in calculations of the bounce action. Therefore, higher momentum terms of the full effective action can give contributions to the bounce action comparable to the contributions from the effective potential. Since it seems unfeasible to calculate even the leading higher momentum terms of the effective action in general BSM models it appears questionable whether using the one-loop effective potential for stability calculations leads to more precise results. For this reason we stick to the tree-level potential which allows us to apply the very concise formulation of Sections 2 and 2.1 and considerably increases the speed and numerical stability of our calculation.
The tunnelling time with respect to the age of the universe as given in Eq. (2.19) depends exponentially on the value of the bounce action B. For any given parameter point, small uncertainties on B are therefore amplified to large uncertainties on the tunnelling time. While this makes precise predictions for the lifetime of individual parameter points very challenging, it is less problematic for constraining the parameter space of BSM models since the bounce action B is also very sensitive to the values of the model parameters. Therefore, a small shift in parameter space typically leads to a change in the bounce action substantially larger than the uncertainties described above. For this reason the resulting constraints on the model parameter space depend only mildly on the precise way B is calculated. This dependence can be estimated from the width of the uncertainty band 390 < B < 440 (see Section 2.2) that we will show in all of our results (see e.g. Figs. 4 and 9 below).

Application to the Minimal Supersymmetric Standard Model
Constraints from vacuum stability can have an important impact on supersymmetric models. The main reason is their abundance of scalar fields with at least two Higgs doublets as well as two scalar superpartners for each SM fermion. The MSSM is the simplest and best studied supersymmetric model where only a second Higgs doublet is added to the Higgs sector [78][79][80]. The most important vacuum stability constraints in the MSSM concern instabilities in directions with sfermion vevs. These are commonly referred to as charge or colour breaking (CCB) vacua . The existing results include both analytic and semianalytic studies of specific directions in field space as well as fully numeric approaches. After specifying our notation for the scalar potential of the MSSM we are going to illustrate our treatment of vacuum stability constraints for an example application to MSSM benchmark scenarios for Higgs searches.

Treatment of the MSSM Scalar Potential
For our discussion, we focus on the third generation of SM fermions and their corresponding superpartners as those have the largest couplings to the Higgs sector. We comment on the impact of the sfermions of the first and second generations below. The superpotential of the MSSM including only third generation fermion superfields is given by , the left-chiral superfields containing the SM quark and lepton doublets Q L = (t L , b L ) and L L = (ν L , τ L ), respectively, as well as the superfields containing the SU (2) L singletst R ,b R andτ R . The Yukawa couplings are denoted by y t,b,τ and the dot product is the SU j with the totally antisymmetric tensor ab , where 12 = −1. The superpotential gives rise to the F -terms contributing to the scalar potential, where the sum runs over all scalar components of the superfields in Eq. (3.1). Note that the F -terms contain quadratic, cubic and quartic terms.
Additional supersymmetric contributions to the scalar potential come from the gauge structure of the model. These D-terms are given by where the sum over φ runs over all scalar components as in Eq.
where y t,b,τ are the Yukawa couplings of Eq. (3.1). We shall express the soft breaking parameter B µ via the mass m A of the CP-odd Higgs boson, using the ratio of the vevs of the two Higgs doublets at the EW vacuum The full scalar potential of the MSSM including all Higgs and third generation sfermion fields is thus given by see Appendix A for the explicit expression. We have so far written V in terms of complex fields, while the method outlined in Section 2.1 relies on a reduction of the field space to a direction parametrised by a single real scalar field. For Eq. (2.14) to be applicable, this field also needs to have a canonically normalised kinetic term. We ensure both requirements by expressing V exclusively through real scalar fields with canonically normalised kinetic terms and by expanding all complex scalar fields as This ensures that ϕ is canonically normalised after expanding ϕ → ϕφ to obtain Eq. (2.6) as long as |φ| = 1. In this notation the EW vacuum is given by GeV is the SM Higgs vev. It would be unfeasible to vary all real scalar degrees of freedom simultaneously since the runtime of the minimisation procedure scales exponentially with the number of fields considered. 2 For our studies of the MSSM we combine all stationary points found by varying the three sets of fields All sets contain the real parts of the neutral Higgs fields that participate in EW symmetry breaking. The first set additionally contains the realt andb fields, the second set thet and τ fields and the third set theb andτ fields. This method will not be able to find stationary points for whicht,b andτ vevs are simultaneously non-zero. The distance in field space between the EW vacuum and another minimum is expected to increase as more fields take non-zero values at this second minimum. We therefore neglect these configurations since the tunnelling time increases with the field-space distance. Moreover, we foundν vevs and vevs of the first and second sfermion generations to have no impact on the observed constraints. Therefore we are not going to show and discuss them in detail in the following, but we will comment on their impact below. We also do not take charged or CP-odd Higgs fields and the imaginary parts of the sfermion fields into account. Ignoring the CP-odd and charged Higgs directions is motivated by the absence of any spontaneous CP or charge breaking in the Higgs sector of the 2HDM [22,23] (and thus the MSSM). While non-zero charged and CP-odd Higgs vevs can in principle develop in the presence of sfermion vevs we found no region of parameter space where these are relevant. We neglect the imaginary parts of the sfermions as they are not expected to add new features in the absence of CP-violation. 3 Note that this discussion is specific to the MSSM. In the NMSSM for example the different kinds of Higgs vevs are expected to be more relevant [81].

The Impact of Yukawa Couplings
The third generation Yukawa couplings are the largest Yukawa couplings and sensitively depend on tan β already at the tree level. Their value is determined via the quark masses 4 2 For example, considering the set of fields in Eq. (3.10a) yields ∼ 10 times longer runtimes than varying the two sets of fields Apart from the CP-violation in the CKM matrix which does not enter our study. 4 We treat the quark masses as running masses at the SUSY scale. To this end we use RunDec [82][83][84] to run the MS quark masses to the SUSY scale assuming SM running. and the vev of the Higgs doublet coupling to them at the tree level: For small tan β, the value of the bottom Yukawa coupling is suppressed with respect to the top Yukawa coupling, while for large tan β they become comparable. For large tan β, the bottom Yukawa coupling is very sensitive to SUSY loop corrections that are enhanced by tan β. The leading corrections can be resummed [85][86][87][88], and it is advantageous to include them despite the fact that the scalar potential is evaluated at the tree level, since they effectively change the value of the bottom Yukawa coupling. The impact of the resummed corrections on the Yukawa coupling can be included by replacing y tree b in Eq. (3.11) by where ∆ b contains SUSY loop corrections. The dominant contributions arise from the gluino-sbottom and higgsino-stop loop, which enter in the sum where mt 1,2 , mb 1,2 are the masses of thet andb mass eigenstates, M 3 denotes the gluino mass, and µ is the higgsino mass parameter. The function C(x, y, z) is given by (3.14) These ∆ b corrections lead to an enhancement of y b especially for large µ < 0 and A t , M 3 > 0 as both contributions are negative in this case and reduce the denominator of Eq. (3.12). For ∆ b → −1 the bottom Yukawa coupling gets pushed into the nonperturbative regime. Taking Eq. (3.12) into account can lead to important effects ofb vevs in the MSSM, see [55], as will be visible in our numerical analysis below. We also take into account a similar but numerically smaller effect for the Yukawa coupling of the τ lepton, y τ [86].

Constraints from Vacuum Stability in MSSM Benchmark Scenarios
In the following we are going to present an example application of our method for obtaining vacuum stability constraints. We will study vacuum stability constraints for some of the MSSM benchmark scenarios defined in [58] to illustrate the impact of the constraints and compare our method to previous approaches.

Vacuum Stability in the M 125
h Scenario The first benchmark scenario defined in [58] is the M 125 h scenario. It features rather heavy SUSY particles with a SM-like Higgs boson at 125 GeV and can be used to display the sensitivity of searches for additional Higgs bosons at the LHC. Its parameters are given by while m A and tan β are varied in order to span the considered parameter plane. The soft SUSY breaking parameters A t,b,τ vary as a function of tan β for fixed X t . Note that the gaugino mass parameters M 1,2,3 only enter our analysis through the ∆ b and ∆ τ corrections. Figure 3 shows the vacuum stability analysis in this benchmark plane.
In the left panel of Fig. 3, the colour code indicates the lifetime of the EW vacuum at each point in the parameter plane. In the dark green region the EW vacuum is the global minimum of the theory, and the EW vacuum is stable in this parameter region. The light green area depicts regions where deeper minima exist but the lifetime of the false EW vacuum is large compared to the age of the universe (see Section 2.2). For these parameter points the EW vacuum is metastable and the parameter points are allowed. For points in the red region, on the other hand, the tunnelling process is fast, and they are excluded as the EW vacuum is short-lived. The small yellow region contains all points in the intermediate region discussed in Section 2.2 where an estimate of the uncertainties gives no decisive conclusion on the longevity of the false vacuum. This plot of Fig. 3 left shows that the M 125 h benchmark plane is hardly constrained by the requirement of vacuum stability. Only a parameter region with small values of tan β 1 can be excluded.
The middle panel of Fig. 3 shows the character of the most dangerous minimum (MDM), i. e. it displays which fields acquire non-vanishing vevs at this vacuum. The MDM is defined as the minimum with the lowest bounce action for tunnelling from the EW vacuum. One can see that for small values of tan β or both moderate tan β and m A the MDM is a CCB minimum witht vevs (yellow). For larger values of tan β and m A a minimum withb vevs takes over (blue). This behaviour is expected as for higher tan β the couplings of the Higgs sector to d-type (s)quarks are enhanced which also increases the impact ofb vevs on our vacuum stability analysis.
The right panel of Fig. 3 displays the character of the global minimum for the M 125 h benchmark plane. It is important to note that the MDM and the global minimum of the scalar potential can in general differ from each other. 5 We see that the global minimum has onlyb vevs for most of the parameter space, while there is a large region where the MDM involvest vevs.

Impact of the Trilinear Terms
We noted in Section 2 that the parameters entering the cubic terms of the potential are expected to be especially important for the stability of the EW vacuum. Since m A is related to a bilinear term in the potential, and tan β mostly affects the quartic Yukawa couplings, we switch to a different slice of the parameter space which is more relevant for the stability studies. We start from a point in the m A -tan β plane of the M 125 h scenario that is absolutely stable and given by tan β = 20 , m A = 1500 GeV .  [58]. Among these, the non-observation of heavy Higgs bosons decaying into τ pairs [89,90] is the most relevant constraint. In contrast to the m A -tan β plane of the M 125 h scenario, we now vary the parameters µ and A ≡ A t = A b = A τ starting from this point. Figure 4 shows the vacuum stability analysis in this new parameter plane.
The left panel of Fig. 4 indicates the lifetime of the EW vacuum. The colour coding is the same as in the corresponding plot of Fig. 3, i. e. red points depict short-lived configurations, while the EW vacuum for light green points is metastable, and the EW vacuum in the dark green area is stable. The thin yellow band indicates the uncertainty band of 390 < B < 440 discussed in Section 2.2. The EW vacuum becomes more and more unstable for larger absolute values of µ and A. For small values of these parameters the potential is absolutely stable with a region of long-lived metastability in between. Note that also in this parameter plane the yellow uncertainty region of 390 < B < 440 corresponds to only a thin band between long-and short-lived regions. The point marked by the × is the starting point in the M 125 h plane depicted in Fig. 3. One can see that in the plane of Fig. 4 this point is close to a region of metastability, but quite far from any dangerously short lived  The character of the MDM, i. e. the fields that acquire non-zero vevs in this vacuum, is shown in the middle panel of Fig. 4. It is dominated by yellowt vevs in this plane, but bluẽ b vevs are also important for large negative values of µ. The ∆ b corrections described in Section 3.2 are enhanced in this parameter region and have a large impact. They are also the cause of the tachyonic region for large negative µ and positive A. Between thet-vev and b-vev regime a region appears (shown in green) wheret andb vevs occur simultaneously. The small blue region for µ > 0 is visible because the more dangerous minima witht vevs only appear for slightly higher values of A and µ, and the globalb-vev minimum is the only other vacuum in this parameter region besides the EW vacuum.
In the right panel of Fig. 4, the fields which acquire non-zero vevs at the global minimum are indicated. In this parameter plane, there are large regions with simultaneoust andτ vevs at the global minimum. Through most of the plane the fields acquiring vevs differ between the MDM and the global minimum. The green region of simultaneoust andb vevs which is visible in the middle panel of Fig. 4 does not correspond to the global minimum of the theory. This is expected as additional large quartic F and D-term contributions appear if multiple kinds of squarks take on non-zero vevs simultaneously. These are positive contributions to the scalar potential that lift up these regions of field space. No such contributions appear in the case of simultaneous squark and slepton vevs which is why the orange regions of simultaneoust andτ vevs are present in the right panel of Fig. 4. Note that the quartic F and D-term contributions do not prevent the minima with mixedt andb vevs from being the MDM as Fig. 4 (centre) shows. However, for the parameter plane considered here these minima featuring simultaneoust andτ vevs have no impact on the stability constraints of Fig. 4 (left).
We finally comment on the impact of the detailed field content, in particular the first and second generation sfermions and theν vevs, on these results. The small Yukawa couplings of these particles tend to push any additional minima to very large field values, which renders these configurations long-lived. As a consequence, the metastability bound (yellow region in Fig. 4, left) and the character of the MDM (Fig. 4, centre) are insensitive to the impact of those fields. On the other hand, the character of the global minimum may be significantly affected by fields with relatively small couplings to the Higgs sector. Indeed, if we were to includec ands vevs they would dominate the global minimum (Fig. 4, right) through most of the parameter plane. They would even cut slightly into the edges of the stable dark green region turning it long-lived. Our analysis therefore shows that neither the investigation of just the region of absolute stability nor of the character of the global minimum yields reliable bounds from vacuum stability. This is due to the fact that both these quantities can sensitively depend on the considered field content, where even very weakly coupled scalar degrees of freedom can have a significant impact. Instead, the correct determination of the boundary between the short-lived and the long-lived region crucially relies on the correct identification of the MDM, which in general can be very different from the global minimum. This boundary, and accordingly the constraint on the parameter space from vacuum stability, is in fact governed by the fields with the largest Yukawa couplings and therefore insensitive to effects from particles with a small coupling to the Higgs sector.

Comparison to Semi-Analytic Bounds and Existing Codes
We now compare our results shown in Fig. 4 with results from the literature. An approximate bound for MSSM CCB instabilities including vacuum tunnelling is given by [36,38] 7.5 long-lived.
is sometimes used to judge whether a parameter point might be sufficiently long-lived (see e. g. the discussion in [91]). The public code Vevacious [50,57] can calculate the lifetime of the EW vacuum in BSM models using the tree-level or Coleman-Weinberg one-loop potential, optionally including finite temperature effects, see e. g. [49]. semi-analytic bounds stable long-lived heuristic  Fig. 4 in [38]). Moreover, the dependence on tan β [52] is not included in the approximate bound. Another reason is that onlyt-related parameters enter Eq. (4.3), while our analysis shows that alsob vevs have important effects in this case.

bound in Eq. (4.3) becomes less reliable for values of m
The heuristic bound Eq. (4.4) (dotted black contour in Fig. 5) should also be compared to the yellow region. While there are clear differences in shape, the size of the long-lived region in our result roughly matches the heuristic bound. This can be qualitatively understood from Eq. (2.8). For λ ∼ O(1) this yields A/m > 2 as a bound for absolute stability. Therefore A/m > 3 as a bound for metastability appears to be a reasonable estimate. While we only show these comparisons for one parameter plane they hold very similarly for every plane we have studied. Our comparison shows that all of these approximate bounds have deficiencies in determining the allowed parameter region, and dedicated analyses are necessary to obtain more reliable conclusions.
Next we compare our results (Fig. 6, left) to the tree-level (Fig. 6, centre) and one loop (right) results of Vevacious. In the Vevacious runs we have taken into account only the fields from Eq. (3.10a) since we found no relevant constraints fromτ vevs in this plane. Vevacious by default considers tunnelling to the minimum which is closest in field space to the EW vacuum. In the newest version (1.2.03+ [92]) one can optionally consider tunnelling to the global minimum instead. In generating Fig. 6 we combined the results from both of these approaches by choosing the option giving the stronger bound at each individual point. One obvious difference between our results and Vevacious are the metastable regions that Vevacious finds for µ ∼ 3 TeV and |A| ∼ 5 TeV. In this region Vevacious considers the wrong minimum to be the MDM. The global minimum (withb vevs, see Fig. 4, right) is closest to the EW vacuum in field space. Therefore, Vevacious can only consider tunnelling into this minimum instead of a slightly further and shallower minimum witht vevs which Vevacious one-loop short-lived long-lived stable gives the stronger constraints shown in our results. A similar issue is responsible for the edge in the Vevacious result around µ ∼ −2.5 TeV and A ∼ 4 TeV. A second kind of visible difference is the absence in the Vevacious result of the bumps in the long-lived region in our result around |µ| ∼ 2 TeV and |A| ∼ 5 TeV. The optimization of the bounce action by CosmoTransitions [71], which is used by Vevacious, leads to a slightly stronger and more reliable metastability bound in this region. 6 Apart from these deviations our results are in good agreement with the tree-level results of Vevacious. The deviations for individual points and the rugged edges of the light green region in the Vevacious result are likely signs of numerical instability. This especially includes the isolated red points in the light green region which result from numerical errors in the calculation of the tunnelling time.
The comparison with the Vevacious results using the Coleman-Weinberg one-loop effective potential at zero temperature (Fig. 6, right) shows that the one-loop effects on the allowed parameter space are small for this scenario. The one-loop result from Vevacious clearly suffers from numerical instabilities. However, the stable region is nearly identical to the tree-level results, and the long-lived region is similarly sized as the tree-level Vevacious result with differences in shape. The long-lived region appearing around A ∼ 5 TeV and µ ∼ 4 TeV as well as the missing region around A ∼ −5 TeV and the spikes around µ ∼ −2.5 TeV are consequences of the same MDM misidentification as in the Vevacious treelevel result (see previous paragraph). Comparing the runtime of our code to the runtime of Vevacious in this parameter plane including only the field set of Eq. (3.10a) we find our tree-level code to be ∼ 5 times faster than the tree-level and ∼ 200 times faster than the one-loop Vevacious run.

Parameter Dependence of the Vacuum Structure and Degenerate Vacua
The dashed line in Fig. 4 is the line where X t has the same value as in the benchmark plane, Fig. 3. The mass m h of the SM-like Higgs boson depends dominantly on the parameters tan β, X t and the stop masses. We therefore expect the Higgs mass to stay close to 125 GeV when moving away from the point × along this line. 7 We use this as motivation to further investigate the vacuum structure along this line. Figure 7 shows the depth of the stationary points of the scalar potential as a function of µ along this line. The constant depth of the EW vacuum is shown in grey while the other colours indicate the CCB stationary points. Note that not only local minima, but all stationary points including saddle points and local maxima are shown in Fig. 7. The dashed line indicates the MDM for each value of µ.
It can be seen from Fig. 7 that for large negative µ simultaneoust andτ vevs (orange) dominate the global minimum for the considered field content until theτ vevs at these stationary points approach zero around µ = −2.2 TeV, and puret vevs take over. From µ ≈ −1.8 TeV onwards the EW vacuum is the global minimum until a CCB vacuum with b vevs appears at µ ≈ 1.6 TeV. The MDM, on the other hand, is the second deepestbvev minimum for µ −3.5 TeV, before switching to thet-vev minimum, followed by the window of absolute stability µ ∈ [−1.8 TeV, 1.5 TeV]. For positive values of µ > 1.5 TeV the instability first develops towards the globalb-vev minimum until thet-vev minimum takes over at µ ≈ 2 TeV.
In Fig. 7 several stationary points with multiple kinds of sfermion vevs appear. Stationary points with mixed squark and slepton vevs can be deeper than the corresponding stationary points with only one type of vev (this can be seen for instance by comparing the deepest stationary point with oranget andτ vevs to the ones with yellowt vevs and red τ vevs). A stationary point with multiple kinds of squark vevs, however, is always higher than one with only one kind of the involved vevs. This is due to the additional positive quartic contributions to the potential for stationary points with both kinds of squark vevs.
Another feature visible in Fig. 7 is theb-vev stationary point approaching the EW vacuum at µ 0 from above. In this regime the ∆ b corrections significantly enhance the bottom Yukawa coupling giving rise to a large mixing in theb sector and a corresponding decrease of one of theb masses. The depth of this stationary point becomes degenerate with the EW vacuum. For even larger negative µ oneb squark becomes tachyonic, and the EW vacuum turns into a saddle point. The plot ends before this happens (corresponding to the white region in Fig. 4) as we require the existence of an EW vacuum.
The scalar potential Eq. (3.7) for our field sets Eqs. (3.10a) to (3.10c) has two accidental Z 2 symmetries. The potential is symmetric under simultaneous sign flips of the left-and right-handed sfermions of a kind and under simultaneous sign flips of all doublet components This results in sets of degenerate and physically equivalent stationary points related by these symmetries. 8 Since the EW vacuum is also invariant under Eq. (4.5) the tunnelling time to minima related by this symmetry is always identical. However, since the EW vacuum breaks Eq. (4.6) 9 the tunnelling time into two stationary points related through this transformation can differ. In most cases, whichever of these two points is closer in field space to the EW vacuum gives the lower value for B. Note that this is not a small effect. The values of B for stationary points related by Eq. (4.6) can differ by more than an order of magnitude. This effect has recently been studied for the simpler case of a 2HDM in [100].

Vacuum Stability in the M 125
h (τ ) Scenario A benchmark scenario with lightτ has been proposed in [58] under the name M 125 h (τ ). It is defined by (4.7) The scenario differs from the M 125 h scenario of Eq. (4.1) only in greatly reduced softτ masses with a correspondingly reduced Aτ . However, µ is not reduced and is now µ ∼ 3mτ .  According to Eq. (4.4) we would therefore expect vacuum stability constraints to be relevant in the M 125 h (τ ) benchmark plane. The authors of [58] used Vevacious to check for vacuum instabilities in this scenario and found a short-lived region in the parameter space for large tan β and small m A .
Our results shown in Fig. 8 confirm these observations. Figure 8 (left) shows a shortlived region for large tan β. This region extends towards smaller values of tan β in the low m A regime as noted in [58] but also in the region of large m A . Also visible is the small region of instability for tan β < 1 noted in Fig. 3. The MDM for the instability at large tan β is a vacuum withτ vevs as can bee seen from Fig. 8 (centre). Compared to Fig. 3 the absolutely stable region is additionally reduced by at-τ -vev minimum appearing around m A ∼ 1 TeV and tan β < 30. The minima withb vevs, which were the MDM for large regions of the M 125 h scenario, are entirely replaced by minima withτ vevs. Only a very small purple region with simultaneousb andτ vevs at the MDM exists. In Fig. 8 (right) the global minimum withb vevs is very similar to Fig. 3 (right). Only at larger m Awhere the EW vacuum in the M 125 h scenario was absolutely stable -global minima withτ vevs are now present. Our results in the M 125 h (τ ) scenario compared to the M 125 h scenario illustrate that constraints from vacuum stability indeed become relevant when the cubic terms in the scalar potential become larger than the quadratic terms. Since µ ∼ 3mτ in this scenario, a further increase of µ or a decrease of mτ could render the M 125 h (τ ) scenario entirely short-lived. In the region of µ ∼ 3mτ chosen in the M 125 h (τ ) scenario, the vacuum stability constraints show a significant dependence on the parameters m A and tan β.

Vacuum Stability in the M 125 h (alignment) Scenario
In this section we turn to another scenario from [58], the M 125 h (alignment) scenario, where constraints from vacuum stability turn out to have a very large impact. The scenario is   with m A ∈ [100 GeV, 1 TeV] and tan β ∈ [1,20]. It is chosen to accommodate light additional Higgs bosons through the so-called alignment without decoupling behaviour [91,[101][102][103][104][105][106][107]. Alignment without decoupling in the phenomenologically interesting region of relatively small tan β requires µ, A t mt. This requirement-according to Eq. (4.4)-has already been noted to be problematic for vacuum stability in [58]. Indeed we find that the EW vacuum in this scenario is short-lived through all of its parameter space with instabilities in directions witht vevs. The corresponding Fig. 11 can be found in Appendix B.
In order to assess how far away this scenario is from a region of metastability we again select a phenomenologically interesting point in its parameter plane, and vary A = A t = A b = A τ and µ starting from this point. The resulting plane is shown in Fig. 9. The colour code and the quantities shown in the sub-plots are the same as in Figs. 3 and 4. In the left panel we have superimposed the contours for the mass of the light Higgs boson calculated with FeynHiggs 2.14.2 [93][94][95][96][97][98][99]. This shows that in order to obtain a long-lived scenario while keeping the correct Higgs mass one could for example change µ = 7.5 TeV → 4 TeV and A = 6.25 TeV → 5 TeV . (4.10) Other choices are possible within the theoretical uncertainty of the Higgs mass prediction. The alignment-without-decoupling behaviour implies that the 125 GeV Higgs boson has SM-like properties. In the MSSM this behaviour arises from a cancellation between treelevel mixing contributions and higher-order corrections in the Higgs-boson mass matrix. Shifting the parameters towards the region where the vacuum is metastable like it is done in Eq. (4.10) could affect these cancellations and therefore spoil the alignment properties. Indeed, with FeynHiggs we obtain ∼ 10% enhanced couplings of h to bottom and tau pairs leading to a ∼ 10% reduced BR(h → γγ) compared to the SM for the parameter point of Eq. (4.10). While the change in the couplings indicates that the cancellation that is present in the alignment without decoupling regime does not fully occur for the shifted point, this parameter point is nevertheless still compatible with the LHC Higgs measurements within the present uncertainties [108].
In the middle panel of Fig. 9 one can see that the MDM directions through most of the parameter space are directions witht vevs. If the instabilities were in directions withb orτ vevs, the requirement of a long-lived vacuum could have been achieved by increasing m D 3 ,L 3 ,E 3 and/or decreasing A b,τ while keeping µ and thet sector unchanged such that the phenomenology would not be much affected. However, our analysis shows that the parameter region associated with the alignment-without-decoupling behaviour gives rise to minima witht vevs that render the electroweak vacuum short-lived. Accordingly, the behaviour of alignment without decoupling and the requirement of a long-lived vacuum are in some tension with each other, since adjusting m Q 3 ,U 3 and A t to ensure vacuum stability in the stop directions would change the mass and phenomenology of the light Higgs boson h. The right plot of Fig. 9 illustrates that the global minimum has non-vanishingτ vevs through significant parts of the parameter plane, while for large values of A, µ 5 TeVb vevs take over. However, the deepest minima appear to be nearly degenerate as can be seen from Fig. 10. The depth of the stationary points of the scalar potential along the line of constant X t = 5 TeV (which is indicated in Fig. 9 centre and right) is shown in Fig. 10. It includes all stationary points in the selected field sets that are at least as deep as the EW vacuum. The EW vacuum is shown in grey, and the CCB stationary points are distinguished by the other colours. The plot illustrates that there is no stable region along this slice of parameter space, in accordance with Fig. 9. As already pointed out, the MDM along this line is a minimum witht vevs through most of the parameter range, with the exception of the region with small values of µ where the MDM hasτ vevs. There exist stationary points that are deeper than the MDM withb vevs 10 ,t-τ vevs andτ vevs with theτ -vev minimum being the global minimum until theb-vev minimum takes over for µ 5.5 TeV. These deeper minima are however very far from the EW vacuum in field space with high barriers and have no impact on the tunnelling rate. For large values of µ stationary points with botht andb vevs develop. In agreement with the previous discussions these stationary points are less deep than the stationary points with only one kind of squark vev. All points along this line of constant X t = 5 TeV would be long-lived if the MDM witht vevs was absent.
A detailed analysis of the question whether adjustments of the proposed scenarios for alignment without decoupling could bring these scenarios into better agreement with the constraints from vacuum stability while retaining their alignment properties goes beyond the scope of this work. We leave this issue for future studies.

Summary and Conclusions
In this work we have presented a fast and efficient method for determining the constraints from vacuum stability on the model parameters of multi-scalar theories beyond the Standard Model. The stability of the EW vacuum at the tree level is investigated using polynomial homotopy continuation for the minimization of the scalar potential and approximating the decay rate of the EW vacuum into a deeper vacuum by an exact solution of the bounce action in the one-field case. The method has been designed to combine a fast evaluation with high numerical stability, enabling a reliable use in large scale parameter scans. The generic approach admits a wide range of applications in many different BSM models. We have argued that this approach is appropriate for applications in multi-dimensional parameter spaces since the dependence of the vacuum stability constraint on the model parameters is typically much stronger than the impact of high precision calculations of the bounce path and the tunnelling action. Therefore, the constraints on the parameter space are relatively insensitive on the details of the calculation of the lifetime, and uncertainties in the classification of long-and short-lived configurations only affect small regions in the parameter space.
In order to evaluate the constraints from vacuum stability it is in principle desirable to simultaneously take into account all possible vevs of all the scalar fields in the considered model. However, this often becomes impractical in models with large scalar sectors. We have illustrated our approach by determining the constraints from vacuum stability for a set of benchmark scenarios for Higgs searches in the MSSM that have recently been proposed in [58]. Our efficient approach has allowed us to consider up to six real degrees of freedom simultaneously. We have included the possible vevs of the particles with the largest Yukawa couplings (namelyt,b andτ ) in addition to the Higgs vevs when searching for unstable directions. We have observed that this yields reliable vacuum stability constraints-especially when including multiple kinds of sfermions simultaneously. However, for the studied scenarios we would have obtained very similar vacuum stability constraints if we had considered each kind of sfermion separately as an approximation.
The result that the sfermions of the first and second generation and the sneutrinos have a minor impact on the determination of the boundary between long-lived and short-lived vacua is related to the fact that their smaller Yukawa couplings lead to additional minima at very large field values, which renders these configurations long-lived. On the other hand, the global minimum-and to some extent also the region of absolute stability-shows a significant dependence on additional field content that couples only weakly to the Higgs sector. We also found several parameter regions where the global minimum is characterised by vevs of different sfermions that are simultaneously non-zero.
As a result, our analysis shows that neither the investigation of just the region of absolute stability nor of the character of the global minimum is sufficient to obtain reliable bounds from vacuum stability. Instead, the determination of the boundary between the short-lived and the long-lived region crucially relies on the correct identification of the most dangerous minimum (MDM), which is the minimum with the shortest tunnelling time from the EW vacuum. The MDM often differs from the global minimum.
For the considered M 125 h MSSM benchmark scenario we have found that the impact of vacuum stability constraints is small in the m A -tan β parameter plane defining the scenario. However, a variation of the parameters µ and A around their values chosen in the benchmark scenario shows an important impact of vacuum stability constraints. In this plane we have illustrated that the most dangerous and the global minimum are in general different. We have furthermore stressed the importance of corrections to the relation between the bottomquark mass and the bottom Yukawa coupling in certain regions of parameter space. These ∆ b corrections can significantly enhance the value of the bottom Yukawa coupling and thus triggerb-vev instabilities.
We have also used this parameter plane to compare our results with existing studies and codes. Comparing our results to approximate analytic vacuum stability bounds we have seen that those approximations can serve as a rough estimate of the effect of vacuum stability constraints. However, they cannot capture the complexity of a detailed numerical analysis. Furthermore, we have compared the results of our code to the public code Vevacious. The tree-level results of the codes show some notable differences. The largest differences arise from the determination of the MDM. By default, Vevacious uses the closest minimum in field space-and can be optionally forced to use the global minimum. In contrast, by the use of an analytic expression for the bounce action, we are able to easily calculate all tunnelling times to deeper minima and select the one with the shortest time as the MDM. We have seen that in the benchmark plane under consideration there are regions where the MDM is neither the global minimum nor the minimum closest to the EW vacuum in field space. We have also seen that the impact of the more sophisticated tunnelling path calculation in Vevacious (using the path deformation of CosmoTransitions) is visible but in most cases does not substantially change the boundary between parameter regions with longand short-lived EW vacuum. Calculating the vacuum stability bounds in Vevacious using the one-loop effective potential leads to results qualitatively similar to the tree level results. The one-loop calculation showed strong indications of numerical instability in addition to the problems in the identification of the MDM. Even ignoring the latter problem we have found the potential improvement by the one-loop effective potential on the vacuum stability constraints to be relatively small in view of the significantly increased amount of numerical instabilities and the much longer runtime. Furthermore, as discussed above it is still an open conceptual question whether the application of a loop-corrected effective potential can at all be expected to yield a systematic improvement of vacuum stability constraints compared to a tree-level analysis. The tree-level constraints of our approach using the straight tunnelling path approximation yield numerically stable results in a fraction of the runtime of Vevacious and CosmoTransitions.
We have exemplarily studied the depths of the different stationary points along a line of m h ∼ 125 GeV through the parameter space. This has shown a typical number of stationary points appearing in an analysis of vacuum stability and illustrated the importance of the different scalar degrees of freedom. The MDM was found to coincide with the global minimum only within a restricted part of the relevant parameter range. We have also discussed the impact of degenerate stationary points that are related to each other by a discrete symmetry and found that-in agreement with recent results in the literature [100]if such a symmetry is broken by the EW vacuum the tunnelling times into the degenerate vacua can be vastly different.
In a benchmark scenario with lightτ we have shown important constraints from vacuum stability arising in the m A -tan β plane defining the scenario. In [58] a parameter region with high tan β and low m A of this benchmark plane was identified as being excluded by vacuum stability constraints. Our results provide a more detailed study of these constraints and show that the excluded region extends to high m A if ∆ τ corrections are included.
A further scenario that we have studied in detail is a benchmark scenario in the alignment without decoupling regime. We have found that all parameter points in the benchmark plane yield short-lived EW vacua. We have shown what kind of shift in parameter space could lead to a sufficiently long-lived vacuum while approximately preserving the correct value of the Higgs mass. We found that this naïve approach would not obviously spoil the alignment without decouplings behaviour. The question to what extent the alignment without decoupling behaviour can be realised in phenomenologically viable scenarios where vacuum stability constraints are taken into account would require a detailed study that is beyond the scope of the present work. The remaining CP-conserving benchmark scenarios 11 of [58] behave similarly to the ones we have studied. The vacuum stability constraints in these benchmark planes are shown in Appendix B.
We plan to make our implementation of the described procedure publicly available, which is meant to serve as an easily applicable tool for evaluating vacuum stability constraints that is suitable for the inclusion in scans over multi-dimensional parameter spaces.

A The MSSM Scalar Potential
In this appendix we give the scalar potential of the MSSM in expanded form. We include all fields appearing in Eqs. (3.10a) to (3.10c). The potential is given by consisting of a quadratic (V 2 ), cubic (V 3 ) and quartic (V 4 ) part. The quadratic part is given by the cubic part is given by It should be noted that we have omitted the CP-violating scenario that was proposed in [58] for simplicity only. Our approach is fully applicable also to models containing CP-violating phases. and the quartic part is given by Re(h 0 u ) 2 Re(t R ) 2 Re(h 0 u ) 2 Re(τ L ) 2 + Re(b L ) 2 Re(τ R ) 2 + g 2 1 + g 2 2

(A.4)
B Constraints from Vacuum Stability in the Other Benchmark Scenarios of [58] In this appendix we present our results for the vacuum stability constraints in the remaining CP-conserving benchmark scenarios defined in [58]. The vacuum stability analysis for the M 125 h (alignment) benchmark plane defined in Eq. (4.8) is displayed in Fig. 11. As noted in Section 4.3, all parameter points in this plane give rise to short-lived EW vacua. The MDM is characterised byt vevs throughout the As shown in Fig. 12 the entire benchmark plane features an absolutely stable EW vacuum. This is unsurprising, since it has a very small µ and relatively large soft sfermion masses while the A parameter is not much larger than the soft masses. The small gaugino mass parameters do not lead to instabilities as they only enter the scalar potential through the ∆ b and ∆ τ corrections. In this scenario the heavy Higgs boson H has a mass of about 125 GeV with SM-like couplings. Since the scenario has small soft SUSY-breaking squark mass parameters and a very large µ it is unsurprising that the entire benchmark plane shown in Fig. 13 has a short-lived EW vacuum. Throughout the plane, the MDM has non-vanishingt vevs while the global minimum hast-τ vevs. Considering how constrained the region of phenomenologically viable parameter space for this scenario is (see Fig. 10 of [58]) we do not expect that a   heavy Higgs alignment scenario with a long-lived vacuum exists in the MSSM. However, a detailed assessment of this possibility would require a dedicated study.