Impact of Vacuum Stability Constraints on the Phenomenology of Supersymmetric Models

TL;DR

This paper introduces a fast, reliable method for analyzing vacuum stability constraints in multi-scalar theories, especially within supersymmetric models, improving computational efficiency and accuracy in large-scale parameter scans.

Contribution

The authors develop a novel approach combining polynomial homotopy continuation and exact bounce action solutions to efficiently evaluate vacuum decay rates in supersymmetric models.

Findings

Method shows good agreement with more precise calculations.

Significantly improves numerical stability and runtime.

Highlights importance of correctly identifying the most dangerous vacuum.

Abstract

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.

Figures

Equations

$$\mathcal{L}=\frac{1}{2}{\left(\partial \varphi \right)}^2-V(\varphi ),$$

$$V(\vec{\varphi })={\lambda }_{abcd}{\varphi }_a{\varphi }_b{\varphi }_c{\varphi }_d+A_{abc}{\varphi }_a{\varphi }_b{\varphi }_c+m_{ab}^2{\varphi }_a{\varphi }_b+t_a{\varphi }_a+c,$$

$${\text{diffp}V{\varphi }_a|}_{\vec{\varphi }={\vec{\varphi }}_v}=0,$$

$$M_{ab}={\text{diffp}V{\varphi }_a{\varphi }_b|}_{\vec{\varphi }={\vec{\varphi }}_v}$$

$$V(\vec{\phi })={\lambda }_{abcd}^{\prime }{\phi }_a{\phi }_b{\phi }_c{\phi }_d+A_{abc}^{\prime }{\phi }_a{\phi }_b{\phi }_c+m_{ab}^{\prime 2}{\phi }_a{\phi }_b,$$

$$V(\phi ,\hat{\phi })=\lambda (\hat{\phi }){\phi }^4-A(\hat{\phi }){\phi }^3+m^2(\hat{\phi }){\phi }^2,$$

$${\left(A(\hat{\phi })\right)}^2>\frac{32}{9}{m}^2(\hat{\phi })\lambda (\hat{\phi })$$

$${\left(A(\hat{\phi })\right)}^2>4m^2(\hat{\phi })\lambda (\hat{\phi }).$$

$${\vec{\nabla }}_{\varphi }V=0$$

$$\frac{\Gamma }{V_S}=K{e}^{-B},$$

$$\text{diff}[2]\varphi \rho +\frac{3}{\varphi }\text{diff}\varphi \rho =\text{diffp}U\varphi$$

$$\varphi (\infty )={\varphi }_v,{\text{diff}\varphi \rho |}_{\rho =0}=0.$$

$$B=2{\pi }^2{\int }_0^{\infty }{\rho }^3\text{d}\rho \left[\frac{1}{2}{\left(\text{diff}\rho {\varphi }_B(\rho )\right)}^2+U({\varphi }_B(\rho ))\right].$$

$$B=\frac{{\pi }^2}{3\lambda }{(2-\delta )}^{-3}\left(13.832\delta -10.819{\delta }^2+2.0765{\delta }^3\right)$$

$$\delta =\frac{8\lambda m^2}{A^2}.$$

$$\phi \to n\phi$$

$$\Rightarrow B\to n^{-4}B.$$

$$K={\mathcal{M}}^4.$$

$$\frac{{\tau }_{\text{decay}}}{t_{\text{uni}}}={\left(\frac{\Gamma }{V_S}\right)}^{-\frac{1}{4}}\frac{1}{t_{\text{uni}}}=\frac{1}{t_{\text{uni}}\mathcal{M}}{e}^{B/4}.$$

$$P=\exp \left(-\frac{\Gamma }{V_S}{\tilde{V}}_{\text{light-cone}}\right)=\exp \left(-{\mathcal{M}}^4{\tilde{V}}_{\text{light-cone}}{e}^{-B}\right),$$

$$W=\mu H_u\cdot H_d+y_t{Q}_L\cdot H_u{\bar{t}}_R+y_b{H}_d\cdot Q_L{\bar{b}}_R+y_{\tau }{H}_d\cdot L_L{\bar{\tau }}_R,$$

$$F=\sum _{\varphi }{\left|{\partial }_{x}W\right|}^2,\varphi \in \{h_u^0,h_u^{+},h_d^0,h_d^{-},{\tilde{t}}_L,{\tilde{b}}_L,{\tilde{\tau }}_L,{\tilde{\nu }}_L,{\tilde{t}}_R^{\ast },{\tilde{b}}_R^{\ast },{\tilde{\tau }}_R^{\ast }\}$$

$$D$$

$$D_{{U(1)}_Y}$$

$$D_{{SU(2)}_L}$$

$$D_{{SU(3)}_c}$$

$$V_{\text{soft}}$$

$$+m_{Q_3}^2{\tilde{Q}}_L^{†}{\tilde{Q}}_L+m_{L_3}^2{\tilde{L}}_L^{†}{\tilde{L}}_L+m_{U_3}^2|{\tilde{t}}_R{|}^2+m_{D_3}^2|{\tilde{b}}_R{|}^2+m_{E_3}^2|{\tilde{\tau }}_R{|}^2$$

$$+\left(y_t{A}_t{\tilde{t}}_R^{\ast }{\tilde{Q}}_L\cdot h_u+y_b{A}_b{\tilde{b}}_R^{\ast }{h}_d\cdot {\tilde{Q}}_L+y_{\tau }{A}_{\tau }{\tilde{\tau }}_R^{\ast }{h}_d\cdot {\tilde{L}}_L+\text{h.c.}\right),$$

$$B_{\mu }=m_A^2\sin \beta \cos \beta ,$$

$$\tan \beta =v_u/v_d.$$

$$V=F+D+V_{\text{soft}},$$

$$\varphi \to \frac{1}{\sqrt{2}}\text{Re}(\varphi )+\frac{i}{\sqrt{2}}\text{Im}(\varphi ).$$

Full text

††preprint:

DESY 18–148

TTP 18–036

1 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 $\lambda$ Sher:1988mj ; Casas:1994qy ; Espinosa1995 ; Isidori2001 ; Espinosa2008 ; Ellis2009 ; Degrassi:2012ry ; Bezrukov:2012sa ; Chetyrkin:2013wya ; Buttazzo:2013uya ; Andreassen:2014gha ; Lalak2014 ; Zoller:2014cka ; Bednyakov:2015sca ; Branchina2015 ; Iacobellis2016a ; Andreassen:2017rzq ; Chigusa2017 ; Chigusa2018 . 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 $\gtrsim 10^{10}\,\text{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 Coleman1977 ; Callan1977 . 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 Ferreira2004 ; Barroso2006 ; Barroso2006a ; Maniatis2006 ; Barroso2007 ; Barroso:2013awa ; Barroso2013 . 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. Muhlleitner:2016mzt ). 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 Frere:1983ag ; Claudson:1983et ; Drees:1985ie ; Gunion:1987qv ; Komatsu:1988mt ; Langacker:1994bc ; Casas:1995pd ; Strumia:1996pr ; Kusenko:1996jn ; Casas:1996de ; LeMouel:1997tk ; Abel:1998wr ; LeMouel:2001sf ; Ferreira:2001tk ; Ferreira:2000hg ; Ferreira:2004yg ; Park:2010wf ; Altmannshofer:2012ks ; Carena:2012mw ; Camargo-Molina:2013sta ; Chattopadhyay:2014gfa ; Camargo-Molina2014 ; Chowdhury:2013dka ; Blinov:2013fta ; Bobrowski:2014dla ; Altmannshofer:2014qha ; Hollik:2015pra ; Hollik:2016dcm .

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 Camargo-Molina2013 ; Camargo-Molina2014 , which is designed to check the stability of the EW vacuum including one-loop and finite temperature effects, is the only dedicated public tool that is applicable to a variety of BSM models.111The code BSMPT Basler:2018cwe —while not designed for vacuum stability studies—can also check for absolute stability in several BSM models including one-loop and finite temperature effects. 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 Bahl2018a . 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 vacuum 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 Bahl2018a . 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 Bahl2018a .

2 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(\phi)$, which describes the potential energy density of a field $\phi$. In general, the Lagrangian for such a real scalar field $\phi$ is given by

[TABLE]

where $V$ can be an arbitrary function of the field $\phi$ 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 $\phi_{\text{cl}}$ that minimise the effective action. For our purpose, the field theoretical potential and the effective potential are the same when replacing field operators $\phi$ by classical commuting field values $\phi_{\text{cl}}$ and defining the effective potential $V(\phi=\phi_{\text{cl}})$ as function of $\phi\equiv\phi_{\text{cl}}$ Coleman:1973jx . Thus, we treat all scalar fields as commuting variables.

We consider now the general case of $n$ real scalar fields $\phi_{a}$ with $a\in\{1,\ldots,n\}$ in a renormalisable quantum field theory at tree-level

[TABLE]

where the sum over repeated indices is implied. The totally symmetric coefficient tensors $\lambda_{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 points222Note that we discard complex solutions here as we consider real fields. out of which an initial vacuum at $\vec{\phi}=\vec{\phi}_{v}$ is selected as minimum with

[TABLE]

and the mass matrix

[TABLE]

is positive definite. After expanding Eq. 2 around the vacuum as $\vec{\phi}=\vec{\phi}_{v}+\vec{\varphi}$, with $\vec{\varphi}={(\varphi_{1},\ldots,\varphi_{n})}^{\mathrm{T}}$, we obtain

[TABLE]

where $t^{\prime}_{a}$ vanishes due to Eq. 3, and we have normalised the potential energy at $\vec{\varphi}=0$ to zero. For particle physics applications, this normalisation plays no role. Note, however, that a constant term yields a non-vanishing cosmological constant Weinberg:1988cp . We rewrite the field-space vector as $\vec{\varphi}\rightarrow\varphi\hat{\varphi}$ with a unit vector $\hat{\varphi}$ and its absolute value $\varphi=\sqrt{\varphi_{1}^{2}+\ldots+\varphi_{n}^{2}}$ and obtain

[TABLE]

where all the dependence on the normalised direction in field space $\hat{\varphi}$ has been absorbed into the coefficients $\lambda$, $A$ and $m^{2}$. The potential has to be bounded from below, so $\lambda>0$ for all directions $\hat{\varphi}$. Furthermore, in order to have a minimum at $\varphi=0$, the condition $m^{2}>0$ has to be satisfied. There is a freedom of sign choice in either $\varphi$ or $A$. We always choose $A>0$ without loss of generality and have defined the minus sign in Eq. 6 such that a possible minimum will be located in $\varphi>0$.

Figure 1 shows the resulting possible shapes of the potential in Eq. 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

[TABLE]

and is deeper than the minimum at the origin if

[TABLE]

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 $\hat{\varphi}$ fulfilling Eq. 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 $\hat{\varphi}$ infeasible, see e. g. Ref. Hollik:2016dcm . We instead use the numerical method of polynomial homotopy continuation (PHC) (see e. g. morgan2009solving or Maniatis:2012ex ) to find all stationary points of Eq. 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

[TABLE]

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

2.1 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. 1. The semi-classical tunnelling and first quantum corrections were calculated in Coleman1977 ; Callan1977 . It was found that the decay rate $\Gamma$ of a metastable vacuum state per (spatial) volume $V_{S}$ is given by the exponential decay law

[TABLE]

where $K$ is a dimensionful parameter that will be specified below, and $B$ denotes the bounce action which gives the dominant contribution to $\Gamma$. The bounce $\phi_{B}(\rho)$ is the solution of the euclidean equation of motion

[TABLE]

with the boundary conditions

[TABLE]

$U$ is the euclidean scalar potential, $\rho$ is a spacetime variable and $\phi_{v}$ is the location of the metastable minimum. The bounce action $B$ is the stationary point of the euclidean action given by the integral

[TABLE]

In the one field case, Eqs. 13 and 11 can be solved numerically by the undershoot/overshoot method (see e.g. Apreda:2001us ). While all of the above equations generalise trivially to the multi-field case $\phi\rightarrow\vec{\phi}$, 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. 11 numerically in multiple field dimensions Kusenko:1996jn ; Dasgupta:1996qu ; Moreno:1998bq ; Cline:1999wi ; John:1998ip ; Cline:1998rc ; Konstandin:2006nd ; Park:2010rh ; Wainwright:2011kj ; Masoumi:2016wot ; Athron:2019nbd 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. 6. For this form of the potential there exists a semi-analytic result for the bounce action Adams:1993zs

[TABLE]

with

[TABLE]

The expression in brackets was obtained in Adams:1993zs by fitting a cubic polynomial in $\delta$ to the numerical result. The coefficients do not depend on any model parameters, and the polynomial approximation agrees with the numerical result within a $0.004$ absolute tolerance for all values of $\delta$. We use this formula to calculate $B$ for all deep directions from the initial vacuum. The deep direction with the smallest bounce action is the dominant decay path.

The value of $B$ is not invariant under rescaling of the field $\varphi$ of Eq. 6

[TABLE]

This dependence on the field normalisation arises from the equation of motion, where Eq. 11 only applies to fields with canonically normalised kinetic terms. A consistent expansion of the form of Eq. 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.

2.2 Lifetime of the Metastable Vacuum

The vacuum lifetime in Eq. 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]=\mathrm{GeV}^{4}$, it can be estimated from a typical scale $\mathcal{M}$ of the theory as

[TABLE]

Comparing the vacuum decay time $\tau_{\text{decay}}$ with the age of the universe $t_{\text{uni}}$ Aghanim:2018eyx yields Andreassen:2017rzq

[TABLE]

Figure 2 shows the relative lifetime $\tau_{\text{decay}}/t_{\text{uni}}$ as a function of $B$ and $\mathcal{M}$. As expected, the threshold of instability where $\tau_{\text{decay}}\sim t_{\text{uni}}$ is highly sensitive to $B$ and only mildly sensitive to $\mathcal{M}$. In Fig. 2 we also show the contours corresponding to a $5\,\sigma$ expected decay or a $5\,\sigma$ expected survival of the vacuum during the evolution of the universe. The survival probability is given by

[TABLE]

where the (spacetime) volume of the past light-cone is $\tilde{V}_{\text{light-cone}}\sim 0.15/H_{0}^{4}$ Buttazzo:2013uya , and $H_{0}$ is the current value of the Hubble parameter Aghanim:2018eyx . 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 $\mathcal{M}$ over a generous range from $10\,\mathrm{GeV}$ to $100\,\mathrm{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 $\mathcal{M}$.

2.3 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 Wainwright:2011kj ; Masoumi:2016wot —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. 14 Adams:1993zs . A comparison between such simple approximations and the multiple-shooting method of Masoumi:2016wot has been performed in Masoumi:2017trx where agreement within $\mathcal{O}(10\%)$ for polynomial potentials has been found. The approximation according to Eq. 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 Sher:1988mj ; Casas:1994qy ; Degrassi:2012ry ; Buttazzo:2013uya ; Andreassen:2014gha ; Bednyakov:2015sca ; Andreassen:2017rzq precision involving the two-loop effective Higgs potential and NNLO running and threshold effects. Even for BSM models the public tool Vevacious Camargo-Molina2013 ; Camargo-Molina2014 can perform vacuum stability calculations using the one-loop Coleman–Weinberg potential. However, it has recently been shown Andreassen:2016cff ; Andreassen:2016cvx 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 $\hbar$ 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. 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).

3 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 Nilles:1983ge ; Haber:1984rc ; Gunion:1984yn . 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 Frere:1983ag ; Claudson:1983et ; Drees:1985ie ; Gunion:1987qv ; Komatsu:1988mt ; Langacker:1994bc ; Casas:1995pd ; Strumia:1996pr ; Kusenko:1996jn ; Casas:1996de ; LeMouel:1997tk ; Abel:1998wr ; LeMouel:2001sf ; Ferreira:2001tk ; Ferreira:2000hg ; Ferreira:2004yg ; Park:2010wf ; Altmannshofer:2012ks ; Carena:2012mw ; Camargo-Molina:2013sta ; Chattopadhyay:2014gfa ; Camargo-Molina2014 ; Chowdhury:2013dka ; Blinov:2013fta ; Bobrowski:2014dla ; Altmannshofer:2014qha ; Hollik:2015pra ; Hollik:2016dcm . The existing results include both analytic and semi-analytic 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.

3.1 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

[TABLE]

with the chiral Higgs superfield ${{SU(2)}_{L}}$ doublets $H_{u}=(H^{0}_{d},H_{d}^{-})$ and $H_{d}=(H_{u}^{+},H_{u}^{0})$, the left-chiral superfields containing the SM quark and lepton doublets $Q_{L}=(t_{L},b_{L})$ and $L_{L}=(\nu_{L},\tau_{L})$, respectively, as well as the superfields containing the ${{SU(2)}_{L}}$ singlets $\bar{t}_{R}$, $\bar{b}_{R}$ and $\bar{\tau}_{R}$. The Yukawa couplings are denoted by $y_{t,b,\tau}$ and the dot product is the ${{SU(2)}_{L}}$ invariant multiplication $\Phi_{i}\cdot\Phi_{j}=\epsilon_{ab}\Phi_{i}^{a}\Phi_{j}^{b}$ with the totally antisymmetric tensor $\epsilon_{ab}$, where $\epsilon_{12}=-1$. The superpotential gives rise to the $F$-terms

[TABLE]

contributing to the scalar potential, where the sum runs over all scalar components of the superfields in Eq. 21. 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

[TABLE]

where the sum over $\phi$ runs over all scalar components as in Eq. 22, and $\Phi_{i},\Phi_{j}\in\{h_{u},h_{d},\tilde{Q}_{L},\tilde{L}_{L}\}$ run over the scalar ${{SU(2)}_{L}}$ doublets. The different prefactor for the ${{SU(3)}_{c}}$ $D$-term arises from the sum over the ${{SU(3)}_{c}}$ generators. The $D$-terms only contain quartic terms.

Finally, another contribution to the scalar potential of the MSSM are the soft SUSY breaking terms

[TABLE]

where $y_{t,b,\tau}$ are the Yukawa couplings of Eq. 21. We shall express the soft breaking parameter $B_{\mu}$ via the mass $m_{A}$ of the CP-odd Higgs boson,

[TABLE]

using the ratio of the vevs of the two Higgs doublets at the EW vacuum

[TABLE]

The full scalar potential of the MSSM including all Higgs and third generation sfermion fields is thus given by

[TABLE]

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

[TABLE]

This ensures that $\varphi$ is canonically normalised after expanding $\vec{\varphi}=\varphi\hat{\varphi}$ to obtain Eq. 6 as long as $|\hat{\varphi}|=1$. In this notation the EW vacuum is given by

[TABLE]

where $v=\sqrt{v_{u}^{2}+v_{d}^{2}}\approx 246\,\mathrm{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.333For example, considering the set of fields in Eq. 30a yields $\sim 10$ times longer runtimes than varying the two sets of fields $\{\text{Re}(h_{u}^{0}),\,\text{Re}(h_{d}^{0}),\,\text{Re}(\tilde{t}_{L}),\,\text{Re}(\tilde{t}_{R})\}$ and $\{\text{Re}(h_{u}^{0}),\,\text{Re}(h_{d}^{0}),\,\text{Re}(\tilde{b}_{L}),\,\text{Re}(\tilde{b}_{R})\}$ separately. For our studies of the MSSM we combine all stationary points found by varying the three sets of fields

[TABLE]

All sets contain the real parts of the neutral Higgs fields that participate in EW symmetry breaking. The first set additionally contains the real $\tilde{t}$ and $\tilde{b}$ fields, the second set the $\tilde{t}$ and $\tilde{\tau}$ fields and the third set the $\tilde{b}$ and $\tilde{\tau}$ fields. This method will not be able to find stationary points for which $\tilde{t}$, $\tilde{b}$ and $\tilde{\tau}$ 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 $\tilde{\nu}$ 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 Ferreira2004 ; Barroso2006 (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.444Apart from the CP-violation in the CKM matrix which does not enter our study. 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 Krauss:2017nlh .

3.2 The Impact of Yukawa Couplings

The third generation Yukawa couplings are the largest Yukawa couplings and sensitively depend on $\tan\beta$ already at the tree level. Their value is determined via the quark masses555We treat the quark masses as running masses at the SUSY scale. To this end we use RunDec Chetyrkin:2000yt ; Schmidt2012 ; Herren2018 to run the $\overline{\text{MS}}$ quark masses to the SUSY scale assuming SM running. and the vev of the Higgs doublet coupling to them at the tree level:

[TABLE]

For small $\tan\beta$, the value of the bottom Yukawa coupling is suppressed with respect to the top Yukawa coupling, while for large $\tan\beta$ they become comparable.

For large $\tan\beta$, the bottom Yukawa coupling is very sensitive to SUSY loop corrections that are enhanced by $\tan\beta$. The leading corrections can be resummed Hall:1993gn ; Carena:1994bv ; Pierce:1996zz ; Carena:1999py , 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_{b}^{\text{tree}}$ in Eq. 31 by

[TABLE]

where $\Delta_{b}$ contains SUSY loop corrections. The dominant contributions arise from the gluino-sbottom and higgsino-stop loop, which enter in the sum $\Delta_{b}=\Delta_{b}^{\text{gluino}}+\Delta_{b}^{\text{higgsino}}$:

[TABLE]

where $m_{\tilde{t}_{1,2}}$, $m_{\tilde{b}_{1,2}}$ are the masses of the $\tilde{t}$ and $\tilde{b}$ mass eigenstates, $M_{3}$ denotes the gluino mass, and $\mu$ is the higgsino mass parameter. The function $C(x,y,z)$ is given by

[TABLE]

These $\Delta_{b}$ corrections lead to an enhancement of $y_{b}$ especially for large $\mu<0$ and $A_{t},M_{3}>0$ as both contributions are negative in this case and reduce the denominator of Eq. 32. For $\Delta_{b}\to-1$ the bottom Yukawa coupling gets pushed into the non-perturbative regime. Taking Eq. 32 into account can lead to important effects of $\tilde{b}$ vevs in the MSSM, see Hollik:2015pra , 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 $\tau$ lepton, $y_{\tau}$ Guasch:2001wv .

4 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 Bahl2018a to illustrate the impact of the constraints and compare our method to previous approaches.

4.1 Vacuum Stability in the $M_{h}^{125}$ Scenario

The first benchmark scenario defined in Bahl2018a is the $M_{h}^{125}$ scenario. It features rather heavy SUSY particles with a SM-like Higgs boson at $125\,\mathrm{GeV}{}$ and can be used to display the sensitivity of searches for additional Higgs bosons at the LHC. Its parameters are given by

[TABLE]

while $m_{A}$ and $\tan\beta$ are varied in order to span the considered parameter plane. The soft SUSY breaking parameters $A_{t,b,\tau}$ vary as a function of $\tan\beta$ for fixed $X_{t}$. Note that the gaugino mass parameters $M_{1,2,3}$ only enter our analysis through the $\Delta_{b}$ and $\Delta_{\tau}$ 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_{h}^{125}$ benchmark plane is hardly constrained by the requirement of vacuum stability. Only a parameter region with small values of $\tan\beta\lesssim 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\beta$ or both moderate $\tan\beta$ and $m_{A}$ the MDM is a CCB minimum with $\tilde{t}$ vevs (yellow). For larger values of $\tan\beta$ and $m_{A}$ a minimum with $\tilde{b}$ vevs takes over (blue). This behaviour is expected as for higher $\tan\beta$ the couplings of the Higgs sector to $d$-type (s)quarks are enhanced which also increases the impact of $\tilde{b}$ vevs on our vacuum stability analysis.

The right panel of Fig. 3 displays the character of the global minimum for the $M_{h}^{125}$ 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.666Out of the minima that are deeper than the EW vacuum, the MDM is usually the one that is closest to the EW vacuum in field space. However, this is not always the case, see Fig. 6 below for a counterexample. We see that the global minimum has only $\tilde{b}$ vevs for most of the parameter space, while there is a large region where the MDM involves $\tilde{t}$ vevs.

4.1.1 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\beta$ 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\beta$ plane of the $M^{125}_{h}$ scenario that is absolutely stable and given by

[TABLE]

This point is indicated with a $\boldsymbol{\times}$ in Figs. 3 and 4. It features a Higgs mass of $m_{h}\approx 125\,\mathrm{GeV}$ and is allowed by all the constraints considered in Bahl2018a . Among these, the non-observation of heavy Higgs bosons decaying into $\tau$ pairs Aaboud:2017sjh ; Sirunyan:2018zut is the most relevant constraint. In contrast to the $m_{A}$–$\tan\beta$ plane of the $M^{125}_{h}$ scenario, we now vary the parameters $\mu$ and $A\equiv A_{t}=A_{b}=A_{\tau}$ 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 $\mu$ 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 $\boldsymbol{\times}$ is the starting point in the $M_{h}^{125}$ 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 parameter regions. The missing points in the top-left corner of the plot are points with tachyonic tree-level $\tilde{b}$ masses where the EW vacuum is a saddle point.

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 yellow $\tilde{t}$ vevs in this plane, but blue $\tilde{b}$ vevs are also important for large negative values of $\mu$. The $\Delta_{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 $\mu$ and positive $A$. Between the $\tilde{t}$-vev and $\tilde{b}$-vev regime a region appears (shown in green) where $\tilde{t}$ and $\tilde{b}$ vevs occur simultaneously. The small blue region for $\mu>0$ is visible because the more dangerous minima with $\tilde{t}$ vevs only appear for slightly higher values of $A$ and $\mu$, and the global $\tilde{b}$-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 simultaneous $\tilde{t}$ and $\tilde{\tau}$ 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 simultaneous $\tilde{t}$ and $\tilde{b}$ 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 simultaneous $\tilde{t}$ and $\tilde{\tau}$ 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 mixed $\tilde{t}$ and $\tilde{b}$ vevs from being the MDM as Fig. 4 (centre) shows. However, for the parameter plane considered here these minima featuring simultaneous $\tilde{t}$ and $\tilde{\tau}$ 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 $\tilde{\nu}$ 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 include $\tilde{c}$ and $\tilde{s}$ 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.

4.1.2 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 Casas:1995pd ; Kusenko:1996jn

[TABLE]

Furthermore, a “heuristic” bound of

[TABLE]

is sometimes used to judge whether a parameter point might be sufficiently long-lived (see e. g. the discussion in Bechtle:2016kui ).

The public code Vevacious Camargo-Molina2013 ; Camargo-Molina2014 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. Chattopadhyay:2014gfa .

Figure 5 displays our results in comparison to the approximate bounds given in Eqs. 37 and 38, where the contours arising from Eqs. 37 and 38 are superimposed on our results. The solid black contour arising from Eq. 37 should be compared with the edge of the dark green region where the vacuum is stable. This comparison shows significant deviations. This is in particular due to the fact that the absolute stability bound from Eq. 37 considers only the $D$-flat direction. The dashed black contour should be compared with the yellow region at the border between the long-lived (light green) and short-lived (red) regions and shows similar deviations. One source of these deviations is that the metastability bound in Eq. 37 becomes less reliable for values of $m_{\tilde{t}_{R}}^{2}+m_{\tilde{t}_{L}}^{2}\gtrsim{(1200\,\mathrm{GeV})}^{2}$ and large $(A_{\tilde{t}}^{2}+3\mu^{2})$ (see Fig. 4 in Kusenko:1996jn ). Moreover, the dependence on $\tan\beta$ Blinov:2013fta is not included in the approximate bound. Another reason is that only $\tilde{t}$-related parameters enter Eq. 37, while our analysis shows that also $\tilde{b}$ vevs have important effects in this case.

The heuristic bound Eq. 38 (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. 8. For $\lambda\sim\mathcal{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. 30a since we found no relevant constraints from $\tilde{\tau}$ 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+ vevaciousweb ) 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 $\mu\sim 3\,\mathrm{TeV}$ and $|A|\sim 5\,\mathrm{TeV}$. In this region Vevacious considers the wrong minimum to be the MDM. The global minimum (with $\tilde{b}$ 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 with $\tilde{t}$ vevs which gives the stronger constraints shown in our results. A similar issue is responsible for the edge in the Vevacious result around $\mu\sim-2.5\,\mathrm{TeV}$ and $A\sim 4\,\mathrm{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 $|\mu|\sim 2\,\mathrm{TeV}$ and $|A|\sim 5\,\mathrm{TeV}$. The optimization of the bounce action by CosmoTransitions Wainwright:2011kj , which is used by Vevacious, leads to a slightly stronger and more reliable metastability bound in this region.777As a cross check, forcing Vevacious to use the direct path approximation yields the same lifetimes as our approach. 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\sim 5\,\mathrm{TeV}$ and $\mu\sim 4\,\mathrm{TeV}$ as well as the missing region around $A\sim-5\,\mathrm{TeV}$ and the spikes around $\mu\sim-2.5\,\mathrm{TeV}$ are consequences of the same MDM misidentification as in the Vevacious tree-level 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. 30a we find our tree-level code to be $\sim 5$ times faster than the tree-level and $\sim 200$ times faster than the one-loop Vevacious run.

4.1.3 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\beta$, $X_{t}$ and the stop masses. We therefore expect the Higgs mass to stay close to $125\,\mathrm{GeV}$ when moving away from the point $\boldsymbol{\times}$ along this line.888We have verified using FeynHiggs 2.14.3 Heinemeyer1999 ; Heinemeyer2000 ; Degrassi2003 ; Frank2007 ; Hahn2014 ; Bahl2016 ; Bahl2018 that $124\,\mathrm{GeV}\lesssim m_{h}\lesssim 126\,\mathrm{GeV}$ indeed holds along this line as long as $|\mu|\lesssim 3\,\mathrm{TeV}$. 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 $\mu$ 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 $\mu$.

It can be seen from Fig. 7 that for large negative $\mu$ simultaneous $\tilde{t}$ and $\tilde{\tau}$ vevs (orange) dominate the global minimum for the considered field content until the $\tilde{\tau}$ vevs at these stationary points approach zero around $\mu=-2.2\,\mathrm{TeV}$, and pure $\tilde{t}$ vevs take over. From $\mu\approx-1.8\,\mathrm{TeV}$ onwards the EW vacuum is the global minimum until a CCB vacuum with $\tilde{b}$ vevs appears at $\mu\approx 1.6\,\mathrm{TeV}$. The MDM, on the other hand, is the second deepest $\tilde{b}$-vev minimum for $\mu\lesssim-3.5\,\mathrm{TeV}$, before switching to the $\tilde{t}$-vev minimum, followed by the window of absolute stability $\mu\in[-1.8\,\mathrm{TeV},1.5\,\mathrm{TeV}]$. For positive values of $\mu>1.5\,\mathrm{TeV}$ the instability first develops towards the global $\tilde{b}$-vev minimum until the $\tilde{t}$-vev minimum takes over at $\mu\approx 2\,\mathrm{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 orange $\tilde{t}$ and $\tilde{\tau}$ vevs to the ones with yellow $\tilde{t}$ vevs and red $\tilde{\tau}$ 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 the $\tilde{b}$-vev stationary point approaching the EW vacuum at $\mu\ll 0$ from above. In this regime the $\Delta_{b}$ corrections significantly enhance the bottom Yukawa coupling giving rise to a large mixing in the $\tilde{b}$ sector and a corresponding decrease of one of the $\tilde{b}$ masses. The depth of this stationary point becomes degenerate with the EW vacuum. For even larger negative $\mu$ one $\tilde{b}$ 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. 27 for our field sets Eqs. 30a, 30b and 30c has two accidental $\mathbb{Z}_{2}$ symmetries. The potential is symmetric under simultaneous sign flips of the left- and right-handed sfermions of a kind

[TABLE]

and under simultaneous sign flips of all doublet components

[TABLE]

This results in sets of degenerate and physically equivalent stationary points related by these symmetries.999Since these minima are degenerate they cannot be distinguished in Fig. 7. Since the EW vacuum is also invariant under Eq. 39 the tunnelling time to minima related by this symmetry is always identical. However, since the EW vacuum breaks Eq. 40101010We can, without loss of generality, choose the EW vacuum with $\text{Re}(h_{u}^{0}),\,\text{Re}(h_{d}^{0})>0$. 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. 40 can differ by more than an order of magnitude. This effect has recently been studied for the simpler case of a 2HDM in Branchina2018 .

4.2 Vacuum Stability in the $M_{h}^{125}(\tilde{\tau}{})$ Scenario

A benchmark scenario with light $\tilde{\tau}$ has been proposed in Bahl2018a under the name $M_{h}^{125}(\tilde{\tau}{})$. It is defined by

[TABLE]

The scenario differs from the $M_{h}^{125}$ scenario of Eq. 35 only in greatly reduced soft $\tilde{\tau}$ masses with a correspondingly reduced $A_{\tilde{\tau}{}}$. However, $\mu$ is not reduced and is now $\mu\sim 3m_{\tilde{\tau}{}}$. According to Eq. 38 we would therefore expect vacuum stability constraints to be relevant in the $M_{h}^{125}(\tilde{\tau}{})$ benchmark plane. The authors of Bahl2018a used Vevacious to check for vacuum instabilities in this scenario and found a short-lived region in the parameter space for large $\tan\beta$ and small $m_{A}$.

Our results shown in Fig. 8 confirm these observations. Figure 8 (left) shows a short-lived region for large $\tan\beta$. This region extends towards smaller values of $\tan\beta$ in the low $m_{A}$ regime as noted in Bahl2018a but also in the region of large $m_{A}$. Also visible is the small region of instability for $\tan\beta<1$ noted in Fig. 3. The MDM for the instability at large $\tan\beta$ is a vacuum with $\tilde{\tau}$ vevs as can bee seen from Fig. 8 (centre). Compared to Fig. 3 the absolutely stable region is additionally reduced by a $\tilde{t}$-$\tilde{\tau}$-vev minimum appearing around $m_{A}\sim 1\,\mathrm{TeV}$ and $\tan\beta<30$. The minima with $\tilde{b}$ vevs, which were the MDM for large regions of the $M_{h}^{125}$ scenario, are entirely replaced by minima with $\tilde{\tau}$ vevs. Only a very small purple region with simultaneous $\tilde{b}$ and $\tilde{\tau}$ vevs at the MDM exists. In Fig. 8 (right) the global minimum with $\tilde{b}$ vevs is very similar to Fig. 3 (right). Only at larger $m_{A}$ — where the EW vacuum in the $M_{h}^{125}$ scenario was absolutely stable — global minima with $\tilde{\tau}$ vevs are now present. Our results in the $M_{h}^{125}(\tilde{\tau})$ scenario compared to the $M_{h}^{125}$ 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 $\mu\sim 3m_{\tilde{\tau}{}}$ in this scenario, a further increase of $\mu$ or a decrease of $m_{\tilde{\tau}{}}$ could render the $M_{h}^{125}(\tilde{\tau})$ scenario entirely short-lived. In the region of $\mu\sim 3m_{\tilde{\tau}{}}$ chosen in the $M_{h}^{125}(\tilde{\tau})$ scenario, the vacuum stability constraints show a significant dependence on the parameters $m_{A}$ and $\tan\beta$.

4.3 Vacuum Stability in the $M_{h}^{125}(\mathrm{alignment})$

Scenario

In this section we turn to another scenario from Bahl2018a , the $M_{h}^{125}(\text{alignment})$ scenario, where constraints from vacuum stability turn out to have a very large impact. The scenario is defined by

[TABLE]

with $m_{A}\in[100\,\mathrm{GeV},1\,\mathrm{TeV}]$ and $\tan\beta\in[1,20]$. It is chosen to accommodate light additional Higgs bosons through the so-called alignment without decoupling behaviour Gunion2003 ; Craig2013 ; Carena2014 ; Dev:2014yca ; Carena2015 ; Bernon2015 ; Bernon2016 ; Bechtle:2016kui ; Haber2017 . Alignment without decoupling in the phenomenologically interesting region of relatively small $\tan\beta$ requires $\mu,A_{t}\gg m_{\tilde{t}}$. This requirement—according to Eq. 38—has already been noted to be problematic for vacuum stability in Bahl2018a . Indeed we find that the EW vacuum in this scenario is short-lived through all of its parameter space with instabilities in directions with $\tilde{t}$ 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,

[TABLE]

and vary $A=A_{t}=A_{b}=A_{\tau}$ and $\mu$ 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.3 Heinemeyer1999 ; Heinemeyer2000 ; Degrassi2003 ; Frank2007 ; Hahn2014 ; Bahl2016 ; Bahl2018 . This shows that in order to obtain a long-lived scenario while keeping the correct Higgs mass one could for example change

[TABLE]

Other choices are possible within the theoretical uncertainty of the Higgs mass prediction. The alignment-without-decoupling behaviour implies that the $125\,\mathrm{GeV}$ Higgs boson has SM-like properties. In the MSSM this behaviour arises from a cancellation between tree-level 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. 44 could affect these cancellations and therefore spoil the alignment properties. Indeed, with FeynHiggs we obtain $\sim 10\%$ enhanced couplings of $h$ to bottom and tau pairs leading to a $\sim 10\%$ reduced $\text{BR}(h\to\gamma\gamma)$ compared to the SM for the parameter point of Eq. 44. 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 Aaboud:2018xdt .

In the middle panel of Fig. 9 one can see that the MDM directions through most of the parameter space are directions with $\tilde{t}$ vevs. If the instabilities were in directions with $\tilde{b}$ or $\tilde{\tau}$ 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,\tau}$ while keeping $\mu$ and the $\tilde{t}$ 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 with $\tilde{t}$ 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 $\tilde{\tau}$ vevs through significant parts of the parameter plane, while for large values of $A,\mu\gtrsim 5\,\mathrm{TeV}$ $\tilde{b}$ 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\,\mathrm{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 with $\tilde{t}$ vevs through most of the parameter range, with the exception of the region with small values of $\mu$ where the MDM has $\tilde{\tau}$ vevs. There exist stationary points that are deeper than the MDM with $\tilde{b}$ vevs111111The $\tilde{b}$-vev stationary point is almost degenerate with the deepest $\tilde{t}$-$\tilde{\tau}$-vev stationary point for most of the parameter range and therefore hidden behind the orange points., $\tilde{t}$-$\tilde{\tau}$ vevs and $\tilde{\tau}$ vevs with the $\tilde{\tau}$-vev minimum being the global minimum until the $\tilde{b}$-vev minimum takes over for $\mu\gtrsim 5.5\,\mathrm{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 $\mu$ stationary points with both $\tilde{t}$ and $\tilde{b}$ 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\,\mathrm{TeV}$ would be long-lived if the MDM with $\tilde{t}$ 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.

5 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 Bahl2018a . 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 (namely $\tilde{t}$, $\tilde{b}$ and $\tilde{\tau}$) 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_{h}^{125}$ MSSM benchmark scenario we have found that the impact of vacuum stability constraints is small in the $m_{A}$–$\tan\beta$ parameter plane defining the scenario. However, a variation of the parameters $\mu$ 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 bottom-quark mass and the bottom Yukawa coupling in certain regions of parameter space. These $\Delta_{b}$ corrections can significantly enhance the value of the bottom Yukawa coupling and thus trigger $\tilde{b}$-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 long- and 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}\sim 125\,\mathrm{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 Branchina2018 —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 $\tilde{\tau}$ we have shown important constraints from vacuum stability arising in the $m_{A}$–$\tan\beta$ plane defining the scenario. In Bahl2018a a parameter region with high $\tan\beta$ 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}$.

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 scenarios121212It should be noted that we have omitted the CP-violating scenario that was proposed in Bahl2018a for simplicity only. Our approach is fully applicable also to models containing CP-violating phases. of Bahl2018a 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.

Acknowledgements.

We thank the authors of Vevacious, especially Jose Eliel Camargo-Molina and Werner Porod, for helpful comments and discussions. We also thank Tim Stefaniak and Stefan Liebler for useful discussions. JW gratefully acknowledges funding from the PIER Helmholtz Graduate School.

Appendix 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. 30a, 30b and 30c. The potential is given by

[TABLE]

consisting of a quadratic ($V_{2}$), cubic ($V_{3}$) and quartic ($V_{4}$) part. The quadratic part is given by

[TABLE]

Appendix B Constraints from Vacuum Stability in the Other Benchmark Scenarios of Bahl2018a

In this appendix we present our results for the vacuum stability constraints in the remaining CP-conserving benchmark scenarios defined in Bahl2018a .

The vacuum stability analysis for the $M_{h}^{125}(\text{alignment})$ benchmark plane defined in Eq. 42 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 by $\tilde{t}$ vevs throughout the parameter plane, while the global minimum has $\tilde{b}$ vevs for $\tan\beta\gtrsim 5.5$ and $\tilde{t}$-$\tilde{\tau}$ vevs for smaller values of $\tan\beta$.

The light neutralino scenario $M_{h}^{125}(\tilde{\chi})$ is defined through

[TABLE]

As shown in Fig. 12 the entire benchmark plane features an absolutely stable EW vacuum. This is unsurprising, since it has a very small $\mu$ 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 $\Delta_{b}$ and $\Delta_{\tau}$ corrections.

Finally, the heavy Higgs alignment scenario $M_{H}^{125}$ is defined through

[TABLE]

In this scenario the heavy Higgs boson $H$ has a mass of about $125\,\mathrm{GeV}$ with SM-like couplings. Since the scenario has small soft SUSY-breaking squark mass parameters and a very large $\mu$ 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-vanishing $\tilde{t}$ vevs while the global minimum has $\tilde{t}$-$\tilde{\tau}$ vevs. Considering how constrained the region of phenomenologically viable parameter space for this scenario is (see Fig. 10 of Bahl2018a ) 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.