Vevacious: A Tool For Finding The Global Minima Of One-Loop Effective Potentials With Many Scalars

Several extensions of the Standard Model of particle physics contain additional scalars implying a more complex scalar potential compared to that of the Standard Model. In general these potentials allow for charge and/or color breaking minima besides the desired one with correctly broken SU(2)_L times U(1)_Y . Even if one assumes that a metastable local minimum is realized, one has to ensure that its lifetime exceeds that of our universe. We introduce a new program called Vevacious which takes a generic expression for a one-loop effective potential energy function and finds all the tree-level extrema, which are then used as the starting points for gradient-based minimization of the one-loop effective potential. The tunneling time from a given input vacuum to the deepest minimum, if different from the input vacuum, can be calculated. The parameter points are given as files in the SLHA format (though is not restricted to supersymmetric models), and new model files can be easily generated automatically by the Mathematica package SARAH. This code uses HOM4PS2 to find all the minima of the tree-level potential, PyMinuit to follow gradients to the minima of the one-loop potential, and CosmoTransitions to calculate tunneling times.


I. INTRODUCTION
A major part of the phenomenology of the incredibly successful standard model of particle physics (SM) is the spontaneous breaking of some (but not all) of the gauge symmetries of the Lagrangian density by the vacuum expectation value (VEV) of a scalar field charged under a subgroup of the SM gauge group. The entire scalar sector of the SM consists of a doublet of SU (2) L which also has a hypercharge under U (1) Y equal in magnitude to that of the lepton SU (2) L doublet. The potential energy of the vacuum is minimized by the scalar field taking a constant non-zero value everywhere. The presence of this VEV radically changes the phenomenology of the theory, and allows for masses for particles that would be forced to be massless if the gauge symmetries of the Lagrangian density were also symmetries of the vacuum state.
Since this scalar field is the only field in the SM that can possibly have a non-zero VEV while preserving Lorentz invariance, finding the minima of the potential energy is straightforward, though of course evaluating it to the accuracy required is quite involved [1][2][3].
Also, with the current measurements for the masses of the top quark and Higgs boson, one finds that the SM potential at one-loop order is actually unbounded from below for a fixed value of the renormalization scale. The value of the Higgs field for which the potential is lower than the desired vacuum is so high that one may worry that large logarithms of the Higgs field over the electroweak scale would render the loop expansion non-convergent. However, the effect of large logarithms can be resummed, and the conclusion that our vacuum is only metastable persists using the renormalization-group-improved effective potential [1][2][3][4][5].
The existence of multiple non-equivalent vacua both raises technical challenges and introduces interesting physics. The technical challenges are now that one has to find several minima and evaluate which is the deepest, as well as calculate the tunneling time from a false vacuum to the true vacuum. However, this is an important ingredient in theories where a first-order phase transition explains the baryon asymmetry of the universe through the sphalerons occuring in the nucleation of bubbles of true vacuum (see [6] and references therein).
Many extensions of the SM introduce extra scalar fields. Sometimes these fields are introduced explicitly to spontaneously break an extended gauge symmetry down to the SM gauge group [7,8], and they are assumed to have non-zero VEVs at the true vacuum of the theory. Other times they are introduced for other reasons, such as supersymmetry [9], and often non-zero VEVs for such fields would be disastrous, such as breaking SU (3) c and/or U (1) EM , which excludes certain parts of the parameter space of the minimal supersymmetric standard model (MSSM) from being phenomenologically relevant.
The technical challenges are much tougher when multiple scalar fields are involved. Even a tree-level analysis involves solving a set of coupled cubic equations, the so-called minimization or tadpole equations. It has generally only been attempted for highly symmetric systems such as two Higgs doublet models (2HDM) [10,11] or with only a minimal amount of extra degrees of freedom such as the (assumed) three non-zero VEVs of the next-to-minimal supersymmetric standard model (NMSSM) [12][13][14].
Since a general solution is usually too difficult, the question of the stability of VEV configurations against tunneling to other minima of the potential is often ignored. Instead, potentials are often engineered to have a local minimum at a desired VEV configuration through ensuring that the tadpole equations are satisfied for this set of VEVs. This approach allows one to go beyond tree level straightforwardly, and one-loop tadpoles are the norm, and in supersymmetric models two-loop contributions are often included [15]. This local minimum is implicitly assumed to be stable or long-lived enough to be physically relevant. Unfortunately, as some examples will show, local minima which are not the global minimum of their parameter point are often extremely short-lived, excluding some benchmark parameter points for some models.
The program Vevacious has been written to address this. Given a set of tadpole equations and the terms needed to construct the one-loop effective potential 1 , first all the extrema of the tree-level potential are found using homotopy continuation (HOM4PS2), which are then used as starting points for gradient-based minimization (PyMinuit) of the (real part of the) one-loop potential, and finally, if requested, the tunneling time from an input minimum to the deepest minimum found is estimated at the one-loop level (CosmoTransitions). The program is intended to be suitable for parameter scans, taking parameter points in the SLHA format [16,17] and giving a result within seconds, depending on the number of fields allowed to have non-zero VEVs and the accuracy of the tunneling time required.

II. THE POTENTIAL ENERGY FUNCTION AT TREE LEVEL AND ONE-LOOP LEVEL
The terminology of minimizing the effective potential of a quantum field theory is rather loaded. Hence first we shall clarify some terms and conventions that will be used in the rest of this article. In the following we consider models where only scalars can get a VEV as required by Lorentz invariance.
In principle, the effective potential is a real-valued 2 functional over all the quantum fields of the model. However, under the assumption that the vacuum is homogeneous and isotropic, for the purposes of determining the vacua of the model, the effective potential can be treated as a function of sets of (dimensionful) numbers, which we shall refer to as field configurations.
Each field configuration is a set of variables which correspond to the classical expectation values for the spin-zero fields which are constant with respect to the spatial co-ordinates.
The example of the SM is relatively simple: the field configuration is simply a set of two complex numbers, which are the values of the neutral and charged scalar fields assuming that each is constant over all space. These four real degrees of freedom can be reduced to a single degree of freedom by employing global SU (2) L and phase rotations, leaving an effective potential that is effectively a function of a single variable.
Henceforth we shall assume that each complex field is treated as a pair of real degrees of freedom, and note that this may obscure continuous sets of physically equivalent degrees of freedom which are manifestly related by phase rotations when expressed with complex fields.
Also we shall refer to the local minima of the effective potential as its vacua, and label the global minimum as the true vacuum, while all the others are false vacua. A potential may have multiple true vacua, either as a continuous set of minima related by gauge transformations as in the SM for example, or a set of disjoint, physically inequivalent minima, each of which may of course be a continuous set of physically equivalent minima themselves. In cases where there is a continuous set of physically equivalent minima, we assume that a single exemplar is taken from the set for the purposes of comparison of physically inequivalent minima. Furthermore, the term vacuum expectation value can be used in many confusing ways.
In this work, VEVs will only refer to the sets of constant values which the scalar fields have at the field configurations which minimize the effective potential. Hence we do not consider the effective potential to be a "function of the VEVs", rather a function of a set of numbers that we call a field configuration.

A. The tree-level potential and tadpole equations
If one merely considers the SM at tree level, minimizing the potential is straightforward.
A global SU (2) L rotation can bring the part of the potential due to the scalar doublet into the form where φ is the neutral component of the SU (2) L doublet. After a little differentiation and algebra, one finds that if λ > 0, µ 2 < 0, then the potential is minimized for |φ| = v = −µ 2 /λ.
However, for a set of N real scalar degrees of freedom φ i (writing complex scalar fields as two separate real scalars), the scalar part of the tree-level potential of a renormalizable quantum field theory in four space-time dimensions is of the form which, when differentiated with respect to the N independent φ i , yields N polynomial equations of up to degree three. We have assumed that any terms linear in the fields have been removed by constant shifts of the fields.
Although we assume renormalized potentials here for simplicity, the methods used by Vevacious are equally applicable to non-renormalizable potentials, as long as V tree is expressed as a finite-degree polynomial. The value of loop corrections to a non-renormalizable potential may be debatable, but Vevacious can be restricted to using just the tree-level potential.
While closed-form solutions for cubic polynomials in one variable exist, solving a coupled system in general requires very involved algorithms, such as using Gröbner bases to decompose the system [20,21], or homotopy continuation to trace known solutions of simple systems as they are deformed to the complicated target system of tadpole equations.

The homotopy continuation method
The homotopy continuation method [22,23] has found use in several areas of physics [24][25][26], in particular to find string theory vacua [27,28] and extrema of extended Higgs sectors [29], where the authors investigated a system of two Higgs doublets with up to five singlet scalars in a general tree-level potential, and [30], where systems of up to ten fields were allowed to have non-zero VEVs. In contrast, the Gröbner basis method is deemed prohibitively computationally expensive for systems involving more than a few degrees of freedom [20].
The numerical polyhedral homotopy continuation method is a powerful way to find all the roots of a system of polynomial equations quickly [31]. Essentially it works by continuously deforming a simple system of polynomial equations with known roots, with as many roots as the classical Bézout bound of the system that is to be solved (i.e. the maximum number of roots it could have). The simple system with known roots is continuously deformed into the target system, with the position of the roots updated with each step. While the method is described in detail in [22,23], a light introduction can be found for example in [29].

B. The one-loop potential
The general form of the renormalized one-loop effective potential [32,33] is where V tree is as above and V counter has the same form as V tree , i.e. a polynomial of the same degree in the same fields, but the coefficients are instead the renormalization-dependent finite parts of the appropriate counterterms [32]. The term V mass has the form, for a given field configuration Φ, where the sum runs over all real scalar, Weyl fermion, and vector degrees of freedom, with s n being the spin of the degree of freedom. Complex scalars and Dirac fermions are accounted for as mass-degenerate pairs of real scalars and Weyl fermions respectively.
For scalar degrees of freedom, theM 2 n (Φ) are the eigenvalues of the second functional derivative of V tree , i.e. the eigenvalues of (M 2 Thus theseM 2 n (Φ) are the eigenvalues of the tree-level scalar "mass-squared matrix" that would be read off the Lagrangian with the scalars written as fluctuations around the field configuration. Unless the field configuration corresponds to a minimum of the effective potential, these do not correspond to physical masses in any way, of course.
Likewise, theM 2 n (Φ) for fermionic and vector degrees of freedom are the eigenvalues of the respective "mass-squared matrices" where the scalar fields are taken to have constant values given by the field configuration. (The fermion mass-squared matrix is given by the mass matrix multiplied by its Hermitian conjugate.) The terms c n depend on the regularization scheme. In the MS scheme, c n is 3/2 for scalars and Weyl fermions, but 5/6 for vectors, while in the DR scheme [33,34], more suitable for supersymmetric models, c n is 3/2 for all degrees of freedom. Since this is a finite-order truncation of the expression, the renormalization scale Q also appears explcitly in the logarithm, as well as implicitly in the scale dependence of the renormalized Lagrangian parameters.
Much of the literature on one-loop potentials (including [33]) assumes a renormalization scheme where V counter is zero; however, such a scheme is often inconvenient for other purposes, such as ensuring tadpole equations have a given solution at the one-loop level, see e.g. the appendix of [35]. (SARAH automatically generalizes this approach to extended SUSY models as explained in [36,37].) Finally, we also note that models can be constructed where spontaneous symmetry breaking does not happen at tree level, but does exist when one takes loop corrections into account [18,38].

Scale dependence
As noted, the one-loop effective potential depends on the renormalization scale. Ideally one would use the "renormalization-group improved" expression for the potential [32] as this is invariant under changes of scale; however, this is often totally impractical except for potentials with only a handful of parameters and a single scalar field.
If one must use a scale-dependent expression for the potential, as is often the case, the renormalization scale should be chosen carefully: if one chooses a scale too high or too low, one may find that with a finite-loop-order (and thus scale-dependent) effective potential, there is no spontaneous breaking of any symmetry, or even that the potential is not bounded from below [39]! This is often simply due to the fact that higher orders become more important in such a case, especially when corrections from the next order would introduce new, large couplings, such as often happens when going from tree level to one loop. It can also be that the scale is so large or small that the loop expansion is no longer a reliable expansion. We also note that rather undermines arguments that radiative effects do not change tree-level conclusions on the absolute stability of vacua such as in [40] (where the argument also fails to take into account that there may not even be a scale at which the renormalization condition used can be satisfied).
Indeed, it is crucial that the scale is chosen so that the loop expansion is valid. Explicitly, large logarithms should not spoil the perturbativity of the expansion in couplings. Loops with a particle n typically come with a factor of ln(M 2 n (Φ)/Q 2 ) along with the factor of α n = [relevant coupling] 2 /(4π), and thus α n ln(M 2 n (Φ)/Q 2 ) should remain sufficiently smaller than one such that the expansion can be trusted [32]. A rough first estimate then of the region of validity, assuming that the dimensionful Lagrangian parameters are all of the order of the renormalization scale to some power, is where ln(v 2 /Q 2 )/(4π) ≤ 1/2, say, for a field configuration with vector length v, so where the VEVs are within a factor of e π 20 of the renormalization scale.
Furthermore, it is in general not valid to compare the potential for different field configurations using a different scale for each configuration, if one is using a scale-dependent effective potential. The reason is that there is an important contribution to the potential that is field-independent yet still depends on the scale 3 [41,42]. Of course, if one knows the full scale dependence of all the terms of the Lagrangian regardless of whether they lead to field-dependent contributions to the effective potential, then one can correctly evaluate different field configurations at different scales.

C. Gauge dependence
The one-loop potential is explicitly gauge-dependent [43,44]. However, as shown in [44,45], the values it takes at its extrema are independent of the gauge chosen, except for spurious extrema of poorly-chosen gauges. The popular R ξ gauges are well-behaved and do not have fictitious gauge-dependent extrema for reasonable choices of ξ [45].
It is also possible to formulate the effective potential in terms of gauge-invariant composite fields [46], though this may not always be practical. One can also verify the gaugeindependence of extrema using more complicated gauges and applying BRST invariance [47].

D. Convexity
As shown in [18], the effective potential, which can be thought of as the quantum analogue of the classical potential energy for constant fields, is real for all values of the fields. However, the loop expansion leads to complex values in regions where the classical potential is nonconvex. While one can take the convex hull of the truncated expansion of the potential when evaluating the potential for configurations of fields in the convex region, it is not particularly helpful for the purpose of computing tunneling transition times. Fortunately, the one-loop truncation of the effective potential as a function of constant values for the fields can consistently be interpreted as a complex number with real part giving the expectation value of the potential energy density for the given field configuration and imaginary part proportional to the decay rate per unit volume of this configuration [19].

E. Comparing two vacua
If there are two or more physically inequivalent minima of a potential, then it is vitally important to know if the phenomenologically desired minimum is the global minimum, or, if not, how long the expected tunneling time to the true vacuum is.
Given the issues raised above in sections II B 1 and II C, it is safe to use a scale-and gaugedependent one-loop effective potential to compare two inequivalent minima provided that the scale is held fixed and that the two minima are within the region of validity determined by the renormalization scale. Of course, an explicitly gauge-and scale-independent expression for the effective potential would obviously be unburdened by such concerns, but unfortunately it is rare to be able to formulate such an expression.

III. TUNNELING FROM FALSE VACUA TO TRUE VACUA
The usual expression for the decay rate Γ per unit volume for a false vacuum is given in [48,49] as where A is a factor which depends on eigenvalues of a functional determinant and B is the bounce action. The A factor is typically estimated on dimensional grounds as it is very complicated to calculate and, because of the exponentiation of B, is far less important than getting the bounce action as accurate as possible. If A is taken to be O((100 − 1000 GeV) 4 ), then for Γ/vol. to be roughly the age of the known Universe to the fourth power, B must be around 400 , and a per-cent variation in B leads to a factor of e variation in the tunneling time.
Given a path through the field configuration space from one vacuum to another with a lower value, which for convenience we shall label as the false vacuum and true vacuum respectively, one can solve the equations of motion for a bubble of true vacuum of critical size in an infinite volume of false vacuum [48,50]. This allows one to calculate the bounce action and thus the major part of the tunneling time.
Unfortunately, this means that to calculate the tunneling time from a false vacuum to a true vacuum, one needs to evaluate the potential along a continuous path through the field configuration space, and even though the extrema of the potential are gauge-invariant as noted above, the paths between them are not. However, it has been proved that at zero temperature, the gauge dependence at one-loop order cancels out [51].
At finite temperature, the situation is not so clear, though the Landau gauge may be most appropriate [52]. While some studies have shown that for "reasonable" choices of gauge, the differences in finite-temperature tunneling times are small [53], it is still possible to choose poor gauges that can even obscure the possibility of tunneling [54].

A. Objectives
The Vevacious program is intended as a tool to quickly evaluate whether a parameter point with a given set of VEVs, referred to henceforth as the input vacuum, has, to one-loop order 4 , any vacua with lower potential energy than the input vacuum, and, optionally, to estimate the tunneling time from the input vacuum to the true vacuum if so.
A typical use envisaged is a parameter scan for a single model. Some effort needs to be put into creating the model file in the first place, though this is straightforward if using

B. Outline
Here we present the steps taken by Vevacious, which are schematically shown in Fig. 1.
(1) An input file in the SLHA format [16] is read in to obtain the Lagrangian parameters defining the parameter point, required to evaluate the potential. We emphasize that even though the SUSY Les Houches Accord is used as the format, the model itself does not need to be supersymmetric, as long as the SLHA file contains appropriate BLOCKs.
(2) All the extrema of the tree-level potential are found using the homotopy continuation (3) The tree-level extrema are used as starting points for gradient-based minimization of the one-loop effective potential. The MINUIT algorithms [55] are used here through the Python wrapper PyMinuit [56]. The points where PyMinuit stops are checked to see if they are really minima 5 . Any saddle point is then split into two further points, displaced from the original in the directions of steepest descent by amounts given by the <saddle nudges> arguments (see sec. VI C), which are then also used as starting points for PyMinuit. By default, PyMinuit is restricted to a hypercube of field configurations where each field is only allowed to have a magnitude less than or equal to one hundred times the renormalization scale of the SLHA file. This is rather excessive by the reasoning of sec. II B 1, which would lead one to take at most maybe ten or twenty times the scale as an upper limit; however, it was considered better to allow the user to decide whether the results of Vevacious are within a trustworthy region. b. Rolls to one-loop minima. Vevacious rolls from the tree-level extrema to the minima of the one-loop effective potential before comparing them, because in general the VEVs get shifted. In addition, extrema that change their nature with radiative corrections, such as the field configuration with zero values for all the fields in the Coleman-Weinberg model of radiative spontaneous symmetry breaking [38], which is a minimum of the tree-level potential but a local maximum of the one-loop effective potential, are found.

c. Calculates tunneling times or upper bounds on them. A parameter point in a model
is not necessarily ruled out on the basis that the desired minimum of the potential is not the global minimum, since a metastable configuration with a lifetime of roughly the observed age of the Universe or longer is compatible with the single data point that we have. to change the A factor to (246 GeV) 4 for example, if one feels that the electroweak scale is a more appropriate choice.
D. Limitations f. Garbage in, garbage out. Vevacious performs very few sanity checks, so rarely protects the user from their own mistakes. For instance, Vevacious does not check if the potential is bounded from below. However, there are some checks, such as those which result in the warning that the given input minimum was actually rather far from the nearest minimum found by Vevacious (though Vevacious carries on regardless after issuing the warning). Another important sanity check that is not performed is to check that the SLHA BLOCKs required by the model file are actually present in the given SLHA input file. The user is fully responsible for providing a valid SLHA file to match the model. As model files are expected to be produced automatically by software such as SARAH, it is expected that the SLHA files for the model will also be prepared consistently with the expected BLOCKs.
Unfortunately it is quite easy to miss this point when using the example model files provided by default with Vevacious: if one does use these files, one must use the correct model file for the input SLHA files, as described in sec. VII. to limit VEVs to be no larger than twenty times the renormalization scale. One could also insert more complicated Python code here, but one should be aware that the PyMinuit object deals with a potential where the field values are in units of the renormalization scale.
h. Not guaranteed to find minima induced purely by radiative effects. While Vevacious does find all the extrema of the tree-level effective potential, there is no guarantee that these correspond to all the minima of the effective potential at the one-loop level. The strategy adopted by Vevacious will find all the minima of the one-loop effective potential that are in some sense "downhill" from tree-level extrema, but any minima that develop which would require "going uphill" from every tree-level extremum will not be found. Vevacious should be used bearing in mind the caveat that it will not find vacua with nonzero VEVs for degrees of freedom which are not allowed non-zero VEVs in the model file.
It is up to the user to decide on the best compromise between speed and comprehensiveness by choosing which degrees of freedom to use.
Importantly, the user is responsible for ensuring that the model file has a tree-level potential which has a discrete number of minima. Mostly this means that the user has to identify unphysical phases and hence remove the associated degrees of freedom from the imaginary parts of such complex fields. Another way that problems can arise is when there are flat directions such as the tan β = 1 direction of the MSSM with unbroken supersymmetry, even though the degeneracy of these directions may be lifted by loop corrections. to save giving several command-line arguments, and this file is described along with these arguments in sec. VI C.
A. Preparing the input file for SARAH SARAH [57][58][59][60] is a tool to derive many analytical properties of a particle physics model, like mass matrices, tadpole equations, vertices and renormalization group equations, from a very short user input. This information can be used, for instance, to write model files for several matrix generators or source code for SPheno [61,62]. While previous versions of SARAH were optimized for supersymmetric models but supported also to some extent nonsupersymmetric models, SARAH 4 will provide a simplified input and new features also for non-supersymmetric models [63]. In addition, SARAH 4 supports also the output of input files for Vevacious. To get this file, run in Mathematica The possible options are: • ComplexParameters, Value: list of parameters, Default: {}: By default, all parameters are assumed to be real when writing the Vevacious input files. However, the user can define those parameters which should be treated as complex.
• IgnoreParameters, Value: list of parameters, Default: {}: The user can define a list of parameters which should be set to zero when writing the Vevacious input.
• OutputFile, Value: String, Default MyModel.vin, where MyModel here is the same name as is given in Start["MyModel"]; above: The name used for the output file.
The first two options allow one to treat parameters differently in the Vevacious output as defined in the SARAH model file. It may be in the user's interest to try to speed up the evaluation by taking out those parameters which on physical grounds play only a subdominant role, but the gain by doing so has yet to be quantified.
Example: MSSM with stau VEVs Here we discuss briefly the main steps to prepare an MSSM version including stau VEVs. For a more general discussion of the format of the SARAH model files, we refer the interested reader to [60,64]. In general, three changes are always necessary to include new VEVs in a model.  where it is important that the VEVs have names that are at least two characters long. The first two lines are the standard decomposition of the complex Higgs scalar There are two new features in SARAH 4 which are shown here: (i) it is possible to give VEVs just to specific generations of a field, such as in this example we only use the third one (vLR [3]) -in the same way, one can allow for smuon and stau VEVs using vLR [2,3]; (ii) VEVs can have real and imaginary parts. Until now, complex VEVs in SARAH had been defined by an absolute value and a phase (ve iφ ); however, it is easier to handle within Vevacious in the form v R + iv I . within SLHA files, and should be ignored by programs that do not recognize them. Vevacious is intended to be used for many other models, so accepts any BLOCKs that are mentioned in its model file and looks for them in the given SLHA file. In this sense, the user is free to define BLOCKs as long as the names are unbroken strings of alphanumeric characters (e.g. BLOCK THISISAVALIDNAME or BLOCK MY BLOCK 0123 6 ). However, we strongly advise against redefining those BLOCKs specified by [16] and [17], or any of their elements, to have meanings other than those given in [16] and [17]. All three can be derived from elements 2, 3, and 4 of HMIX, but are much more suited to conversion between renormalization schemes and different gauges (as v d and v u , and thus tan β, are gauge-dependent quantities, with values which also depend on the renormalization conditions).

Scale dependence in the SLHA file
Vevacious is a tool for finding minima for a one-loop effective potential evaluated at a single scale, and thus it is important that the Lagrangian parameters are provided con-sistently at this scale. The scale is also required to be given explicitly. Since the SLHA convention specifies that running parameters are given in BLOCKs each with their own scale, at first glance this may seem problematic. Even worse, the format allows for multiple instances of the same BLOCK, each with its own scale. However, the default behaviour of SPheno, SoftSUSY, SuSpect, and ISAJET when writing SLHA output is to give all running parameters consistently at a single scale. This is the behaviour that Vevacious assumes.
The explicit value of the scale Q used in V mass in eq. (4) is that given by the BLOCK GAUGE. All other BLOCKs are assumed to be at this same scale Although it is not the default behaviour of any of the popular spectrum generators to give the same BLOCKs at different scales, if Vevacious finds multiple instances of the same BLOCK, the BLOCK with the lowest scale is used and the others ignored. In addition, Vevacious performs a consistency check that all the BLOCKs used have the same scale, aborting the calculation if not.

SLHA expression of parameters at different loop orders
The output BLOCKs enumerated in the SLHA papers are specified to be in the DR renormalization scheme 7 , but some users may prefer a different renormalization scheme. The SLHA does not insist on private BLOCKs adhering to the same standards of those explicitly part of the accord, so Vevacious allows for a certain pattern of private BLOCKs to give values for a different renormalization scheme. (Again, we strongly advise against using the BLOCKs explicitly mentioned in [16,17] to convey values that do not adhere to the definitions in [16,17].) The additional renormalization schemes that Vevacious allows are those where the finite parts of Lagrangian parameters are themselves apportioned into loop expansions, e.g. µ + δµ, where δµ is considered to be a parameter already of at least one order higher than µ.
To allow for different renormalization conditions of this type, Vevacious first looks for extra ("private") SLHA BLOCKs that specify particular loop orders. Since Vevacious deals with one-loop effective potentials, it has two categories of parameters: "tree-level" and "one-loop". When writing the minimization conditions for the tree-level potential, it uses exclusively the "tree-level" values. When writing the full one-loop effective potential, it uses both sets appropriately to avoid including spurious two-loop terms. With reference to eqs. (2), (3), and (4), Vevacious writes combines the sum V tree + V counter by inserting the "one-loop" parameter values into V tree as part of V 1-loop . (The tree-level potential function that is also written automatically for convenience, as mentioned in sec. IV C, uses the "treelevel" values, of course.) The term V mass is already a loop correction, so "tree-level" values are used in theM 2 n (Φ) functions. If only a single value for any parameter is given, it is assumed to be in a scheme where it has a single value which is to be used in all parts of the effective potential.
As an example, within the renormalization used by SPheno3.1.12, at the point SPS1a , µ has the value 374.9 GeV at tree level, and 394.4 GeV at one loop. Vevacious inserts the value 374.9 for µ in the minimization conditions (as the units are assumed to be in GeV, as per the SLHA standard), and also into the "mass-squared" matrices that are part of the evaluation of the V mass contributions to the one-loop effective potential. Vevacious inserts the value 394.4 for µ in the polynomial part of the potential, accounting for the contributions of both V tree and V counter together.
In detail, when inserting a "tree-level" value, Vevacious looks for an SLHA BLOCK with the prefix "TREE" first, and uses that value as its "tree-level" value. If there is no BLOCK with that prefix, or if there is such a BLOCK but it does not specify the appropriate element, then the BLOCK without the prefix is assumed to have the appropriate value. Likewise, when inserting a "one-loop" value, BLOCKs with the prefix "LOOP" have priority. Hence for the SPS1a example above, Vevacious could be given an SLHA file with 394.4 as element 1 of HMIX 8 and 374.9 as element 1 of TREEHMIX. When Vevacious is writing the "tree-level" value of µ, it would first look for element 1 of TREEHMIX, and since it would find 374.9, this value would be used, and element 1 of HMIX would not be looked for. When writing the "one-loop" value, element 1 of LOOPHMIX would be looked for, but not found, and because of this, element 1 of HMIX would then be looked for, and 394.4 would be found and used. One can specify element 1 of both TREEHMIX and LOOPHMIX, and then Vevacious would never use element 1 of HMIX, which could give the DR value, for example, without worrying that it might mix schemes in the calculation. (If one prefers to use other prefixes, both "TREE" and "LOOP" can be replaced by other strings in the model file in the <block prefixes> element; e.g. <block prefixes tree="ZEROLOOP" loop="ONELOOP" / >, so that ZEROLOOPHMIX and ONELOOPHMIX would be looked for appropriately).
Those aspects are taken into account in the SPheno [61,62] output of SARAH 4 [63]. For this purpose, a new flag has been introduced which can be used in the SLHA input file In that case, the new tree and one-loop level block will be present.

C. Setting up and running Vevacious
There are many options that can be passed to Vevacious. If values other than the defaults are required, they can either be passed by command-line arguments or with an initialization file in XML format. If an option is given by both commmand-line argument and in an initialization file, the command-line argument takes precedence.
We enumerate the options as they would be as command-line arguments setting the options to the defaults. All floating point numbers may be given as standard decimals, such as 0.1234, or in scientific E notation, such as 1.234E-1 (uppercase 'E' or lowercase 'e', with or without preceeding '0' characters, with or without '+' in the exponent, e.g. 987 or 9.874E+002 or 0098.7e001).
1. --hom4ps2 dir=./HOM4PS2/ This is a string giving the path to the directory where the HOM4PS2 executable is.

--homotopy type=1
This is an integer used to decide which mode HOM4PS2 should run in: 1 is used for polyhedral homotopy and 2 for linear homotopy.  (where some linebreaks have been inserted so that the output fits the width of the page). We can see that the numerical minimization did not quite roll properly to a zero VEV forτ L at the input minimum, but stopped extremely close to it. Here, an upper limit of the lifetime has been calculated using a direct path between the input and the global minimum. The same information is also written to the SLHA file using the new block VEVACIOUSRESULTS,   Vevacious first deletes any files with the same name as was given as the output file, and if Vevacious was unable to end properly, no output is produced. In the case that there were problems during the run which did not crash the program, Vevacious prints warnings, and creates an additional XML element <warning> in the results file, and also appends these warnings in a BLOCKVEVACIOUSWARNINGS. Possible warnings are (with <...> standing for strings describing field configurations, potential depths, or the string given by PyMinuit when it throws an exception) • No tree-level extrema were found. This might happen for example if HOM4PS2 did not find any real solutions (recalling that complex fields have already been written as pairs of real scalars) because it was given a system with continuous degenerate solutionshowever, it is not necessarily the case that this is indicative of this problem, and also it is not necessarily guaranteed to result from such problematic systems.
No t r e e −l e v e l extrema were found .
• PyMinuit threw exceptions: To our knowledge, there is only one publicly-available program that purports to find the global minimum of a potential of a quantum field theory: ScannerS [66]. However, at the time of writing, the routines to actually find the global minimum are still under development and are not available.
The model files generated by SARAH use the internal expressions that have already been cross-checked in [59,64]. In addition, the potential was always found to have a local minimum at the input field configuration, as expected, for a wide range of test points. Furthermore, using the model file SARAH-SPhenoNMSSM JustNormalHiggsAndSingletVevs.vin described below and with the modification to Vevacious to use only a tree-level analysis as described in sec. IV C, we confirmed the results of [20].
Several example model and parameter files are provided with the download of Vevacious.
Three model files are given for the MSSM with different allowed non-zero VEVs: • SARAH-SPhenoMSSM JustNormalHiggsVevs.vin with only the normal real neutral Higgs VEVs allowed; • SARAH-SPhenoMSSM RealHiggsAndStauVevs.vin with real VEVs allowed for the neutral Higgs components and for the staus; • SARAH-SPhenoMSSM RealHiggsAndStauAndStopVevs.vin with real VEVs allowed for the neutral Higgs components, for the staus, and for the stops.
These model files were generated automatically with SARAH 4. They assume that the SLHA Three model files are also provided for a constrained version of the NMSSM (where the dimensionful superpotential parameters µ, µ , and ξ F , and the soft SUSY-breaking parameters m 2 3 , m 2 S , and ξ S , in the notation of [17], are set to zero) with different allowed non-zero VEVs: • SARAH-SPhenoNMSSM JustNormalHiggsAndSingletVevs.vin with only the normal real neutral Higgs and singlet VEVs allowed; • SARAH-SPhenoNMSSM RealHiggsAndSingletAndStauVevs.vin with real VEVs allowed for the neutral Higgs components, for the singlet, and for the staus; • SARAH-SPhenoNMSSM RealHiggsAndSingletAndStauAndStopVevs.vin with real VEVs allowed for the neutral Higgs components, for the singlet, for the staus, and for the stops. In addition, several example parameter points have been provided to demonstrate the existence of metastable points: • SPS1a, the CMSSM point SPS1a from [65] (which is stable); • CMSSM CCB, which corresponds to the CMSSM best-fit point including LHC and m h = 126 GeV constraints from [67] (which has a global charge-and color-breaking minimum); • NUHM1 CCB, which corresponds to the NUHM1 best-fit point ("low") from [68] (which also has a global charge-and color-breaking minimum); • CNMSSM wrong neutral, which corresponds to benchmark point P1 from [69] (which also has a neutral global minimum which however is not the input minimum).
These parameter files are given in the SLHA output of SPheno as SPS1a.slha.out and so on, with example output as SPS1a.vout and so on. (The SLHA input files are SPS1a.slha.in and so on.) We note that the CCB vacua found for CMSSM CCB, for example, have VEVs of the order of five times the renormalization scale (shown in sec. VI D). This should not cause much concern, as ln(5 2 )/(4π) 0.256 so the one-loop effective potential should still be reasonable around this minimum. However, even if one restricts the VEVs to be less than twice the scale, as described in sec. IV D, Vevacious still finds a tunneling time (upper bound) of 10 −6 times the age of the known Universe to a CCB configuration at the edge of the bounding hypercube.
The MSSM model and parameter files are in the MSSM subdirectory and those of the NMSSM are in the NMSSM subdirectory. To tackle this problem, we have presented the program Vevacious as a tool to quickly evaluate one-loop effective potentials of a given model. It finds all extrema at tree level, which allows for a first check for undesired minima. Starting from these extrema, it calculates the one-loop effective potential to obtain a more reliable result. This is important as loop contributions can potentially change the nature of an extremum. In the case that the original minimum turns out to be merely a local minimum rather than the global minimum, the possibility is given to evaluate the tunneling time. As test cases we have considered supersymmetric models, but the program can be used for a general model provided that the tree-level potential is polynomial in the scalar fields. If you want to create new input files using SARAH, you need at least Mathematica 7. To compile SPheno to create the numerical input for Vevacious, a Fortran compiler such as gfortran or ifort is needed.
Vevacious makes use of several public tools. We give here only a brief introduction to the installation of these tools but refer to the corresponding references and authors for more information.
We assume in this descrption that all codes are downloaded and extracted in the same directory. The placeholder for the path to this directoy is called $VPATH in the following.      MyModel can be for instance MSSM or NMSSM.  Vevacious is not restricted to a specific model, but any necessary information to evaluate the one-loop effective potential in a given model, assuming a particular set of VEVs, must be provided by the user as a model file. This file can usually be produced by SARAH as described in sec. VI A. However, in the case that the user wants to modify things manually, we give here more information about the format. The XML element <input vevs> serves the dual purpose of enumerating the allowed non-zero VEVs along with specifying what is considered to be the input minimum. In the MSSM example here, the Higgs fields of course are allowed VEVs, and the input minimum is specified by vdR (v d ) being provided in the SLHA file in the BLOCKHMIX as entry 102, and vuR (v u ) by entry 103. The model file is also allowing non-zero stau VEVs (vLR3, vER3), and also at the same time specifying that their value at the input minimum is 0. (One can actually put valid Python code here within the quote marks for the input minimum values, though it is not recommended, and surrounding the code with brackets is advised if one insists on going through with it.)

SPheno
In addition, the is usually some redundancy in the sets of VEVs which minimize the tree-level scalar potential, since different minima can be related by phase rotations.
To reduce the number of redundant solutions, it can be explicitly defined that some VEVs have to be positive using <taken positive> ... </taken positive>. To calculate the one-loop effective potential all field-configuration-dependent "masses" M 2 n (Φ) have to be specified. This happens by using for each mass-squared matrix the XML element <mass-squared matrix ...>. Two attributes must be given: an overall constant (factor) which takes into account the degrees of freedom of the particle, including the spin s n from eq. 4, and also saves reproducing identical matrices due to color factors or pairs of charge conjugates, for example. This factor is given by r · c F · (−1) s (2s + 1), where r = 1 holds for real bosons or Majorana fermions, and r = 2 for complex bosons and Dirac fermions, and c F is the number of degenerate states (e.g. 6 = 3 colors ×2 charge-conjugate states for quarks, if SU (3) c is unbroken in the model file).
The body of each block contains all entries of the mass-squared matrix. Here, the