Decay Rate of Electroweak Vacuum in the Standard Model and Beyond

We perform a precise calculation of the decay rate of the electroweak vacuum in the standard model as well as in models beyond the standard model. We use a recently-developed technique to calculate the decay rate of a false vacuum, which provides a gauge invariant calculation of the decay rate at the one-loop level. We give a prescription to take into account the zero modes in association with translational, dilatational, and gauge symmetries. We calculate the decay rate per unit volume, $\gamma$, by using an analytic formula. The decay rate of the electroweak vacuum in the standard model is estimated to be $\log_{10}\gamma\times{\rm Gyr~Gpc^3} = -582^{+40~+184~+144~+2}_{-45~-329~-218~-1}$, where the 1st, 2nd, 3rd, and 4th errors are due to the uncertainties of the Higgs mass, the top quark mass, the strong coupling constant and the choice of the renormalization scale, respectively. The analytic formula of the decay rate, as well as its fitting formula given in this paper, is also applicable to models that exhibit a classical scale invariance at a high energy scale. As an example, we consider extra fermions that couple to the standard model Higgs boson, and discuss their effects on the decay rate of the electroweak vacuum.


Introduction
In the standard model (SM) of particle physics, it has been known that the Higgs quartic coupling may become negative at a high scale through quantum corrections, so that the Higgs potential develops a deeper vacuum. The detailed shape of the Higgs potential depends on the Higgs and the top masses; with the recently observed Higgs mass of ∼ 125 GeV, it has been known that the electroweak (EW) vacuum is not absolutely stable if the SM is valid up to ∼ 10 10 GeV or higher. #1 In such a case, the EW vacuum can decay into the deeper vacuum through tunneling in quantum field theory. The lifetime of the EW vacuum has been one of the important topics in particle physics and cosmology.
The decay rate of the EW vacuum has been discussed for a long time. The calculation of the decay rate at the one-loop level first appeared in [15] and was also discussed in other literature [16][17][18][19][20][21][22][23][24][25]. However, there are subtleties in the treatment of zero modes related to the gauge symmetry breaking, which make it difficult to perform a precise and reliable calculation of the decay rate. The lifetime of a vacuum can be evaluated through a rate of bubble nucleation in unit volume and unit time as formulated in [26,27]. The rate is expressed in the form of where B is the action of a so-called bounce solution, and prefactor A is quantum corrections having mass dimension 4. Bounce solution is an O(4) symmetric solution of the Euclidean equations of motion, connecting the two vacua. Although the dominant suppression of the decay rate comes from B, the prefactor A is also important. This is because of large quantum corrections from the top quarks and the gauge bosons. Thus, it is essential to calculate both A and B to determine the decay rate precisely. In the SM, there are infinite bounce solutions owing to (i) the classical scale invariance at a high energy scale, (ii) the global symmetries corresponding to SU(2) L × U(1) Y /U(1) EM , as well as (iii) the translational invariance. For the calculation of the prefactor A, a proper procedure to take account of the effects of the zero modes related to (i) and (ii) were not well understood until recently. In addition, the previous calculations of A were not performed in a gauge-invariant way, which made the gauge invariance of the result unclear.
Recently, a prescription for the treatments of the gauge zero modes was developed [28,29], based on which a complete calculation of the decay rate of the EW vacuum became possible. The calculation has been performed by the present authors in a recent publication [25] and also by [24]. The purpose of this paper is to give more complete and detailed discussion about the calculation of the decay rate. In [25], we have numerically evaluated the functional determinants of fluctuation operators, which are necessary for the calculation of the decay rate. Here, we perform the calculation analytically; a part of the analytic results was first given in [24]. The effects of the zero modes and the modes with the lowest angular momentum are carefully taken into account, on which previous works had some confusion. We give fitting formulae of the functional determinants based on analytic results, which are useful for the numerical calculation of the decay rate. We also provide a C++ package to study the ELectroweak VAcuum Stability, ELVAS, which is available at [30].
In this paper, we discuss the calculation of the decay rate of the EW vacuum in detail. We first present a detailed formulation of the calculation of the decay rate at the one-loop level. We derive a complete set of analytic formulae that can be used for any models that exhibit classical scale invariance at a high energy scale like the SM. Then, as one of the important applications, we calculate the decay rate of the EW vacuum in the SM. We find that the lifetime of the EW vacuum is much longer than the age of the universe. There, we see that one-loop corrections from the top quark and the gauge bosons are very large although there is an accidental cancellation. It shows the importance of A for the evaluation of a decay rate. We also evaluate the decay rates of the EW vacuum for models with extra fermions that couple to the Higgs field. In such models, the EW vacuum tends to be destabilized compared with that of the SM since the quartic coupling of the Higgs field is strongly driven to a negative value. (For discussion about the stability of the EW vacuum in models with extra particles, see [22,.) We consider three models that contain, in addition to the SM particles, (i) vector-like fermions having the same SM charges as the left-handed quark and the right-handed down quark, (ii) vector-like fermions with the same SM charges as left-handed lepton and right-handed electron, and (iii) a right-handed neutrino. We give constraints on their couplings and masses, requiring that the lifetime of the EW vacuum be long enough.
This paper is organized as follows. In Section 2, we summarize the formulation for the decay rate at the one-loop level, where we provide an analytic formula for each field that couples to the Higgs boson. The detail of the calculation is given in Appendices A − D. In Section 3, we evaluate the vacuum decay rate in the SM. Readers who are interested only in the results can skip over the former section to this section. In Section 4, we analyze decay rates in models with extra fermions. Finally, we conclude in Section 5.

Formulation
We first discuss how we calculate the decay rate of the EW vacuum. In the SM, the EW vacuum becomes unstable due to the renormalization group (RG) running of the quartic coupling constant of the Higgs boson, which makes the quartic coupling constant negative at a high scale. In the SM, the instability occurs when the Higgs amplitude becomes much larger than the EW scale. Since the typical field value for the bounce configuration is around that scale, we can neglect the quadratic term in the Higgs potential.
In this section, we use a toy model with U(1) gauge symmetry to derive relevant formulae. The calculation of the decay rate of the EW vacuum is almost parallel to that in the case with U(1) gauge symmetry; the application to the SM case will be explained in the next section.

Setup
Let us first summarize the setup of our analysis. We study the decay rate of a false vacuum whose instability is due to an RG running of the quartic coupling constant of a scalar field, Φ. We assume that Φ is charged under the U(1) gauge symmetry (with charge +1); the kinetic term includes where A µ is the gauge field and g is the gauge coupling constant, while we consider the following scalar potential: The quartic coupling, λ, depends on the renormalization scale, µ, and is assumed to become negative at a high scale due to the RG effect. As we have mentioned before, we neglect the quadratic term assuming that λ becomes negative at a much higher scale. In this setup, the scalar potential has scale invariance at the classical level. In the application to the case of the SM, Φ corresponds to the Higgs doublet and λ corresponds to the Higgs quartic coupling constant.
Hereafter, we perform a detailed study of the effects of the fields coupled to Φ on the decay rate of the false vacuum. We consider a Lagrangian that contains the following interaction terms: where σ is a real scalar field, and ψ L and ψ R are chiral fermions (with relevant U(1) charges). #2 We take y real and κ > 0. We neglect dimensionful parameters that are assumed to be much smaller than the typical scale of the bounce. In addition, gauge fixing is necessary to take into account the effects of gauge boson loops. Following [29,55], we take the gauge fixing function of the following form: Then, the gauge fixing term and the Lagrangian of the ghosts (denoted asc and c) are given by where ξ is the gauge fixing parameter. #3 #2 We assume that there exist other chiral fermions that cancel out the gauge anomaly. #3 In the non-Abelian case, one of ∂ µ in eq. (2.6) is replaced by the covariant derivative. The interaction of the ghosts with the gauge field does not affect the following discussion.
The bounce solution is an O(4) symmetric object [56,57], and transforms under a global U(1) symmetry. Thus, choosing the center of the bounce at r = 0 (with r ≡ √ x µ x µ ), we can write the bounce solution as with our choice of the gauge fixing function, without loss of generality. Here, θ is a real parameter. The function,φ, obeys with boundary conditions ∂ rφ (0) = 0 andφ(∞) = 0. For a negative λ, we have a series of Fubini-Lipatov instanton solutions [58,59]: which is parameterized byφ C (i.e., the field value at the center of the bounce). We also define R, which gives the size of the bounce, as The action of the bounce is given by Notice that the tree level action is independent ofφ C owing to the classical scale invariance.
Once the bounce solution is obtained, we may integrate over the fluctuation around it. We expand Φ as where h and ϕ are the physical Higgs mode and the Nambu-Goldstone (NG) mode, respectively. At the one-loop level, the prefactor can be decomposed as where V 4D is the volume of spacetime, and A (X) is the contribution from particle X. Each of the factors has a form of where M (X) and M (X) are the fluctuation operators around the bounce solution and around the false vacuum, respectively. Here, w (X) = 1 for Dirac fermions and Faddeev-Popov ghosts, and w (X) = −1/2 for the other bosonic fields. The fluctuation operator is defined as second derivatives of the action:

Functional determinant
For the evaluation of the functional determinants, we first decompose fluctuations into partial waves, making use of the O(4) symmetry of the bounce [15]. The basis is constructed from Y J,m A ,m B (Ω), the hyperspherical function on S 3 with Ω being a coordinate on S 3 . The decomposition of each fluctuation is given in Appendix A. Here, J = 0, 1/2, 1, · · · is a non-negative half-integer that labels the total angular momentum in four dimension, and m A and m B are the azimuthal quantum numbers for the A-spin and the B-spin of so(4) ≃ su(2) A × su(2) B , respectively. The four-dimensional Laplacian operator acts on the hyperspherical function as With the above expansion, the functional determinants can be expressed as DetM .

(2.27)
Here, where g is the gauge coupling constant, and Γ(z) is the gamma function.

Zero modes
In the calculation of the decay rate of the EW vacuum with the present setup, there show up zero modes in association with dilatation, translation, and global transformation of the bounce solution. Consequently, M have zero eigenvalues. Their determinants vanish as shown in eqs. (2.24) and (2.27) (see also eq. (A.80)); a naive inclusion of those results gives a divergent behavior of the decay rate, which requires a careful treatment of the zero modes. In the present case, we can consider the effect of each partial wave (labeled by J) separately. Thus, in this subsection, we consider the case where the fluctuation operator M where M (X) The above argument can be applied to the zero modes of our interest. In Appendix B, we obtain the following replacements to take care of the dilatational, translational, and gauge zero modes:

Renormalization
After taking the product over J, we have a UV divergence and thus we need to renormalize the result. In this paper, we use the MS-scheme, which is based on the dimensional regularization. In this subsection, we explain how the divergences can be subtracted using counter terms in the MS-scheme.
Since the dimensional regularization cannot be directly used in this evaluation, we first regularize the result by using the angular momentum expansion as (2.49) We call this regularization as "angular momentum regularization." Here, ε X is a positive number, which will be taken to be zero at the end of the calculation. Since the divergence is at most a power of J, the regularized sum converges. In Appendix C, we calculate the sum analytically, and obtain where the functions, S σ and S ψ , are given in eqs. (C. 13) and (C.16), respectively. In addition, we define Then, the primed quantities are given by where A G ≃ 1.282 is the Glaisher number. Next, we relate the above results with those based on the dimensional regularization in D dimension, which contain the regularization parameter,ε D , defined as with γ E being the Euler's constant. We convert the results based on two different regularizations by calculating the following quantity: Det M (X) div has the same divergence as ln DetM (X) Det M (X) does when δM (X) does not have a derivative operator.
The first line of eq. (2.57) can be calculated by directly evaluating the traces in a momentum space with using the dimensional regularization; the result is denoted as ln On the other hand, the second line of eq. (2.57) can be evaluated as  for p = 1 and 2, and Ψ (0) =Ψ. Then, we calculate We relate the results based on two regularizations by replacing The expressions of [ln A (X) ] div,ε X and [ln A (X) ] div,ε D for each field are given in Appendix D, where we further simplify the relation so that the left-hand side only includes terms that are divergent at the limit of ε X → 0. We summarize the results below: • Scalar field:

65)
• Gauge and NG fields: Subtracting 1/ε D , we obtain the renormalized prefactor in the MS-scheme.

Dilatational zero mode
In the calculation of the decay rate of the EW vacuum, we have an integral over R in association with the classical scale invariance, as we saw in eq. (2.46). So far, we have performed a one-loop calculation of the decay rate, based on which the decay rate is found to behave as where β (1) λ is the one-loop β-function of λ. (Here, we only show the R-and µ-dependencies of the one-loop corrections.) Thus, using the purely one-loop result, the integral does not converge.
We expect, however, that the integration can converge once higher order effects are properly included. To see the detail of the path integral over the dilatational zero mode, let us denote the decay rate as where B eff fully takes account of all the effects of higher order loops. In order to discuss how B eff should behave, it is instructive to rescale the coordinate variable asx as well as the fields asΦ Using the rescaled fields, all the explicit scales disappear from the action as a result of scale invariance: where L is the total Lagrangian,˜ In addition, the rescaled bounce solution is given by (2.78) Based onL and˜ , we expect: • Only positive powers of κ |λ| , y √ |λ| and g √ |λ| appear in the decay rate since there is no singularity when any of these goes to zero. In paticular, they cannot be in a logarithmic function.
• When we renormalize divergences using dimensional regularization, we introduce a renormalization scaleμ. It is always in a logarithmic function and is related to the original renormalization scale asμ = µR.
• Quantum corrections have˜ ℓ−1 at the ℓ-th loop since the loop expansion is equivalent to the˜ expansion.
Based on the above arguments, B eff is expected to be expressed as where P ℓ is the contribution at the ℓ-loop level, and n zero is the number of zero modes. If the effects of higher order loops are fully taken into account, B eff should be independent of µ because the decay rate is a physical quantity; in such a case, we may choose any value of the renormalization scale µ. In the perturbative calculation, the µ-dependence is expected to cancel out order-by-order; as shown in eq. (2.67), we can explicitly see the cancellation of the µ-dependence at the one-loop level [64]. In our calculation so far, however, we only have the one-loop result, in which µ dependence remains. As indicated in eq. (2.81), the µ dependence shows up in the form of ln p µR with p = 1, 2, · · · . If | ln µR| ≫ 1, the logarithmic terms from higher order loops may become comparable to the tree-level bounce action and the perturbative calculation breaks down. In order to make our one-loop result reliable, we should take µ ∼ O(1/R), i.e., we use R-dependent renormalization scale µ. #4 With such a choice of µ (as well as with the use of coupling constants evaluated at the renormalization scale µ), the integration over the size of the bounce is dominated only by the region where |λ(1/R)| becomes largest. In the case of the SM, the integration over the size of the bounce converges with this prescription as we show in the following section.

Final result
Here, we summarize the results obtained in the previous subsections and Appendices. The decay rate with a resummation of important logarithmic terms is given by

86)
#4 This is equivalent to summing over large logarithmic terms appearing in higher loop corrections if we work with a fixed µ. Since we have calculated the decay rate at the one-loop level, it is preferable to use, at least, the two-loop β-function to include the next-to-leading logarithmic corrections.
Here, V G is the volume of the group space generated by the broken generators. The definitions of S σ (z), S ψ (z) can be found in Appendix C. Here, we note that the analytic results for the scalar, Higgs, and fermion contributions were first given in [24] with different expression. We emphasize that the final result does not depend on the gauge parameter, ξ, and hence our result is gauge invariant. The above result is also applicable to the case where the U(1) symmetry is not gauged as explained in Appendix E.
We have also derived fitting formulae of the functions necessary for the calculation of the decay rate; the result is given in Appendix F. The fitting formulae are particularly useful for the numerical calculation of the decay rate. In addition, a C++ package to study the ELectroweak VAcuum Stability, ELVAS, is available at [30], which is also applicable to various models with (approximate) classical scale invariance.

Decay rate
Now, we are in a position to discuss the decay rate of the EW vacuum in the SM. As we have discussed, the decay of the EW vacuum is induced by the bounce configuration whose energy scale is much higher than the EW scale. Thus, we approximate the Higgs potential as #5 where H is the Higgs doublet in the SM and λ is the Higgs quartic coupling constant. Notice that λ depends on the renormalization scale µ; in the SM, λ becomes negative when µ is above ∼ 10 10 GeV with the best-fit values of the SM parameters. In addition, the relevant part of the Yukawa couplings are given by where q L is the left-handed 3rd generation quark doublet, t c R is the right-handed anti-top quark, and y t is the top Yukawa coupling constant. #5 We assume that the Higgs potential given in eq. (3.1) is applicable at a high scale. In particular, we assume that the effect of Planck suppressed operators, which may arise from the effect of quantum gravity, is negligible. For the discussion about the effect of Planck suppressed operators, see [65][66][67][68][69][70][71][72][73].
Assuming that λ < 0, the bounce solution for the SM is given by where σ a is the Pauli matrices and functionφ is given by eq. (2.9). In particular, remember thatφ contains a free parameter, which we choose R, because of the classical scale invariance. The results given in the previous section can be easily applied to the case of the SM. Taking account of the effects of the (physical) Higgs boson, top quark, and weak bosons (as well as NG bosons), #6 the decay rate of the EW vacuum in the SM can be written in the following form: As we have mentioned, the relevant renormalization scale of the integrand is µ ∼ O(1/R); in the following numerical analysis, we take µ = 1/R unless otherwise stated. If λ(µ) is positive, there is no bounce solution; the integrand is taken to be zero in such a case. In addition, since we neglect the mass term in the Higgs potential, 1/R should be much larger than the EW scale. This condition is automatically satisfied in the present analysis because λ becomes negative at the scale much higher than the electroweak scale. The Higgs contribution A ′(h) is given in eq. (2.83), while the top-quark contributions is given by where the factor of 3 is the color factor. As for the gauge contributions, we have SU(2) L × U(1) Y /U(1) EM broken symmetries, instead of U(1) in our previous example. Thus, we have a different volume of the group space, V G . To calculate V G , we first expand H around the bounce solution with θ a = 0 as Here, ϕ 1 and ϕ 2 are NG bosons absorbed by charged W -bosons while ϕ 3 is that by Z-boson.
With the change of θ a , the NG modes are transformed as We checked that the effect of the bottom quark is numerically unimportant.
The integration over the zero modes in association with the gauge transformation of the bounce solution can be replaced by the integration over the parameter θ a : with the above definition of θ a . Then, following the argument given in Appendix B, the gauge contribution is evaluated as where A ′(Aµ,ϕ) is given in eq. (2.87), and with g 2 and g Y being the gauge coupling constants of SU(2) L and U(1) Y , respectively.

Numerical results
Now, let us evaluate the decay rate of the EW vacuum in the SM. The decay rate of the EW vacuum is very sensitive to the coupling constants in the SM. In our numerical analysis, we use [74]: (3.14) where m h and m t are the Higgs mass and the top mass, respectively, while α s is the strong coupling constant. Following [75], the gauge couplings, the top Yukawa coupling, and the Higgs quartic coupling are determined at µ = m t ; the calculations are done in the on-shell scheme at NNLO precision. In addition, we use three-loop QCD and one-loop QED βfunctions [76][77][78] together with values in [74] in order to determine the bottom and the tau Yukawa couplings at µ = m t . First, we show the RG evolution of the SM coupling constants in fig. 1. We use mainly three-loop β-functions summarized in [75] and the central values for the SM parameters. The black dotted line indicates whereφ C reaches the Planck scale M Pl ≃ 2.4 × 10 18 GeV. We also show the running above the Planck scale, assuming there are no significant corrections from gravity. For 10 10 GeV µ 10 30 GeV, λ becomes negative; for such a region, we use a dashed line to indicate λ < 0.
In order to understand the µ dependence of λ, let us show one-loop RG equations of λ and y t (although, in our numerical calculation, we use RG equations including two-and In the shaded region, λ is positive and the integrand is zero. three-loop effects and the contribution from the bottom and the tau Yukawa couplings): At a low energy scale, the term proportional to y 4 t drives λ to a negative value. As the scale increases, y t decreases while g Y increases, which brings λ back to a positive value. Notice that λ is bounded from below in the SM.
We show the integrand of γ in the bottom panel of fig. 1, together with that of They are also shown in a linear scale in the top panel of fig. 2. There are some remarks on the integral over R.
• As indicated by the top panel of fig. 2, the integral is dominated by the interval 10 17 GeV 1/R < 10 18 GeV, corresponding to 10 18 GeV φ C < 10 19 GeV, which is close to the Planck scale. We may formally perform the R-integration up to the scale where λ becomes positive again; the result of such an analysis is denoted as γ ∞ . Otherwise, we may stop the integration atφ C ∼ M Pl , expecting that the SM breaks down at the Planck scale due to an effect of quantum gravity; we also perform such a calculation terminating the integral atφ C = M Pl , assuming that the bounce solution is unaffected by the effect of quantum gravity. The result is denoted as γ Pl .
• As one can see in the bottom panel of fig.1, there is an artificial divergence of the integrand at 1/R ≃ 10 10 GeV. This is due to a breakdown of perturbative expansion owing to a too small |λ|, which makes the one-loop effect larger than the tree-level one. We expect that the effect of such a bounce configuration is unimportant because the bounce action for such a small |λ| suppresses the decay rate significantly. Thus, we exclude such a region from the region of integration. In our numerical calculation, eff is the one-loop contribution to B eff , and [ln A (X) ] MS is a contribution from particle X; the region that does not satisfy these conditions is excluded from the integration. #7 By numerically integrating over R, we obtain #8 where the 1st, 2nd, 3rd, and 4th errors are due to the Higgs mass, the top mass, the strong coupling constant, and the renormalization scale, respectively. (In order to estimate the uncertainty due to the choice of the renormalization scale, we vary the renormalization scale from 1/2R to 2/R.) Currently, the largest error comes from the uncertainty of the top mass. With a better understanding of the top quark mass at future LHC experiment [79][80][81][82][83][84][85], or even with at future e + e − colliders [86], more accurate determination of the decay rate will become possible. One can see that the predicted decay rate per unit volume is extremely small, in particular, compared with H −4 0 ≃ 10 3 Gyr Gpc 3 (with H 0 being the expansion #7 The numerical result is insensitive to the cut parameter, 0.8, as far as only the region where the perturbation breaks down is removed from the integration. In the SM, with the central values of the couplings, the numerical result is not affected much even if we change the number from 0.04 to 1.2. #8 In our previous analysis [25], we used a different renormalization scale, i.e., µ =φ C instead of µ = R −1 . With µ =φ C , the decay rate becomes log 10 γ × Gyr Gpc 3 = −570 for the best-fit values of the SM parameters. Difference between this result and that in [25] is due to the correction of an error in Eq. (29) of [25] (see Eq. (D.38)). The uncertainty related to the choice of µ should be regarded as theoretical uncertainty.
eff . rate of the present universe). Such a small decay rate is harmless for realizing the present universe observed. #9 In fig. 3, we show the decay rate in m h vs. m t plane. In the red region, γ becomes larger than H 4 0 , which we call unstable region. In the yellow region, the EW vacuum is metastable, meaning that 0 < γ < H 4 0 . In the green region, the EW vacuum is absolutely stable because λ is always positive. The dashed, solid, and dotted lines correspond to α s = 0.1192, 0.1181, and 0.1170, respectively. The black dot-dashed contours show log 10 [γ × Gyr Gpc 3 ] = 0, −100, −300, and −1000 with the central value of α s . We also show 68, 95, and 99 % C.L. constraints on the Higgs mass vs. top mass plane assuming that their errors are independently Gaussian distributed. In fig. 3, we terminate the integral atφ C = M Pl , but it does not change the figure as far as the cut-off is not so far from the Planck scale. #10 The value ofφ C at the #9 Cosmologically, the Higgs field may evolve into the unstable region due to the dynamics during and after inflation . We do not consider such cases. It is well known that, currently, our universe is (almost) dominated by the dark energy. If it is a cosmological constant, then our universe will eventually become de Sitter space with the expansion rate of about 56.3 km/sec/Mpc [108]. Then, based on γ Pl , the phase transition rate within the Hubble volume of such a universe is estimated to be 10 −580 Gyr −1 , which may be regarded as a decay rate of the EW vacuum in the SM.
For comparison, we also perform a "tree-level" calculation of the decay rate using eq. (3.17). The results are log 10 γ Thus, the difference between γ and γ (tree) turns out to be rather small. This is a consequence of an accidental cancellation among the contributions of several fields. In the bottom panel of fig. 2, we show individual quantum corrections separately, as well as the total oneloop contribution. We can see that the large quantum correction from the top quark is cancelled by those from the gauge bosons. We have also checked that the unstable region on the m h vs. m t plane shifts upward by ∆m t ≃ 0.2 GeV if we use γ (tree) .

Models with Extra Fermions
So far, we assumed that the SM is valid up to the Planck or some higher scale. However, the decay rate of the EW vacuum may be affected if there exist extra particles. In particular, extra fermions coupled to the Higgs boson may destabilize the EW vacuum because the new Yukawa couplings tend to drive λ to a negative value through RG effects [32-43, 46, 48-52, 54]. Consequently, the decay rate of the EW vacuum becomes larger than that in the SM. Potential candidates of such fermions include vector-like fermions as well as right-handed neutrinos for the seesaw mechanism [109][110][111].
In this section, we consider several models with such extra fermions. We perform the RG analysis of the runnings of coupling constants with the effects of the extra fermions. We include two-loop effects of the extra fermions into the β-functions, which can be calculated using the result in [112][113][114][115]. We also take account of one-loop threshold corrections due to the extra fermions, which are summarized in Appendix G. #11 For the integration over R, we follow the procedure in the SM case, as well as the following treatments: • We terminate the integration if any of the coupling constants (in particular, Yukawa coupling constants of extra particles) exceeds √ 4π.
• In order to maintain the classical scale invariance at a good accuracy, we require 1/R > 10M ex , where M ex is the mass scale of the new particles.

Vector-like fermions
where L (Yukawa) SM is the SM part. We also add the following mass terms: For simplicity, we assume M Q = M D . We take the new Yukawa coupling constants and mass parameters real and positive. We also take As we have mentioned before, the scale dependence of the new Yukawa coupling constants is evaluated by using two-loop RG equations and one-loop threshold corrections (see Appendix G). The calculation of the decay rate is parallel to the SM case, and the decay rate is given in the following form: In fig. 4, we show the contours of the constant decay rate on M D vs. y D plane. Here, we use the central values of the SM parameters. The meanings of the shading colors are the same as in the SM case. The left and the right panels show the results with and without imposing the conditionφ C < M Pl in integrating over R, respectively. As we can see, the effect of such a cut-off is significant. This is becauseφ C at the maximum of the integrand, φ max C , can become much larger than the Planck scale in the case with extra fermions; we showφ max C for the case with vector-like colored fermions in the left panel of fig. 6. (In the upper-left corner of the figure, the value ofφ max C becomes smaller; this is because, in such a region, the Yukawa coupling constants become non-perturbative at a lower RG scale, which gives an upper bound onφ C in the integration over R.). To see the cut-off dependence of the decay rate, we show the constraint with terminating the integration atφ C = 0.1M Pl in the left panel (dashed line). In addition, when M D and y D are small, we have a region of absolute stability. This is because the addition of colored particles makes the strong coupling constant larger than the SM case. It rapidly drives y t to a small value, which makes λ always positive. Requiring that the lifetime of the EW vacuum should be longer than the age of the universe, we obtain y D 0.35 − 0.5 for 10 3 GeV M D 10 15 GeV.
The second example is non-colored extra fermions, having the same SM quantum numbers as leptons. We introduce L (1, 2, − 1 2 ), L (1, 2, 1 2 ), E (1, 1, −1), and E (1, 1, 1), and the Yukawa and mass terms in the Lagrangian are given by respectively. For simplicity, we take M L = M E and adopt the following renormalization conditions The decay rate of the EW vacuum is given by In fig. 5, we show the contours of constant decay rate. Since the extra fermions are not colored, we do not have a region of absolute stability. We observe a larger effect of the cut-off at the Planck scale. This is becauseφ max C is typically large in a wider parameter space, as indicated in the right panel of fig. 6.
Requiring that the lifetime of the EW vacuum should be longer than the age of the universe, we obtain y E 0.4 − 0.7 for 10 3 GeV M E 10 15 GeV. The constraint becomes significantly weaker for larger M E owing to the cut-off at the Planck scale.

Right-handed neutrino
Next, we consider the case with right-handed neutrinos, which is responsible for the active neutrino masses via the seesaw mechanism [109][110][111]. For simplicity, we concentrate on the case where only one mass eigenstate of the right-handed neutrinos, denoted as N, strongly couples to the Higgs boson (as well as to the third generation lepton doublet). Then, the Yukawa and mass terms in the Lagrangian are where, in this subsection, L denotes one of the lepton doublets in the SM. We define Assuming that, for simplicity, the neutrino Yukawa matrix is diagonal in the mass-basis of right-handed neutrinos, the following effective operator shows up by integrating out N: In the gray region, the EW vacuum is absolutely stable.
One of the active neutrino masses is related to the value of C at the EW scale as with v ≃ 246 GeV being the vacuum expectation value of the Higgs boson. In our numerical calculation, we use the following one-loop RG equation to estimate the neutrino mass [116]: In the SM with right-handed neutrinos, the decay rate of the EW vacuum is evaluated with In fig. 7, we show the contour plots of the decay rate. Since it does not have any SM charges, the decay rate goes to the SM value when y N goes to zero. The effect of the cut-off at the Planck scale is again large, which is because of a largeφ max C as shown in fig. 8. The purple solid lines show the left-handed neutrino mass. Requiring that the decay rate should be smaller than the age of the universe, we obtain y N 0.65 − 0.8 for 10 12 GeV M N 10 15 GeV.

Conclusion
In this paper, we have calculated the decay rate of the EW vacuum in the framework of the SM and also in various models with extra fermions. We included the complete one-loop corrections as well as large logarithmic terms in higher loop corrections. We used a recently developed technique to calculate functional determinants in the gauge sector, which not only gives a prescription to perform a gauge invariant calculation of the decay rate but also allows us to calculate the functional determinants analytically. In addition, in calculating the decay rate of the EW vacuum, zero modes show up in association with the dilatational and gauge symmetries. We have properly taken into account their effects, which was not possible in previous calculations. We have given an analytic formula of the decay rate of the EW vacuum, which is also applicable to models that exhibit classical scale invariance at a high energy scale. The decay rate of the EW vacuum is sensitive to the coupling constants in the SM and their RG behavior. We have used three-loop RG equations for the study of the RG behavior of the SM couplings. The result is used for the precise calculation of the decay rate of the EW vacuum. The decay rate of the EW vacuum is estimated to be log 10 γ Pl × Gyr Gpc 3 = −582 +40 +184 +144 +2 −45 −329 −218 −1 , where the errors come from the Higgs mass, the top mass, the strong coupling constant, and the renormalization scale, respectively. Here, only the bounce configurations with its amplitude smaller than the Planck scale is taken into account; for the decay rate of the EW vacuum, such a cut-off of the bounce amplitude does not significantly affect the result as far as it is not so far from the Planck scale. Since H −4 0 ≃ 10 3 Gyr Gpc 3 , the lifetime is long enough compared with the age of the universe.
We have also considered models with extra fermions. Since they typically make the EW vacuum more unstable, the constraints on their masses and couplings are phenomenologically important. We have analyzed the decay rate for the extensions of the SM with (i) vectorlike quarks, (ii) vector-like leptons, and (iii) a right-handed neutrino. We have obtained a constraint on the parameter space for each model, requiring that the lifetime be long enough. The results constrain the Yukawa couplings that are larger than about 0.3 − 0.5 if we do not consider the cut-off at the Planck scale. The effect of the cut-off was found to be rather large and the constraints on the Yukawa couplings become weaker, at most, by 0.3 after including the cut-off. been corrected in the recent revision. The differences between our numerical results and those in [24] come mainly from the difference of the threshold corrections to the MS top Yukawa coupling constant, which is regarded as a theoretical uncertainty. In addition, there is a difference in the treatment of the integration over the bounce size, although it has little effect on the numerical results.

A Functional Determinant
In this Appendix, we present analytic formulae for various functional determinants. For simplicity, we consider the case with U(1) gauge interaction. The charge of the scalar field that is responsible for the instability, Φ, is set to be +1. Application of our results to the case of general gauge groups is straightforward. We are interested in the case where the Lagrangian has (classical) scale invariance; the potential of Φ is given in eq. (2.2), and the bounce solution is obtained as eq. (2.7).

A.1 Scalar contribution
We first consider a real scalar field, σ, which couples to Φ as where κ is a positive coupling constant. The contribution to the prefactor is given by First, we expand σ into partial waves; Here and hereafter, α J,m A ,m B denotes the radial mode function of X. For notational simplicity, the summations over J, m A , and m B are implicit.
Since the fluctuation operator for partial waves does not depend on m A and m B , we have (2J + 1) 2 degeneracy for each J. Summing up all the contributions, Then, using eq. (2.21), we have where the function f (σ) satisfies and lim r→0 f (σ) (r)/r 2J = 1. The solution of the above differential equation is given by where 2 F 1 (a, b; c; z) is the hypergeometric function and Taking the limit of r → ∞, we have . (A.10)

A.2 Higgs contribution
Using the bounce solution given in eq. (2.9), the Higgs mode fluctuation is parametrized as Then, the Higgs contribution to the prefactor is given by The functional determinant of the Higgs mode can be obtained with the same procedure as the case of the scalar contribution, taking κ → −3|λ|: . (A.13) As we can see, the above ratio vanishes for J = 0 and J = 1/2, which are due to the scale invariance and the translational invariance, respectively.

A.3 Fermion contribution
Let us consider chiral fermions ψ L and ψ R with the following Yukawa term 14) The contribution to the prefactor is given by Taking the basis given in [117], we expand it as Solutions of eq. (A.19) can be expressed by using two functions χ (ψ) and η (ψ) as where i = 1 and 2 are for two independent solutions of eq. (A. 19), and χ (ψ) and η (ψ) obey For the first solution, we take where the function f (ψ) satisfies and The analytic formula of f (ψ) is given by and f (ψ) behaves as For the second solution, we take Using these solutions, we get and hence (A.33)

A.4 Gauge contribution
We consider the contributions from gauge bosons, NG bosons, and Faddeev-Popov ghosts. The Lagrangian is given in the following form: where F µν is the field strength tensor, and where For the partial wave expansions, we use the following basis; where V and M (T ) can be constructed from the functions, χ i , η i , and ζ i , as [28,29] where χ i , η i , and ζ i obey We also note useful relations: The first solution is obtained by setting ζ 1 = 0 and η 1 = 0 as For the second solution, we set ζ 2 = 0 and where the function f (η) satisfies and We can find an analytic formula of f (η) as Then, we get .
The last solution can be obtained with The asymptotic form of η 3 is given by We also need three independent solutions around the false vacuum. We take Then, using eq. (2.21) with combining the above expressions, we obtain the ratio of the functional determinants for S, L, and NG modes as .
The functional determinant for T modes can be obtained by using the result for the scalar contribution. With the replacement κ → g 2 in eq. (A.10): .
For J = 0, the solutions of M (S,ϕ) 0 Ψ = 0 can be decomposed as [28,29] with i = 1 and 2 for two independent solutions, where the functions χ i and ζ i satisfy Notice that there are useful relations: 1 The first solution is given by The second solution can be obtained with giving rise to The vanishing of the above ratio is due to the existence of zero mode. The treatment of the zero mode is discussed in Appendix B.

B Zero Modes
In this Appendix, we discuss the zero modes associated with the dilatation, translation, and global transformation.

B.1 Dilatational zero mode
In J = 0 mode of the Higgs fluctuation, there exists a dilatational zero mode. The dilatational transformation is parametrized byφ C ; withφ C →φ C + δφ C , such a change can be absorbed by where we neglect higher order terms in δφ C , and The second term in the right-hand side of eq. (B.1) is nothing but the change of the amplitude of the dilatational zero mode; one can easily check M 0G D = 0. In order to translate the path integral over the dilatational zero mode to the integration overφ C , we calculate Notice that the condition (2.40) is satisfied with the above choice of ρ D , i.e., The ratio of the determinant can be obtained from eq. (A.10) by replacing κ → −3|λ| + 15 16π |λ| 2φ2 C ν. Expanding the result with respect to ν, we get Thus, we obtain Here, we take an absolute value since there is a negative mode [26,27].

B.2 Translational zero mode
In J = 1/2 mode of the Higgs fluctuation, translational zero modes exist. The translation is parametrized by the center of the bounce. The shift of the center of the bounce to a µ can be absorbed by the transformation of the Higgs mode as h → h + a µGT Y 1/2,µ + · · · , (B.8) #13 A prescription used in [25] is consistent with our argument, where ρ D is a constant. whereG and withx µ = x µ /|x|. Notice that Y 1/2,µ is given by a linear combination of Y 1/2,m A ,m B with the same normalization as Y 1/2,m A ,m B . Thus, as noted in [27], the path integral over the translational zero modes can be converted to the integration over the spacetime volume.
In order to take care of the translational zero mode, we calculate One can see that which is consistent with eq. (2.40). Following the argument in Subsection 2.3, we have The functionf (B.17) Thus, we obtain where V 4D is the spacetime volume.

B.3 Gauge zero mode
In J = 0 mode of the gauge field, we have a gauge zero mode. For the case of the U(1) gauge symmetry, the bounce solution is parameterized as eq. (2.7) with the parameter θ. The path integral over the gauge zero mode can be understood as the integration over the variable θ, as we see below.
The effect of the shift θ → θ + δθ can be absorbed by the transformation of the NG mode as Using the equation of motion of the bounce solution, the following relation holds: In order to deal with gauge zero mode, we calculate We can decomposeΨ 1 asΨ

C Infinite Sum over Angular Momentum
In this Appendix, we perform various infinite sums appearing in the calculation of functional determinants. We first evaluate the following sum: which can be used for the calculation of the scalar contribution to the prefactor. Notice that J = 0, 1/2, 1, · · · is half-integer. In addition, here, z is a complex number satisfying For ε σ > 0, the sum converges thanks to the factor 1/(1 + ε σ ) 2J .
First, we rewrite the log-gamma functions with integrals of digamma functions as where ψ Γ (z) is the digamma function. Then, we use the following relation: which is valid for ℜ(x) > 0 and ℜ(y) > 0, and we obtain Notice that this is verified only in region (C.2). Then, we interchange the two integrals, #14 and integrate over x first. Consequently, it becomes Notice that the integration over u is convergent.
Since we have regularized the sum, the result should be finite. Thus, we can take the sum first: We can see that new poles appear at u = 1 + ε σ , but the integral is still convergent for positive ε σ . Then, using we obtain #14 This is justified when Interchanging the integrals and integrating over u, we get where 2 F 1 (a, b; c; z) is the hypergeometric function. Since we do not need higher order terms in ε σ , we expand I σ as where H(z) is the harmonic number. After the final integral, we obtain We can repeat a similar calculation for (C.14) which can be used for the calculation of the fermionic contribution to the prefactor. Here, −2 < ℜ(z) < 2. The result is Next, we consider which is for the calculation of the Higgs contribution. Using a similar technique as in the previous cases, we obtain Then, performing the sum first, I h is obtained as Finally, for the gauge and NG contributions, let us consider (C.20) Similary to I h , it is expressed as which results in

D Renormalization with MS-scheme
In this Appendix, we relate the regularization based on the angular momentum expansion, which we call "angular momentum regularization," and the dimensional regularization. In particular, we derive the relation between ε X , which shows up in the angular momentum regularization, andε D , which is for dimensional regularization.

D.1 Scalar field
Let us start with the scalar contribution. For this purpose, let us consider ln A (σ) div , defined in (2.57).
First, we calculate ln A (σ) div with angular momentum regularization with using ε σ as a regularization parameter. The expansion of eq. (2.59) is exactly the same as that with respect to κ. Thus, we get ln where ψ (n) (z) is the polygamma function. Summing over J, we get As is expected, it has the same divergence as ln A (σ) εσ does. Next, we directly calculate ln A (σ) div by using the dimensional regularization. Using the ordinary Feynman rules, we obtain where F [· · · ] is the Fourier transform of the argument. With performing the integration, we obtain Comparing eqs. (D.2) and (D.4), we obtain the relation between the two regularizations as

D.2 Higgs
For the Higgs case, the relation between ε h andε D can be obtained from eq. (D.5) with the replacement of κ → −3|λ|: Notice that the zero modes have nothing to do with the UV divergence.

D.3 Fermion
Next, we derive the relation for the fermionic contribution. As discussed in [24], it is convenient to expand with respect to y. The expansion in eq. (2.57) is equivalent to Then, together with eq. (A.33), we have y 2 |λ| + 1 18 We can also calculate ln A (ψ) div with dimensional regularization: (D.12) Thus, we obtain (D.13)

D.4 Gauge and NG fields
Finally, we consider the gauge and NG contributions. Although we may use eq. (2.57) to subtract the divergence, it is more convenient to use the expression of the prefactor with a different choice of the gauge fixing function from that in eq. (A.37).
Here, we use the result with the following choice of the gauge fixing function: which we call the background gauge. We may perform a calculation of the prefactor A with this choice of the gauge fixing function. We have checked that, irrespective of the choice of ξ, the J > 0 contribution from the gauge and NG sectors agrees with the one that we have obtained using the gauge fixing function in (A.37). We also comment here that, in the background gauge, the treatment of the gauge zero mode is highly non-trivial. However, such an issue is unimportant for the following discussion because the zero mode does not affect the behavior of the divergence, which we will discuss below. For simplicity, we take ξ = 1 in the following. Then, using the same basis as eq. (A.42), we define ln A where the fluctuation operator in the background gauge is given by With the angular momentum decomposition, we have ln A  Following [28], we can show that ln DetM (2) J (C) can be also expressed as .

E Vacuum decay with global symmetry
For completeness, we discuss the case where the field that is responsible for the decay transforms under a global symmetry, although it is not the case of the SM Higgs field. In such a case, we need to take into account quantum corrections from the associated NG bosons. Similarly to the gauge contributions, the NG fluctuation operator has zero modes in association with the breaking of the global symmetry. Let us consider a U(1) symmetry, for simplicity. The contribution from the NG boson, ϕ, is given by As we can see, there is a zero mode for J = 0. Since it is obtained in the limit of g → 0 in eq. (2.27), ln A (ϕ) MS is given by eqs. (2.86) and (2.87) with taking g = z g = 0.

F Numerical Recipe
In this Appendix, we give fitting formulae of the prefactors at the one-loop level. Contrary to the analytic formulae including various special functions with complex arguments, which may be inconvenient for numerical calculations, the fitting formulae give a simple procedure to perform a numerical calculation of the decay rate with saving computational time. Compared to the analytic expressions, the errors of the fitting formulae are 0.05% or smaller. A C++ package based on our fitting formulae for the study of the ELectroweak VAcuum Stability, ELVAS, can be found at [30].

G Threshold Corrections
In this Appendix, we summarize the one-loop threshold corrections for the coupling constants in the models with extra fermions discussed in Section 4. We parameterize the threshold corrections as where c(below) is a coupling constant below the matching scale, while c is that above the scale. For the notational simplicity, we only show the quantity ∆ c for each coupling constant.
Notice that ∆ c depends on the extra fermion mass M X (X = D, E, N). In our analysis, we take the matching scale to be equal to M X .