Scale Invariant Instantons and the Complete Lifetime of the Standard Model

In a classically scale-invariant quantum field theory, tunneling rates are infrared divergent due to the existence of instantons of any size. While one expects such divergences to be resolved by quantum effects, it has been unclear how higher-loop corrections can resolve a problem appearing already at one loop. With a careful power counting, we uncover a series of loop contributions that dominate over the one-loop result and sum all the necessary terms. We also clarify previously incomplete treatments of related issues pertaining to global symmetries, gauge fixing and finite mass effects. In addition, we produce exact closed-form solutions for the functional determinants over scalars, fermions and vector bosons around the scale-invariant bounce, demonstrating manifest gauge invariance in the vector case. With these problems solved, we produce the first complete calculation of the lifetime of our universe: 10^139 years. With 95% confidence, we expect our universe to last more than 10^58 years. The uncertainty is part experimental uncertainty on the top quark mass and on ${\alpha}s$ and part theory uncertainty from electroweak threshold corrections. Using our complete result, we provide phase diagrams in the $mt/mh$ and the $mt/{\alpha}s$ planes, with uncertainty bands. To rule out absolute stability to $3{\sigma}$ confidence, the uncertainty on the top quark pole mass would have to be pushed below 250 MeV or the uncertainty on ${\alpha}s(mZ)$ pushed below 0.00025.


Introduction
Tunneling through a barrier is a quintessentially quantum phenomenon. In quantum mechanics (QM), tunneling has been studied analytically, numerically and experimentally, leading to a consistent and comprehensive picture of when and how fast tunneling occurs. In quantum field theory (QFT), much less is known. In QFT, one cannot solve the Schrödinger equation, even numerically, due to the infinite dimensionality of the Hilbert space. The only approach to calculating tunneling rates in QFT seems to be through the saddle point approximation of the path integral [1][2][3][4][5]. This approach involves analytic continuation in an essential way. Because tunneling in QFT has important implications, such as for the stability of the Standard Model vacuum [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] and because QFT tunneling rates are nearly impossible to measure experimentally, it is critical to make sure the rather abstract formalism is actually capable of calculating something physical.
A number of the subtleties in going from QM to QFT were resolved long ago, some more recently, and some challenges still exist. For example, while tunneling rates are physical and therefore should be gauge invariant, it has been challenging to check directly that this is the case. Although exact non-perturbative proofs of gauge-invariance exist [21,22] and there have been many investigations into gauge-dependence [18,[23][24][25][26][27][28][29][30][31][32][33][34][35][36] it has not been shown that gaugeinvariance holds order-by-order in perturbation theory, as it does for S-matrix elements. For some context, recall that for the simpler question of whether a state is absolutely stable in the quantum theory, it was found that the corresponding bound was gauge-dependent with then-current perturbative methods [17,37,38]. The problem was traced to an inconsistent power counting and improper use of the renormalization group equations. A consistent method was recently developed in [17,38], with non-negligible implications for precision top and Higgs-boson mass bounds in the Standard Model. Recently, progress was made in understanding the gauge invariance of tunneling rates by Endo et al. [39,40]; these authors showed explicitly that the rate is gauge-invariant to one loop for general massive scalar scalar field theory backgrounds and we build upon their results.
In fact, gauge invariance is the least of our worries. In order to produce a precision calculation of the tunneling rate -or even the leading order rate with the correct unitsone must understand a whole slew of subtleties not relevant for the absolute stability bound. First of all, there are suspicious elements in the common derivations [41][42][43][44][45][46] of the Callan-Coleman decay rate formula [3,4]. The leading-order confusion is that the rate is said to be determined, even in QM, by taking the imaginary part of a|e −HT |a , a manifestly real expression. The resolution of this paradox involves not analytic continuation of the potential, as is often cited, but the specification of complex paths to be integrated over in the path integral [47,48]. A more physical derivation of a decay rate in QFT was presented recently in [49,50]. Some elements are reviewed in Section 2.
Even if we ignore gauge-dependence and trust the decay rate formulas, we encounter a new roadblock in trying to evaluate tunneling rates in QFTs like the Standard Model, due to classical scale invariance. The basic issue with scale invariance can be seen in the Gaussian 3 approximation to the path integral around a reference field configuration φ b : One typically evaluates the right-hand-side by expanding the fluctuations φ in a basis of eigenfunctions of the operator S [φ b ]. If the action has a symmetry spontaneously broken by φ b , then there will be fluctuation directions φ 0 with zero eigenvalue, that is, for which S [φ b ]φ 0 = 0. Integrating over dξ 0 in the field direction φ = ξ 0 φ 0 then leads to an infrared divergence in Eq. (1.1). Examples include the zero modes associated with translation invariance where φ 0 ∝ ∂ µ φ b or scale invariance where φ 0 ∝ (1 + x µ ∂ µ )φ b . For translations, the infrared divergence is expected -it generates a factor of V T so that the rate is extrinsic, a decay rate per unit volume. For scale invariance, the infrared divergence has no natural volume cutoff and so the decay rate is apparently infinite. Anyone with even minimal familiarity with QFT would immediately guess that the resolution to the scale-invariance divergence is related to dimensional transmutation [51], that the classical scale invariance is broken by quantum effects. Unfortunately, connecting the β-functions to the decay rate calculation within a consistent perturbative framework has remained elusive. In fact, there are two related technical difficulties.
First of all, to integrate over a zero-mode fluctuation, one must use a collective coordinate [1,[52][53][54][55] rather than an infinitesimal fluctuation. For example, with translations, one must integrate over x µ 0 parametrizing fields φ b (x µ + x µ 0 ) in the path integral before the Gaussian approximation is applied (in the middle expression in Eq. (1.1) not the rightmost one). The difference between x µ 0 and coefficients ξ µ of fluctuations in the ∂ µ φ b direction is a Jacobian factor J = d 4 x(∂ µ φ b ) 2 . For translations, this Jacobian is finite. For scale transformations, one wants to move from linearized fluctuations φ = ξ d φ d proportional to the dilatation mode φ d = (1 + x µ ∂ µ )φ b to a collective scale coordinate R. Unfortunately, in this case, the Jacobian factor J = d 4 xφ 2 d is infinite. Related Jacobians for the spontaneously broken SU (2) × U (1) symmetry of the Standard Model are also infinite [11,19,56].
The second problem is that even if one could regularize the Jacobian and go to collective coordinates, the resulting integral dR over scales R would still be infinite. While quantum corrections do break scale invariance at some order, they do not resolve the infinity in the one loop approximation. Indeed, the R dependence of the integrand can be deduced from renormalization group invariance. As we review in Section 3 at one loop, the integral is still infinite. While there is R dependence at higher-loop order, for the higher-loop effects to cancel the infinity form the one loop integrand would require a diversion from the usual loop power counting. This is certainly possible, as the unusual power counting of the Coleman-Weinberg model [51] is often required to extract physical predictions from the effective potential [17,23,38], but a proper power counting for the deca rate does not seem to have been explored in the literature.
A number of unsatisfying approaches to resolve the two problems with the dilation mode have been used in the literature. One method is to impose a scale on the bounce by hand, by demanding a constraint be satisfied [57][58][59][60], such as φ 3 = Λ 3 for some fixed Λ. Then one can try to split the path integral into integrations around the constrained instanton and integrations over Λ. This approach seems impractical, as explicit constrained instantons are hard to find [59] and the Jacobian to go between R and Λ is no simpler than between ξ d and R. While constrained instantons are helpful in understanding how a scalar mass can be a small perturbation, as we discuss in Section 7, they are irrelevant to resolving the integral over R.
In practice, for the decay rate in the Standard Model, people always just invoke dimensional analysis [11,19,61,62]: cut off the divergence in the Jacobian by the Higgs mass and assume the integral over R is dominated by the bubble size R m with the maximal rate. This seems to us a bit cavalier. After all, the fate of the universe is on the line.
In this paper, we provide definitive resolutions to both challenges associated with the dilatation mode. To regularize the Jacobian, a powerful approach has been known for some time but has not been widely appreciated [58,[63][64][65]. The approach is a based on a powerful stereographic projection into five dimensions, where symmetries are manifest and the natural inner product is the conformal metric on the 4-sphere. By undoing this transformation, it becomes clear that the projection is not actually necessary. The key is simply that a zero mode satisfying S [φ b ]φ d = 0 is also a zero mode of any rescaled operator f (x)S [φ b ]φ d . Since the rescaling factor f (x) changes the norm on the eigenfunctions, one can choose f (x) so that φ d is normalizable and the Jacobian is finite. More explicitly, for which allows the spectrum of eigenfunctions to be found in closed form, and the same basis to be used for fluctuations around φ b and around φ = 0. We explain this procedure in Section 3.1.
Once collective coordinates have been invoked, we address the issue of the proper power counting required to evaluate the integral over R. Indeed, at one loop, the integral over R is infinite. At two loops, there is an exp(− ln 2 R) term which makes the integral over R finite. Despite the suppression of this term, the integral scales as −1/2 compared to the one loop integral, making it divergent as → 0. In other words, the two loop result is parametrically more important than the one loop result, a scaling essential to regulating the divergence. At three loops, the integral scales like −3/2 , so that three loop is more important than two loops. Conveniently, for four loops and higher, the integral stills scales like −3/2 compared to one loop. We show that the entire series of leading contributions can be summed in closed form.
With these scale invariance problems solved, we proceed to compute the functional determinants around the bounce over real scalar, complex scalar, vector boson and fermionic fluctuations. We produce for the first time exact formulas for the path integrals in each case. In the gauge boson case, we work in a general 1-parameter family of Fermi gauges and show the result is gauge invariant and that it agrees with the result in R ξ gauges as well.
Applying our exact formulas to the Standard Model, we update the famous stability/metastability phase diagram. For the first time, we can give an exact NLO prediction for the instability/metastability phase boundary. We find that with current data, the dominant uncertainties are from the top quark mass and α s , and these are both comparable to the theory uncertainty from electroweak threshold corrections, currently known to NNLO.
Although this paper is rather long, we have tried to compartmentalize it into more-or-less self-contained sections. Section 2 reviews how tunneling rates are computed. This section is very brief and the interested reader is referred to [50] for more details. Section 3 contains new results about resolving the problems associated with scale invariance of the classical action. Section 4 introduces the methods we will use in later sections to generate exact expressions for functional determinants. The longest section is Section 5 which computes one-by-one the functional determinants for scalars, vectors and fermions. For a reader just interested in the final formulas, these are summarized in Section 5.5. The application to the Standard Model is in Section 6. We tie up one lose end about the finite Higgs boson mass in Section 7. Our results are digested, including a summary of the SM bounds and limits in our conclusions, in Section 8.

Tunneling Formulas and Functional Determinants
In this section we review how to compute a decay rate in quantum field theory. We give some formal expressions for the rate in Section 2.1 and show how to use the saddle point approximation and how to evaluate functional determinants around non-trivial backgrounds in Section 2.2.

Defining the Decay Rate
Suppose our QFT has a metastable extremum localized around the constant classical field configuration φ(x) = φ FV . We would like to compute the lifetime of φ FV or equivalently the rate to tunnel form φ FV through the energy barrier to any other field configuration. The basic tunneling formula was introduced into high-energy theory by Coleman and Callan in 1977 [4], although it has roots in earlier condensed matter treatments (e.g. [1]). The formula is explained at length in Coleman's famous Erice lectures [66], as well as in numerous textbooks [44][45][46]67]. Most of these treatments start with the premise that the decay rate can be computed by evaluating Unfortunately, this formula cannot be correct as written, since the matrix element, and path integral over real paths, are purely real. We would like the T → ∞ limit to pick out the energy of an state localized near the false vacuum whose imaginary part is to give the decay rate. This is certainly the intuition behind Eq. (2.1). Instead, Eq. (2.1) picks out the true ground state energy E 0 which is real. To see this, we note that there are three relevant paths through field space satisfying the boundary conditions φ(− T 2 ) = φ( T 2 ) = φ FV : 1) The constant state φ = φ FV ; 2) the bounce, interpolating from φ FV at Euclidean time τ = ±T /2 to a bubble of shape φ b ( x) in a small time window (hence the name instanton) near τ = 0; and 3) the shot, φ shot matching φ FV at τ = ±T /2 but hovering near the true vacuum φ 0 over most of Euclidean time [49,50]. While Re[S[z]] Figure 1: (Left) The real part of the action along family of field configurations φ(z) parameterized by a complex parameter z. z is chosen so that real z passes through φ FV (green dot), φ b (red dot), and φ shot (blue dot). (Right) shows a top-down view. Integrating along real field configurations only (black dashed contour) makes the path integral real. The decay rate must be calculated by integrating along steepest descent contours (red and green contours) which involve complex field configurations.
the hope is for the path integral to pick out the bounce configuration at large T , instead it picks out the shot since the shot has smaller action, with the result that the path integral is real. In Fig. 1 we sketch of the real part of the action along paths z passing through φ FV , φ b and φ shot . In order for the path integral and energies to be complex we must introduce a unitaryviolating unphysical deformation of the theory. This deformation should prevent flux from returning to the false vacuum so that the strict T → ∞ limit can be taken. For example, we could impose Gamow's outgoing-wave-only boundary conditions to solve the Schrödinger equation [68]. A more formal deformation commonly used is analytic continuation of the classical potential. For example, with a potential V (φ) = m 2 φ 2 + λφ 4 if λ is positive then energies are real, but if λ is negative energies are complex. The analog of this potential in non-relativistic quantum mechanics has been studied for decades [69][70][71][72][73][74]. Unfortunately, the analytic continuation method only works for (unphysical) situations in which the potential is unbounded from below. When the potential is unbounded, no flux can return, but also the path integral over real field configurations is divergent so one has no choice but to change the integration domain to over complex paths. If the path integral is finite, as in the Standard Model, then it is analytic in a domain around real paths and analytic continuation simply reproduces the original (real) result [49,50].
A useful way to define the rate for physical situations, where the action is bounded from below, is to change the integration domain from real field configurations to field configurations associated with steepest descent contours (i.e. those for which ImS[φ] = 0) [47,48]. To 7 Figure 1: (Left) The real part of the action along family of field configurations φ(z) parameterized by a complex parameter z. z is chosen so that real z passes through φ FV (green dot), φ b (red dot), and φ shot (blue dot). (Right) shows a top-down view. Integrating along real field configurations only (black dashed contour) makes the path integral real. The decay rate must be calculated by integrating along steepest descent contours (red and green contours) which involve complex field configurations. the hope is for the path integral to pick out the bounce configuration at large T , instead it picks out the shot since the shot has smaller action, with the result that the path integral is real. In Fig. 1 we sketch of the real part of the action along paths z passing through φ FV , φ b and φ shot .
In order for the path integral and energies to be complex we must introduce a unitaryviolating unphysical deformation of the theory. This deformation should prevent flux from returning to the false vacuum so that the strict T → ∞ limit can be taken. For example, we could impose Gamow's outgoing-wave-only boundary conditions to solve the Schrödinger equation [68]. A more formal deformation commonly used is analytic continuation of the classical potential. For example, with a potential V (φ) = m 2 φ 2 + λφ 4 if λ is positive then energies are real, but if λ is negative energies are complex. The analog of this potential in non-relativistic quantum mechanics has been studied for decades [69][70][71][72][73][74]. Unfortunately, the analytic continuation method only works for (unphysical) situations in which the potential is unbounded from below. When the potential is unbounded, no flux can return, but also the path integral over real field configurations is divergent so one has no choice but to change the integration domain to over complex paths. If the path integral is finite, as in the Standard Model, then it is analytic in a domain around real paths and analytic continuation simply reproduces the original (real) result [49,50].
A useful way to define the rate for physical situations, where the action is bounded from below, is to change the integration domain from real field configurations to field configurations associated with steepest descent contours (i.e. those for which ImS[φ] = 0) [47,48]. To be precise, the rate per unit volume for the formation of bubbles with the shape φ b ( x, 0) is given by 1 2 Here S is the Euclidean action, V the volume of space, and T is a time for which the transition rate has exponential behavior (see discussion in [49,50]). The contour C FV (green contour in Fig. 1) is the steepest descent trajectory through field space passing through φ FV . The contour C b (red contour in Fig. 1) is the steepest descent trajectory passing through the bounce. Note that if we ignore these contour prescriptions and just integrate over real field configurations, along the black dashed line in Fig. 1, passing through φ FV , φ b and φ shot , there is no imaginary part and the rate defined this way is zero, similar to Eq. (2.1). For additional perspective, and insight into the factor of 1 2 , we can alternatively write the decay rate as 1 2 The contour C FV passes through real field configurations until the saddle point φ = φ b is reached when it veers into complex field space (even for real φ) traveling along C b as in Fig. 1. In contrast, if we reverse the trajectory C b , as it passes through φ b it does not head towards φ FV , since Re(−S[φ]) increases in that direction, rather it continues into the conjugate complex field space. Thus integrating along C b gives twice the imaginary part of the integral along C FV . The doubling of the contour is the origin of the factor of 1 2 in Eq. (2. 2) The explanation of why these arcane contour prescriptions produce the decay rate is given in [49,50]. Briefly, the idea is to start by relating the tunneling rate to the time derivative of the probability R d 3 x|ψ(x)| 2 for a state to be found in a region R on the other side of the energy barrier. This leads to the formula

Functional Determinants and Zero Modes
Since the decay rate is defined by path integrals along steepest descent contours, we can compute these path integrals in the saddle-point approximation. To quadratic order, Eq. (2.2) reduces to 1 2 To evaluate the these path integrals, we must be precise about the integral measure. We do this by expanding the fields in some basis φ j : The path integral measure can then be defined as Dφ = j dξ j . An orthogonal basis is naturally provided by eigenfunctions of an operator. It is often convenient to take the operator to be S [φ b ], so that To find the inner product on these basis functions, we note that where integration-by-parts has been used in the last step. So functions with different eigenvalues are orthogonal according to the inner product φ i |φ j = d 4 xφ i φ j . It is also convenient to normalize the fluctuations so that Then we find The point of the normalization convention in Eq. (2.9) is to make removing a normalized fluctuation equivalent to removing its eigenvalue from the product in Eq. (2.10). If one of the eigenvalues is negative, then this expression (after analytic continuation) will have an an imaginary part, as desired. There is at most a single negative eigenvalue [75]. It corresponds to the bounce being a local maximum of the action on the direction from φ FV (see Fig. 1) but a local minimum in all other directions.
If there are zero eigenvalues, then Eq. (2.10) is infinite. Examples are the translation modes, which are proportional to ∂ µ φ b . To check, using S [φ b ] = 0 and that S[φ] has no explicit position dependence, we find confirming that ∂ µ φ are zero modes. To set the normalization of these modes according to our convention, we note that Thus the rescaled modes Separating out the translation modes Eq. (2.6) becomes To integrate over translations, we use collective coordinates [1,[52][53][54][55], parametrizing fields with By expanding Eq. (2.14) for small x µ 0 and comparing to Eq. (2.13) we see that the Jacobian to go from ξ µ to x µ 0 is Then the path integral can be written as where det refers to the functional determinant with the zero eigenvalues taken out by hand and N some (infinite) constant. Noting that the integral over d 4 x 0 gives the volume of euclidean space time, we find

Scale Invariance
Any classically scale invariant action will admit an infinite family of bounces related by scale transformations. To be explicit, we take the potential V (φ) = 1 4 λφ 4 and assume throughout this paper that λ < 0. Then there is a 5-parameter family of bounces given by These Fubini-Lipatov instantons [76,77] all satisfy φ b − λφ 3 b = 0 and all have the same Euclidean action There are four normalizable fluctuations around the bounce corresponding to translations These are handled using collective coordinates as discussed above. We therefore take x 0 = 0 without loss of generality. We also use r = √ x µ x µ as our radial coordinate so that the bounce is The (unnormalized) dilatation mode is Like the translation modes, it is an eigenfunction of the second variation of the action around the bounce with zero eigenvalue: We would like to proceed, as with translations, by going from linearized fluctuations φ = ξ d (N d φ d ), with N d a normalization factor, to the collective coordinate R so that we can write with det now having the φ d dilatation mode removed and J d the Jacobian.
There are two problems with this. First, the Jacobian is infinite: Second, the integral over R is divergent. Even including the one-loop R-dependence from dimensional transmutation, as required at this order, is not enough to remove this infrared divergence.
In the literature, for the first problem J d is often assumed to be made finite by natural infrared cutoff of the scalar mass [11,19]. Unfortunately, the scalar mass adds more problems than it solves -adding a m 2 φ 2 term to the Lagrangian removes all bounces from the solution to the equations of motion. Moreover, adding an infrared cutoff seems to miss the point. Why is the Jacobian infinite in the first place? Going from small linear fluctuations around a bounce to fluctuations corresponding to exact scale transformations seems perfectly reasonable and therefore should be non-singular. In fact, the mass term is irrelevant to the problem (see Section 7).
For the second problem, of the IR divergent integral over R, it is common to pick the scale R for which the leading order result Γ = exp( 8π 2 3λ(R −1 ) ) is maximal, with λ(µ) the running coupling, and evaluate the R integral by dimensional analysis [11,19,61,62]. Although this ad hoc solution does produce an answer, we cannot assess its accuracy, since dimensional analysis has ignored rather than solved the problem. To get a dimensionally correct answer, one could try choosing µ = R −1 or µ = R −1 before doing the R integral, or doing the Rintegral before choosing R at all. None of these attempts are consistent with perturbation theory, and in any case they all give a divergent answer.
An alternative approach that is discussed in the literature invokes constrained instantons [57][58][59][60]. The idea of constrained instantons is to fix R by demanding that some operator have a given expectation value, such as φ 3 = Λ 3 for some Λ. Fixing the scale in this way merely swaps the R integral for an integral over Λ −1 , and the problem is still there.
One might hope these infinities from J d and the R integral would cancel, but they do not (and should not). From a physical point of view, unless something makes the rate to produce different sizes bubbles different, the net decay rate should be infinite. It is quantum corrections which break the scale invariance, but the Jacobian is determined by the bounce from the classical theory. We will explain how to deal with the IR divergent integral over R and what scale the couplings are evaluated at in Section 3.2 after we have solved the Jacobian problem in the next section.

Solving the Jacobian Problem
The Jacobian is infinite because the dilatation fluctuation φ d is not normalizable. Of course, the path integral is basis independent; changing the normalization of a fluctuation N d φ d → φ d in the expansion in Eq. (2.13) can be compensated for by rescaling ξ j → 1 N d in the path integral measure. The problem with having a non-normalizable zero mode is that one cannot see clearly what happens when it is removed in computing det . In fact, the infinite Jacobian is secretly compensated by an infinity in det (see Appendix A). Here, we cleanly resolve the Jacobian problem by choosing a judicious basis in which the numerator and denominator path integrals can be computed exactly.
What basis would allow us to diagonalize fluctuations around the bounce and around the false vacuum simultaneously? Taking eigenfunctions of for the scalar fluctuations around the bounce and eigenfunctions of for fluctuations around the false vacuum. Note that even for an arbitrary potential these operators will differ only by a constant and therefore have the same eigenfunctions.
An important feature of eigenfunctions of these operators is the inner product by which they are orthogonal. Similarly to Eq.
so that eigenfunctions are orthogonal according to Since we are using the same basis for both path integrals, we can normalize the eigenfunctions however we like. For example, we could choose φ j |φ j V = 2π, and indeed even the dilatation mode will be normalizable according to this metric. Furthermore, this basis still lets us evaluate the path integral, since Thus the path integrals are still Gaussian in the fluctuations. The integral over a fluctuation normalized with φ j |φ j V = 2π then gives the usual factor of 1 λ j . Note that these observations apply to any theory, not just a scale invariant one: one can always simultaneously diagonalize fluctuations around the bounce and fluctuations around the false vacuum.
Now we restrict to the scale-invariant case, with V (φ) = 1 4 λφ 4 . Explicitly, our eigenfunctions should satisfy and be orthogonal with respect to the inner product Remarkably, we can find the solutions in closed form. For x 0 = 0 in Eq. (3.1), they are with P m l (x) the associated Legendre polynomial and Y slm (α, θ, φ) are the 3D spherical harmonics: These spherical harmonics satisfy and are normalized as The full eigenfunctions in Eq. (3.15) satisfy Eq. (3.13), i.e.
There are 5 modes with n = 1, with λ φ 1 = 0. The spherically symmetric one is This is proportional to the dilatation mode: φ 1000 = − −λ 4π Rφ d . The other n = 1 modes, which also have λ 1 = 0, are the zero modes for translations.
The modes with n > 1 are not particularly interesting: (3.28) and so on.
Since the dilation and translation modes are both normalizable, computing the Jacobian is straightforward: The Jacobian for the translation modes with this metric differs from Eq. (2.15): Note the important factors of R in both Jacobians -these are expected by dimensional analysis but obscure without the rescaling (cf. Appendix A). So we find Note that all the eigenvalues of O φ and O φ are dimensionless, so this expression has the correct units. We have shown that by rescaling the operators for fluctuations around the bounce and the false vacuum, the natural basis for field fluctuations changes, and the Jacobian for going between this basis and the basis containing a collective coordinate for dilatations is finite. Since the final result in Eq. (3.31) should be independent of this rescaling, there must be something that compensates for the infinite Jacobian if we do not rescale. In Appendix A we show that in fact without rescaling det is infinite as well.

Solving the Scale Invariance Problem
The next problem is that the integral over R in Eq. (3.31) is infinite. Even without evaluating the functional determinants, we can determine the R dependence of the integrand in Eq. (3.31) completely by exploiting renormalization group invariance of Γ. To see this, and to resolve the infrared divergence issue, it is critical to be consistent in power counting the loop expansion, or equivalently, orders of . A similar consistency was essential to resolve the gauge invariance problem of the ground state energy density in [17,38]. In the following, we insert appropriate factors of . Powers of will always correspond to powers of couplings such as λ in this scalar field theory or g 2 in a gauge theory.
To leading order (LO) in , the rate is determined entirely by the exponential factor in Eq. (3.31). Expanding this factor out explicitly we have where λ(µ) is the MS coupling at the scale µ and · · · refer to the rest of Eq. (3.31). Everything after the exponential comes from a one loop calculation and is subleading in . It is commonly said that the leading order prediction for the rate is Γ/V = e 8π 2 3λ(µ) . However, such a claim does not really make sense -not only is this equation dimensionally inconsistent, there is no indication at what scale µ to choose -so it is really no prediction at all. Indeed, the leading prediction must start at one loop. And, as we will see the leading prediction actually involves terms at two loops and higher. We will refer to the leading finite prediction with correct units as the NLO rate. Now, Γ is physical, so µ d dµ Γ = 0. This implies that the implicit µ-dependence of λ(µ) must be compensated by explicit µ-dependence in the NLO contribution. In turn, the µ dependence of λ(µ) is fixed by its RGE. Thus we know the exact µ dependence of the integrand to one-loop-higher order than we know the µ-independent part. By dimensional analysis, the only scale around to compensate µ is R, and therefore we also know the full R dependence of Eq. (3.32) to one loop: At one loop, the terms in (· · · ) have no explicit dependence on µ, by RG invariance, and therefore no dependence on R either, by dimensional analysis. Here β(µ) is the β-function coefficient in the RGE for λ, µ d dµ λ = β(µ). Now we see clearly the IR divergence problem. All of the R dependence in the one loop rate is explicit and the integral over R is infinite.
The only hope is for two-loop and higher-order contributions to come in and resolve the infinity. At first pass, this seems impossible, simply by counting factors of : terms in (· · · ) at two loops and higher are necessarily suppressed compared to the terms we have written. The resolution is that after the integral, superleading dependence is generated, as we will now see.
First of all, let us assume the MS coupling λ(µ) has a minimum at some scale µ = µ , so β(µ ) = 0. If this is not true, then the running coupling λ(µ) is unbounded from below and rate is actually infinite. In fact, in this quartic scalar field theory, λ(µ) is monotonic, so we are going to have to assume there are other fields in the theory to continue. For a more general theory, we can perform the path integral over all other fields around the bounce, leading to a decay rate formula of exactly the same form as Eq. (3.33), but with the βfunction for λ depending all all other couplings in the theory. In this case, β(µ) can vanish, as for example it does in the Standard Model for the Higgs quartic at the scale µ ∼ 10 17 GeV. 1 Since Γ is independent of µ we are free to choose µ = µ . Let us do so. Then the exponential in Eq. (3.33) has no R dependence (since β(µ ) = 0) at all and the integral is surely infinite. With µ = µ , the leading R dependence in the exponential factor comes in at two loops, and has the form and β 0 = β (µ ) using the one loop βfunction coefficients only (cf. Eq. (6.13) for its SM expression). At two loops, there is an additional single log term in the exponent scaling like ln(µ R). This will contribute subleading in after the integral so we have dropped it. Now we observe that since β 0 > 0 and λ < 0 at the minimum, the integral over R is finite: Note that this contribution is parametrically more important as → 0 than the one-loop correction which scales like 0 . Indeed, Γ 2 blows up as → 0, as it must to reproduce the divergence of the one loop integral. Thus, even though the two loop result is formally higher order in , we cannot justify expanding the exponential to provide only corrections to (· · · ) in Eq. (3.33). Also note that the divergence returns if β = 0. Thus the scale invariance is not regulated by just dimensional transmutation (i.e. by β = 0) but requires in addition that the β-function have a minimum.
A natural concern is that since the two-loop result parametrically dominates over the oneloop result as → 0, the three-loop result might dominate over two loops, and so on. To see if this happens, we examine possible terms in the exponent, as allowed by RG invariance. At each logarithmic order, the coefficient of ln n Rµ is of order n−1 plus terms suppressed by additional factors of . That is, we have ln Rµ , ln 2 Rµ , 2 ln 3 Rµ and so on. Using the two-loop term to set up a Gaussian around which we perform a saddle-point approximation in , we find a generic term becomes The · · · are all terms subleading in . So we see that in fact three-loop and higher-order contributions are more important than the one-and two-loop terms, but all terms at three loops and beyond are the same order in . One might worry that since the n = 3 term is more important than the n = 2 term, the saddle point approximation cannot be justified. Note however that expanding the exponential of the n = 3 term to next order gives a term scaling like 4 ln 6 Rµ . This term is subleading by a factor of to the terms we keep from the expansion of the n = 6 term. The same justification explains why we can ignore R dependence coming from the RG invariance of the non-logarithmic one-loop terms; these are also subleading in .
Since an infinite number of terms are relevant, we have to sum the series. Fortunately, this is possible since all of these terms depend on only the leading order β-function coefficients. In a pure scalar field theory, the one loop RGE is easy to solve in closed form. In a general theory, this will not be true. In any theory, we can write Note that we now sum over n ≥ 2 since λ (µ ) = 0 by definition of µ and, matching our previous notation, κ 2 = β 0 . The κ n coefficients are all determined by the β-functions in the theory evaluated at one loop. The terms denoted · · · all depend on β-function coefficients beyond one loop and are contribute to subleading order in to the final answer, so we drop them. Eq. (3.37) represents a perturbative solution of the coupled RGEs which is always possible to work out order-by-order in the couplings. It implies Here the · · · have terms of the same logarithmic order but subdominant in compared to the terms we keep. For the integrand to be RGE invariant, we know the rate can be written as To repeat, there are corrections to this exponent that are subleading in both from the higher-order terms in Eq. (3.37) and from those fixed by the RG invariance of one-loop (or higher) non-logarithmic terms in the integrand. All of these terms necessarily make subleading contributions to the rate. Now, every term in the exponent in Eq. (3.41) is proportional to a−1 ln a (µ R) for some a. There is only one term with a = 2, corresponding to n = 2, m = 1. This term generates a factor of Γ 2 after integration, as in Eq. (3.35). All the other terms, including the cross terms, have a > 3 and contribute to the same order in after integration by Eq. (3.36). We can perform the Gaussian integrals over all n and m by adding and subtracting the n = 2, m = 1 term. Then, as in Eq. (3.36), we get The 1 is subleading in , so we can drop it. The geometric series is easily resummed, giving Finally, we can express the answer in terms of the running coupling. The result is In this expression λ 1-loop (μ) means solve the coupled RGEs using the one loop β-function coefficients only and evaluate at the scalê For example, in a complex scalar theory,μ = µ exp(− π 2 6λ ). Note that for small coupling, the scalesμ and λ can be very far apart. Keep in mind however that all the couplings in λ 1-loop (μ) are evaluated at µ , so this resummation does not indicate sensitivity to high scales; it is merely shorthand for a series of terms all of the same order and the couplings.
Putting everything together and resetting = 1, we find Here the extra determinants come from integrating over fluctuations of fields other that φ in the theory around the bounce and false vacuum backgrounds. As long as µ exists (meaning λ(µ) has a minimum), this is a finite expression derived with consistent power counting. All of the singularities associated with scale invariance have been completely resolved. Finally, we point out that one does not have to choose µ = µ . For µ = µ there are terms linear in ln(µR) in the exponent proportional to β(µ), which generate a slew of additional terms in Eq. (3.46). For example, Eq. (3.35) becomes with β 0 the one loop β-function coefficient for λ. The general expression including terms like this can be used to calculate the scale uncertainty on the final prediction.

19
With the divergences associated with scale invariance understood, we now proceed to evaluating the functional determinants. We will be evaluating determinants for a number of cases: real scalar fluctuations, Goldstone boson fluctuations, gauge boson fluctuations (in general gauges), and Dirac fermion fluctuations. These are similar enough that is helpful to work out some results first that we can then apply to the different examples of interest.
In this section we consider the general operator and will evaluate Comparing to Eq.

Regularized Sum
The key to calculating D(x) in Eq. (4.2) is that we know the spectrum exactly. Although we know it exactly in d dimensions [63,64], regularizing the eigenvalues does not necessarily correspond to a well-understood subtraction scheme. Instead, we will work in 4 dimensions and remove the divergences using Feynman diagrams. Defining we know the spectrum of O(x) exactly, as in Section 3.1; Then the determinant is where the degeneracies d n are in Eq. (3.24). This sum is UV divergent at large n. We regularize the sum by subtracting the terms of order x and x 2 and then we will add those terms back in through dimensionally-regularized Feynman diagrams. Expanding at small x, we find Then we can perform the sum. That is, we compute and ζ (−1) ≈ −0.165 is the derivative of the ζ-function at −1.

Divergent Parts
To the subtracted part, we must add in dimensionally-regularized MS-subtracted divergent contributions. The subtractions were determined by removing the terms to second order in x. These terms can be reproduced by computing contributions to second order in x to the effective action using Feynman diagrams. The Euclidean action whose second variation gives We want to treat the mass term, proportional to x, as an interaction to compute the divergent contribution to the effective action.
To compute Feynman diagrams in Euclidean space, we expand e −S . The − sign in front of S affects all the Feynman rules, and Feynman diagrams produce contributions to −S eff ; that is, in Euclidean space −1 serves the role that the i prefactor does in Minkowski space. Thus, the interaction Feynman rule is Then we can perform the sum. That is, we compute and ζ (−1) ≈ −0.165 is the derivative of the ζ-function at −1.

Divergent Parts
To the subtracted part, we must add in dimensionally-regularized MS-subtracted divergent contributions. The subtractions were determined by removing the terms to second order in x. These terms can be reproduced by computing contributions to second order in x to the effective action using Feynman diagrams. The Euclidean action whose second variation gives We want to treat the mass term, proportional to x, as an interaction to compute the divergent contribution to the effective action.
To compute Feynman diagrams in Euclidean space, we expand e −S . The − sign in front of S affects all the Feynman rules, and Feynman diagrams produce contributions to −S eff ; that is, in Euclidean space −1 serves the role that the i prefactor does in Minkowski space. Thus, the interaction Feynman rule is In our notation solid lines are propagating φ fields and dashed lines are background field insertions. The injected momentum is distributed according to the Fourier transform of the bounce-squared [11]: In our notation solid lines are propagating φ fields and dashed lines are background field insertions. The injected momentum is distributed according to the Fourier transform of the bounce-squared [11]: At order x, there is only one graph, a tadpole Here, the 1 2 in front is a symmetry factor. This graph is scaleless and vanishes in dimensional regularization. Note that in 4D, the graph is quadratically divergent, in agreement with the O(n) term in Eq. (4.6) which is quadratically divergent when summed over n.
There is one graph with two x insertions: Here the 1 4 is the symmetry factor. The k integral can be done first giving Then we use to find that and therefore, rescaling Note that the integrals in Eqs. (4.16) and (4.17) are proportional to S[φ b ]. This is not a coincidence, as the divergences must be canceled by renormalizing λ in the tree-level action. 2 The full result in MS for a bosonic functional determinant this combined with Eq. (4.7) with S fin (x) in Eq. (4.8).
2 Technically, one should do the integrals in Eqs. (4.16) and (4.17) in d dimensions, using a d-dimensional bounce, generating O(ε) terms which contribute additional finite parts to Eq. (4.19). However, these finite parts must exactly cancel the finite contributions from the d-dimensional action renormalizing λ, since both multiply the same 1 ε terms. Thus we can ignore both and use the 4D integrals as written.

22
Here, the 1 2 in front is a symmetry factor. This graph is scaleless and vanishes in dimensional regularization. Note that in 4D, the graph is quadratically divergent, in agreement with the O(n) term in Eq. (4.6) which is quadratically divergent when summed over n.
There is one graph with two x insertions: At order x, there is only one graph, a tadpole Here, the 1 2 in front is a symmetry factor. This graph is scaleless and vanishes in dimensional regularization. Note that in 4D, the graph is quadratically divergent, in agreement with the O(n) term in Eq. (4.6) which is quadratically divergent when summed over n.
There is one graph with two x insertions: Here the 1 4 is the symmetry factor. The k integral can be done first giving Then we use to find that and therefore, rescaling Note that the integrals in Eqs. (4.16) and (4.17) are proportional to S[φ b ]. This is not a coincidence, as the divergences must be canceled by renormalizing λ in the tree-level action. 2 The full result in MS for a bosonic functional determinant this combined with Eq. (4.7) with S fin (x) in Eq. (4.8).
2 Technically, one should do the integrals in Eqs. (4.16) and (4.17) in d dimensions, using a d-dimensional bounce, generating O(ε) terms which contribute additional finite parts to Eq. (4.19). However, these finite parts must exactly cancel the finite contributions from the d-dimensional action renormalizing λ, since both multiply the same 1 ε terms. Thus we can ignore both and use the 4D integrals as written.

22
Here the 1 4 is the symmetry factor. The k integral can be done first giving Then we use to find that and therefore, rescaling Note that the integrals in Eqs. (4.16) and (4.17) are proportional to S[φ b ]. This is not a coincidence, as the divergences must be canceled by renormalizing λ in the tree-level action. 2 The full result in MS for a bosonic functional determinant this combined with Eq. (4.7) with S fin (x) in Eq. (4.8).

Angular Momentum Decomposition
We will also find it helpful to do the above sum in a different order, summing over n first exactly and then regularizing the sum over angular momentum modes s. The operators whose determinants we are trying to calculate are spherically symmetric. Thus their eigenfunctions are separable and can be written as The 4D Laplacian then reduces to a 1D operator, and there is a (1 + s) 2 -fold degeneracy for each s.

23
The appearance of the digamma function ψ(z) makes this subtraction more complicated than the subtraction in Eq. (4.6). Note that ψ (s + 1) ∼ 1 s at large s so there is a logarithmic divergence encoded in this expression. Performing the sum, we find in exact agreement with Eq. (4.8).
For some cases, the s = 0 modes will be absent. For these cases, it is helpful to have a form for the finite contribution with the s = 0 modes removed. To that end, we define

The Gelfand-Yaglom Method
There is a very powerful way of computing functional determinants, that does not require knowing the exact spectrum of the operators, called the Gelfand-Yaglom method [78]. Reviews and derivations of the method can be found in [40,45,[79][80][81], here we just summarize its application to the scale-invariant potential of interest in this paper. The Gelfand-Yaglom method says that functional determinants can be calculated by finding zero modes of the operators and evaluating their asymptotic behavior. For example, for 1-dimensional operators, like the ratio R s in Eq. (4.27), the method says that where φ(∞) and φ(0) are shorthand for the limits in this equation and the functions φ s The boundary conditions are that the functions be regular at r = 0. Note how powerful the method is: instead of finding an infinite number of eigenvalues and taking their product, we simply solve the differential equations in Eq. (4.32), which can be done numerically, and evaluate the solutions as r → ∞ and r → 0.
To see how the Gelfand-Yaglom method works, consider the real scalar case (x = −1) first. First of all, we already know the answer from Eq. (4.27): To use the Gelfand-Yaglom method, we find the exact solution to Eq. (4.32) regular at r = 0. It is The regular solution to Eq. (4.32) that reduces to φ s 0 at small r is Thus Gelfand-Yaglom predicts that in exact agreement with Eq. (4.33). Note that for n = 1, Eq. (4.33) gives R s = 0. This is because for n = 1 there are zero modes. Indeed, the n = 1, s = 0 zero mode is the dilatation mode and the n = 1, s = 1 modes are the translations. For the s = 1 case, the determinant ratio with eigenvalues removed is Similarly, R 0 = − 1 5 by multiplying the eigenvalues. To use Gelfand-Yaglom to calculate the zero eigenvalues, we shift the operator by order . That is, we replace Eq. (4.32) by Shifting the free-theory operator by is not necessary since all of its eigenvalues are nonzero. Note that the zero mode φ 11 is an eigenfunction of the shifted operator with eigenvalue . The function with eigenvalue 0 is φ 1, −1 (r) = φ n ,1 (r) with φ ns (r) in Eq. (4.25) and n ε = − 3 2 + 5 2 1 − 24 25 . Then, using φ 1 0 (r) = r, in agreement with Eq. (4.37). The s = 0 functional determinant can be computed in exactly the same way without any additional complication, finding R 0 = − 1 5 , in agreement with the direct calculation.

Functional Determinants
In this section we compute the functional determinants for the fluctuation of scalars, Goldstone bosons, vector bosons and Dirac fermions around the scale-invariant bounce configuration. We produce analytic formulas for all the cases. In the vector boson case, we check explicitly that the result is gauge invariant by using a generic value of ξ in Fermi gauges, and also show agreement between Fermi and R ξ gauges.

Real Scalars
The case of a single scalar field was introduced in Section 3. The Euclidean Lagrangian is Thus real scalar fluctuations correspond to the case studied in Section 4 with x = −1. For x = −1, the finite contribution from the sum over n ≥ 0 is singular, (S fin (x) in Eq. (4.8) is singular as x → −1). This is due to the zero modes at n = 1 corresponding dilatations and translations around the bounce. To compute the determinant with zero modes removed, we must first rescale the operator. We therefore define Recall from Section 3.1 that this rescaling allows the change to collective coordinates to have a finite Jacobian. The Jacobians for dilatation and translations are given in Eqs. (3.29) and (3.30).
To compute det O φ we must remove these modes from the sum in Eq. (4.7) and add in only the n = 1 contributions to the false vacuum fluctuations. We note the only zero mode for n = 1 in Eq. (4.7) is in which is singular at x = −1. Removing this term from the sum, we find a smooth limit as x → −1: Note that we should leave the x 2 terms in the subtraction at n = 1 to avoid overcounting, since these are included in the loops. Also note that we only remove the n = 1 term from det O φ , while the n = 1 term from det O φ is finite and should not be removed.

26
Combining with the divergent part from Eq. (4.19), we then have Here, O φ means the operator with φ b = 0, corresponding to fluctuations around the false vacuum. Note how the factors of γ E have dropped out. The remaining task is to renormalize. In MS the Z-factor for λ at one loop is The renormalized action on the bounce then becomes Combining with Eq. (5.6), we get

Complex Scalars and Global Symmetries
Next, we discuss the case when the false vacuum admits a continuous global symmetry that is spontaneously broken by the bounce. For concreteness, we take the simplest example, a field theory of a complex scalar field Φ with a global U (1) symmetry. The Euclidean Lagrangian density is where V (φ) = λ|Φ| 4 . We expand the field as With this normalization for a complex field, the bounce is the same as Eq. (3.3) and still satisfies − φ b + λφ 3 b = 0. Expanding around the bounce background to quadratic order, the scalar and Goldstone modes satisfy Both of these equations are special cases of Eq. (4.1), with x = −1 for φ and x = − 1 3 for G. The scalar fluctuations φ we have already discussed: there are 5 zero modes with n = 1 (since λ 1 (−1) = 0 with λ n (x) in Eq. (4.4)), corresponding to translations and dilatations.
The Jacobians for removing these zero modes in conformal coordinates are given in Eq. (3.29) and (3.30) and the functional determinant with zero modes removed is in Eq. (5.9) For G with x = − 1 3 there is a single zero mode (λ 0 (− 1 3 ) = 0) corresponding to phase rotations Φ → e iα Φ. The n = 0 mode has no degeneracy and the eigenfunction is as in Eq. (3.26). Infinitesimally along this direction, Φ = 1+iα √ 2 φ b which has the same action as Φ up to order α 2 .
As with φ, we have to remove the fluctuations along the zero-mode direction exactly using collective coordinates. With the measure determined by the operator O G = 1 To calculate the determinant with zero mode removed we follow the procedure in Section 5.1. The n = 0 mode in Eq. (4.7) contributes Note the singularity as x → − 1 3 . Removing this from the sum, we find a smooth limit as The full contribution from the Goldstone fluctuations is therefore where dθ = 2π gives the volume of U (1). For other gauge groups, the procedure is identical up to the group volume factor. Indeed, at quadratic order, all of the Goldstone boson directions decouple and the path integral over each direction gives a factor of Eq. (5.18) and a Jacobian. All that needs to be changed is the group volume factor. For SU (2), this is 2π 2 .
Putting the results for the Goldstone fluctuations together with the scalar fluctuations, we get for the complex scalar theory The R integral is would be cutoff by higher order effects if λ(µ) were bounded from below (which it is not in this theory). The UV divergence is canceled by the rernormalization of λ, as in the real scalar theory. In this case, the action on the bounce in renormalized perturbation theory is (cf. Eq. (5.8)) which is UV finite. We note for future reference that at each s the contribution to the functional determinant for the Goldstone modes follows from Eq. (4.27) We also note that for s = 0 there is a zero mode. Removing the zero mode we find R G 0 = 1. We also note that since n = 0 implies s = 0 if we remove the s = 0 modes, we automatically remove the zero modes. Thus, if we sum over only s > 0, the Goldstone contribution is simply

Vector Fields and Local Symmetries
Next, we discuss the contribution of gauge bosons to the decay rate. We continue the U (1) case, but now with Euclidean Lagrangian where V = λ|Φ| 4 as before and we expand Φ = 1 √ 2 (φ b + φ + iG) as in Eq. (5.11). For the gauge fixing term, we can consider the R ξ gauges, as in [11,19], where So that at quadratic order, While these gauges have some convenient features, particular for ξ = 1, they have a very serious drawback: they break the global U (1) symmetry. As a consequence there is no zero mode for the Goldstone fluctuations. Because of the missing zero mode, we are unable to reproduce the results of [11,19]. Although we are able to get sensible results in R ξ gauges using a numerical implementation of the the Gelfand-Yaglom method, we find Fermi gauges [40] more convenient for deriving exact analytic results. In Fermi gauges, the gauge fixing term is So that at quadratic order Fermi gauges leave the global U (1) symmetry of the Lagrangian intact (the action is invariant under Φ → e iα Φ, A µ → A µ ). Note that since the ghost Lagrangian is independent of the bounce, the functional determinant over ghosts normalized to the false vacuum is just 1.
In Fermi gauges, the equations of motion for A µ and G are coupled. At quadratic order Following [11,19,40], we then exploit the spherical symmetry, expanding A µ as (2) µ two independent generic vectors and Y slm (α, θ, φ) the 3D spherical harmonics in Eq. (3.16). In this basis, and writing G(x) = G(r)Y slm (α, θ, φ) the fluctuation operators decouple for each s, l, m and the resulting operators depend only on s. After some algebra (see [11] for some details), we find for ξ = 1 that the a S and a L modes couple to G, through the operator where the gauge-dependent piece is and ∆ s is in Eq.
For the full contribution of the transverse modes we must take the product over all s > 0. Note that there is no s = 0 contribution since Y 000 is constant and therefore the transverse contribution at s = 0 is absent in Eq. (5.33). Subtracting the s = 0 contributions from our previous results for the sum over all s, we then get for the transverse polarizations where S loops (x) is in Eq. (4.19) and S + fin (x) in Eq. (4.30). Explicitly: Note that there may be s = 0 contributions included in S loops (x). We will account for this by including the full gauge-invariant loop corrections for all modes at once when we discuss renormalization below.
Fluctuations with s = 0 For s = 0, Y 000 is constant and so the transverse and longitudinal modes decouple. In this case, only the scalar vector boson polarization and the Goldstone mode remain. The fluctuation operator is The corresponding operator with φ b = 0 is Note that M SG 0 has two zero modes regular at r = 0: The first zero mode corresponds to the global U (1) invariance we saw already in the g = 0 case. The two zero modes for M SG 0 regular at r = 0 are We need to remove the zero modes by going to collective coordinates, just as in Section 5.2. Since we do not know the eigenfunctinos of M SG 0 exactly we will use the Gelfand-Yaglom method. After rescaling our operator we add to it a shift of order . Then we need to compute We can compute the zero modes of −1 3λφ 2 b M SG 0 + ε · 1 perturbatively. Since only Ψ 1 goes to zero at r → ∞, we only need the corrections to it. So we expand, following [40], We can find the solution to this differential equation exactly by integration, following [40]: The result is Note that the result is gauge invariant, and its (trivial) g → 0 limit agrees with R G 0 = 1 computed at the end of Section 5.2. where det Ψ = det(Ψ j i ) where Ψ i are the 3 solutions and Ψ j i are the components of those solutions. Here and in the following, when we write Ψ(0) or Ψ(∞) we mean the leading behavior as r → 0 or r → ∞ respectively.
The functions of Ψ are easy to find. They are So det Ψ(r) = 2s(s + sξ + 2ξ)r 3s (5.54) Note that the ξ dependence of det Ψ is only in the normalization, so it will drop out in the ratio of the determinant at r = 0 and r = ∞.
To find the solutions Ψ, as discussed in [40], an immensely useful observation is that they can be expressed in terms of 3 auxiliary functions η, χ and ζ as where the auxiliary functions satisfy We define Ψ 1 as the solution with ζ = η = 0, Ψ 2 is the solution with ζ = 0 and η = 0 and Ψ 3 as the solution with ζ = 0. Note that only Ψ 3 can be gauge-dependent. The exact form of Ψ 1 is easy to find. With ζ = η = 0 we find χ = r s and so For Ψ 2 which has ζ = 0 but η = 0, we can solve Eq. (5.57) exactly. We find that the non-zero solution regular at r = 0 is with κ in Eq. (5.40). At small and large r Note the reappearance of the ratio in Eq. (5.39). Now given η 2 , we can solve for χ using Eq. (5.56). Conveniently, we do not need the full solution for all r, only its small r and large r behavior. Eq. (5.56) simplifies in these limits and we find To these we could add a homogeneous solution of the form χ = r s . However, this is exactly the Ψ 1 solution which is orthogonal, so adding a Ψ 1 component to Ψ 2 will not affect the functional determinant. Dropping the homogeneous solutions is extremely important -it is the essential simplification allowed by using these auxiliary functions. Using the limiting forms of η 2 and χ 2 , following the procedure outlined in [40], we find Here we have written only the terms that contribute at leading non-vanishing order to the determinant.
For Ψ 3 , defined to have ζ = 0, we can solve Eq. (5.58) exactly for ζ = r s . Proceeding as for Ψ 2 , we find Putting these solutions together we find where R AG s means all the gauge boson (A µ ), ghost, and Goldstone boson contributions are included. We have shown this result to be manifestly gauge invariant (ξ-independent in Fermi gauges). We also checked through a numerical implementation of the Gelfand-Yaglom method that the same formula emergences in R ξ gauges for each s.
The full functional determinant requires summing over s. We note that at large s, thus there are quadratic, linear and logarithmic divergences in the sum.

Renormalization
To regulate the sum over s, we will subtract the divergent terms and add in dimensionally regulated Feynman diagrams, as explained in Section 4. An important cross check on the result is that the UV divergences should cancel those from −S[φ b ] using the renormalized coupling. In scalar QED, the one loop Z-factor for λ is The action on the bounce then becomes Thus we need the UV divergences in Eq. (5.73) to be matched by the functional determinant over scalar, gauge and Goldstone modes.
To proceed, we want to compute the divergent contributions with Feynman diagrams in d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.30) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.73) as expected.
Fortunately, we can compute the regularized contribution in any gauge. Indeed we have checked through a numerical implementation of the Gelfand-Yaglom method that our result for the Fermi-gauge functional determinant R AG s in Eq. (5.70) is reproduced in R ξ gauges for any s. In R ξ gauges, for s = 0, the transverse and longitudinal model decouple, and we are left with the scalar, Goldstone and ghost fluctuations. The ghost contribution in R ξ gauges corresponds to x = − g 2 λ (like the transverse modes in Fermi gauges). At s = 0 we have checked numerically that this ghost contribution exactly cancels the contribution of the scalar and Goldstone contributions to R AG 0 , giving R AG 0 = 1 in R ξ gauges in agreement with our analytic results in Fermi gauges.
In R ξ gauge with ξ = 1 with Lagrangian in Eq. (5.28), the Feynman rules are d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected. Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction Then we find The sum, after taking µ → d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected. Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction Then we find The sum, after taking µ → d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected. Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction Then we find The sum, after taking µ →

74) In addition, in R ξ gauge the ghosts do not decouple and have an interaction
To proceed, we want to compute the divergent contributions with Feynman diagrams in d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected.
Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction The sum, after taking µ → To proceed, we want to compute the divergent contributions with Feynman diagrams in d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected.
Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction The sum, after taking µ → To proceed, we want to compute the divergent contributions with Feynman diagrams in d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected.
Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction The sum, after taking µ → To proceed, we want to compute the divergent contributions with Feynman diagrams in d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected.
Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction The sum, after taking µ → To proceed, we want to compute the divergent contributions with Feynman diagrams in d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected.
Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction The sum, after taking µ → The sum, after taking µ → Note that the divergent terms agree with those in Fermi gauge, Eq. (B.16), and when the scalar contribution is added (with divergence 3 2ε ), the poles exactly cancel those in Eq. (5.73).
To perform the subtraction, we need to compute the contribution to the functional determinants in 4D from terms to second order in the couplings. For s = 0, since R AG 0 = 1, there are no second-order contributions (or contributions at any order). For s > 0, we note that in R ξ gauge the transverse modes have the same quadratic fluctuations as the ghosts and the two contributions exactly cancel. For the other photon polarizations and Goldstones, the fluctuation matrix with ξ = 1 is Changing basis, following [11] we find a convenient almost diagonal form In this form, we see that if we turn off the off-diagonal couplings, each diagonal term is a 1D operator and the exact result can then be read off, using Eq. (4.27): The required subtractions to second order in the diagonal interactions then come from the expansions of these function to second order in their arguments: ln R SLG,R ξ s, diag sub = 2 (s 3 + 4s 2 + 4s + 2) s 2 (s + 1) 2 − 4ψ (s) + 2 (s 4 + 9s 3 + 20s 2 + 18s + 8) s 2 (s + 1) 2 (s + 2) − 8ψ (s) g 2 λ + 2(2s + 5) (3s 4 + 12s 3 + 18s 2 + 12s + 4) with ψ(z) = Γ (z) Γ(z) the digamma function. The remaining required subtractions involve the off diagonal couplings in Eq. (5.82). Since these couplings are linear in g there are no contributions to first order in g (corresponding to no diagrams with one g insertion). The contributions to second order in g can then be computed turning all the diagonal interactions (mass term) to zero. We do this with the Gelfand-Yaglom method perturbative in g, following a similar procedure to the one used in Section 5. 3   At fixed s the determinant is finite, but there is a UV divergence when summing over s at large s. Subtracting the leading large s behavior lets us do the sum, giving S fin . We regulate the UV divergences in d dimensions and compute the loops and counterterms in MS using Feynman gauge, giving S loops . The loops give the the correct large s behavior, but have a different finite part from what we subtracted to get S fin , so we have to add back in the difference between the loops and our subtraction. This is S diff . Finally, there are IR divergences form the zero modes at s = 0. We remove these from the functional determinant, finding R AG 0 = 1. Since the s = 0 contribution is 1, there is no contribution at s = 0 to the UV subtractions we did. Thus we simply have to remove the s = 0 part from all of our finite sums. This changes S fin to S + fin giving our final answer. 39

Fermions
Next, let us consider the addition of Dirac fermions. The Euclidean Lagrangian for a real scalar interacting with a Dirac fermion is Around the bounce configuration, φ = φ b , the fermion fluctuation operator is To calculate the determinant of this operator, we expand in a basis of half-integer spin spherical harmonics. Including angular momentum, Dirac spinors transform in the direct sum of (k + 1 2 , k) and (k, k + 1 2 ) representations of the Lorentz algebra su(2) ⊗ su(2) ∼ = o(4). In a particular representation of the Euclidean Dirac algebra, the half-integer spherical harmonics take the form of hypergeometric functions (see Appendix A of [82]). Expanding in this basis, M ψ reduces to a form which depends only on the radial coordinate r and only on two components of the Dirac spinor To match the literature, the first matrix here corresponds to Eq. (3.17) of [82] with K → k and the second to (3.18) with L = K − 1 2 → k. The multiplicity of the (a, b) representation is (2a + 1)(2b + 1), so we have Next, we reduce the product of these operators to a quadratic form by conjugating with the unitary matrices U = diag( 1 r 3/2 , 1 r 3/2 ) and V = diag( 1 r 3/2 , − 1 r 3/2 ): This simplifies slightly by writing k = j 2 − 3 4 . Then j = 2k + 3 2 = 3 2 , 5 2 , 7 2 , · · · , the multiplicity becomes j 2 − 1 4 and in agreement with [11]. The matrix for fluctuations around the false vacuum, M jψ ψ is the same as this one with φ b = 0.
Following similar techniques to those described in previous sections, we find, For the subtractions, it is helpful to have the result for the determinant when the off-diagonal terms in M jψ ψ are set to zero. When the matrix is diagonal, the fermion case is a special case of the general formula in Section 4. The result is with R s (x) given in Eq. (4.27).
The subtractions required to sum over j are given by the expansions of Rψ ψ j to second order in the diagonal couplings and second order in the off-diagonal couplings: Using this, we find ∞ j=± 3 2 ,± 5 2 ,··· This function is real for imaginary z and contains only even powers of z when expanded around z = 0. The UV divergent part is added back in through a dimensionally regulated calculation quadratic in the interactions. We do this by evaluating the functional determinant as in [11]: are the interactions from Eq. (5.97), and the subscript W 2 indicates that we are truncating the expansion in W to second order. The traces can be rewritten in momentum space The single W -trace is zero in dimensional regularization, and we can evaluate the other one using W (q) = y 2 As a check, we note that the action on the bounce in this Yukawa theory in terms of the renormalized couplings is The UV divergences in this action exactly cancel those in Eq. (5.111) and Eq. (5.6).

Summary of Results
To provide a convenient reference, we summarize here the main results of this section: the functional determinants for scalars, vectors and fermions. For a real scalar, there are zero modes corresponding to dilatations and translations, with Jacobian factors given in Eqs. (3.29) and (3.30): The fluctuation determinant with zero modes removed is in Eq. (5.6): In collective coordinates, we must integrate over d 4 x dR.
For a complex scalar field, there is a global U (1) invariance spontaneously broken by the bounce. There is a zero mode corresponding to phase rotations. The Jacobian for changing to collective coordinates is given in Eq. (5.15) and the fluctuation operator for the Goldstone bosons is in Eq. (5.18) In collective coordinates this must be integrated over the volume V = 2π of U (1). For a U (1) gauge theory with a complex scalar, namely the Coleman-Weinberg model, the dilatation, translation and phase rotation modes are still present. The functional determinant over gauge and Goldstone fluctutations with zero modes removed is in Eq. (5.90)

Vacuum Stability in the Standard Model
Now we have all the ingredients necessary to compute the next-to-leading order decay rate in the Standard Model. The relevant part of the Standard Model Lagrangian is where H is the Higgs doublet, H = iσ 2 H, W µ are the SU (2) gauge bosons, B µ is the hypercharge gauge boson, Q is the 3rd generation left-handed quark doublet, and t R and b R are the right-handed top and bottom quarks. Contributions from other fermions are negligible and gluons have no effect at next-to-leading order. We have set the Higgs mass parameter m 2 to zero; m 2 = 0 corrections will be discussed in Section 7. From this Lagrangian we see that there are only five parameters relevant to the NLO decay rate: λ, y t , y b and the SU (2)×U (1) couplings g and g . All these parameters depend on scale. As explained in Section 3.2, for a consistent power counting the tunneling calculation has to be done near the scale µ where β λ (µ ) = 0. In the SM this scale is µ ≈ 10 17 GeV. The five parameters are determined at much lower scales, µ ∼ 100 GeV (or µ ∼ m b for y b ). In determining the parameters matching conversions (also known as threshold corrections) are made from a physical scheme (like the pole-mass scheme where the W and Z masses are measured). The ingredients for this matching step are known at NNLO and depend on additional SM parameters, such as α s . After matching, one must run the couplings up to µ ∼ 10 17 GeV. The RG equations for this running are known at three or four loops and involve additional parameters like α s as well. To be clear, the goal of the matching and running is to get λ(µ ), y t (µ ), y b (µ ), g(µ ) and g (µ ) in MS. Thus, it is perfectly consistent to match at NNLO, run at three or four loops and compute the decay rate at NLO -each step is a separate well-defined calculation.

NLO Tunneling Rate Formula
To compute the NLO tunneling rate in the SM we need to combine the formulas from Section 5.
The bounce spontaneously breaks translation and scale invariance as well as SU (2) × U (1) Y → U (1) EM . The zero modes for translations and dilatations must be integrated over with collective coordinates with appropriate Jacobian factors. The three broken internal generators produce three zero modes which must be integrated over the volume of the broken gauge group. As we work to NLO only, we only need the action quadratic in the fluctuations around the bounce. For gauge bosons, this means the non-Abelian interactions are irrelevant and each gauge boson can be treated independently. Thus, each gauge group collective coordinate produces a factor of J G in Eq. (5.15) as for a U (1). In addition, since the U (1) representing electromagnetism is unbroken, fluctutations of the photon are the same around the false vacuum and the bounce and therefore do not contribute to the rate. We can therefore compute the gauge-boson fluctuations by integrating over W ± and Z boson fluctuations and their associated Goldstone bosons.
Resolving the integral over instanton size R through the technique described in Section 3.2 and using Eq. (3.46), the NLO rate formula in the SM is therefore This formula is valid for With SU (2), the group theory volume factor is 3 For the real Higgs scalar, the determinant with zero modes removed is in Eq. (5.114).
For gauge bosons, we note that the W and Z bosons couple to φ b with strengths g W = 2 m W v = g and g Z = 2 m Z v = g 2 + (g ) 2 respectively. Including the conventional factor of 1 2 normalizing Abelian versus non-Abelian generators, the gauge bosons and Goldstone fluctuations give the result summarized in Eq. (5.117) The top quark contributes as in Eq. (5.118) with a factor of N C = 3 for color The bottom quark contribution is identical with y t → y b and we omit y b for simplicity in the next set of formulas. The UV divergences from the product of these functional determinants in 4 − 2ε dimensions is These are exactly canceled by the renormalized tree-level action on the bounce In the SM, the one loop β function for λ is 12) This one loop β function is of course linearly related to the 1 ε poles in Eq. (6.11). The derivative of β λ , required for Eq. (6.2) is Note that although β λ is formally two-loop order, it depends only on one loop β-function coefficients. Thus in the consistent NLO calculation all that is needed is one loop results.

Absolute stability
Absolute stability means that Γ = 0 and our electroweak vacuum will never decay. A naive criterion for Γ = 0 is that λ < 0. That is β λ (µ ) = λ(µ ) = 0 (naive absolute stability) (6.14) This criterion has been been used in many treatments to establish the stability boundaries of phase space. For example, in [18,83], this boundary is fixed by λ(µ cri ) = β λ (µ cri ) = 0 (their µ cri is our µ ). We call this criterion "naive" because is not systematically improvable: it only depends on the running λ. For example, if λ is positive but very very small, the rate can still be nonzero due to loop corrections but the naive criterion would miss this possibility. 4 To compute the absolute-stability phase space boundary, a gauge-invariant systematicallyimprovable procedure was developed in [17,38]. The starting point is that absolute stability is equivalent to the electroweak vacuum being the absolute minimum of the effective potential. Although the exact value of the potential at a minimum is known to be gauge invariant [21,22], some care has to be exercised in extracting this minimum value in perturbation theory. Because the effective potential having a minimum requires tree-level and one-loop effects to be comparable, the standard loop power counting cannot be used to establish stability (it violates gauge-invariance). A self-consistent gauge-invariant procedure for establishing absolute stability was developed in [38]. Briefly, one starts with the leadingorder effective potential The first term of this potential is tree-level and rest comprises all of the one-loop corrections consistent with the power counting established in [17,38]. The point of the power counting is that since the one-loop contribution must overwhelm the tree-level contribution to turn the potential over, λ must be the size of the one loop corrections. Remarkably, one must impose this power counting consistently for gauge invariance to hold order-by-order in perturbation theory. The minimum of V LO is where the couplings satisfy λ = 1 256π 2 g 4 + 3g 4 + 2g 2 g 2 + 3(g 2 + g 2 ) 2 ln 4 g 2 + g 2 + 6g 4 ln 4 g 2 − 16N C y 4 t ln 2 y 2 t + 1 (6.16) We denote by µ X the MS renormalization-group scale where this equation holds. The NLO effective potential with this consistent power counting is then computed by combining oneloop, two-loop and an infinite set of higher-loop daisy diagrams. Since the stability bound with this procedure is gauge invariant (as checked explicitly in scalar QED in [38]), we can choose any gauge. Landau gauge (ξ = 0) is particularly convenient as all the daisy diagrams vanish. The NLO effective potential in Landau gauge is extracted from [14,15]. We present it in Appendix C for completeness.
One cannot be certain that our universe is absolutely stable, as quantum gravity or new physics coming in at an arbitrary high scale can open up new tunneling directions that can destabilize the universe [16,19,50,62,[85][86][87]. So there is no sensible way of estimating a lower bound on the lifetime of our vacuum including new physics. The best one can do is to put an upper bound on the lifetime, and the only question we can reasonably ask about new physics is at what scale, Λ NP , it could come in to stabalize our vacuum? That is, how strong would it have to be to raise the upper bound on the lifetime to make it absolutely stable? To determine this scale, we add to the effective potential a gauge-invariant operator This operator contributes to V LO and modifies the equation for µ X , Eq. (6.16). Then we ask for given SM couplings, what value of Λ NP will lift the minimum of V eff to zero. The curves for this condition in the SM are shown in Fig. 3.

Numerical Results
For numerical calculations, we take as inputs and α s (m Z ). These inputs are converted to MS at a scale µ 0 = m pole t using threshold corrections known to two loops in all SM couplings [14,15,83], including mixed strong/electroweak contributions, and partially to three and four loops in α s . The couplings are then run to high energy using the three-loop renormalization group equations with four-loop running included for α s [88][89][90]. All of these threshold and running calculations are conveniently performed using the MR package of Kniehl, Pikelner and Veretin [91].
The numerical values are taken from the 2017 Particle Data Group [92]. We take as inputs These uncertainties will be propagated through to the final results.
And, as needed for Eq. (6.2), The action on the bounce is The terms in Eq. (6.2) evaluate to  The first three uncertainties are from variation of m t , m h and α s respectively according to Eq. (6.19). The fourth uncertainty is theory uncertainty from varying the threshold matching scale µ thr = ξm pole t with 1 2 < ξ < 2 used in converting observables to MS and as the starting point for RGE evolution. The final uncertainty marked NNLO represents the unknown two loop contributions to the functional determinant around the bounce. We estimate this error by scale variation around µ by a factor of 1 2 or 2. Noting that the NLO tt functional determinant contributes in the exponent at around 3% of the tree-level bounce action; therefore our NNLO estimate of 7% compared to NLO seems reasonable.
The variations in the first line of Eq. (6.26) are not independent and the dependence of Γ on the masses and scales is highly non-linear. Nevertheless, since we can compute the effect on Γ for any combination of their variations, we can determine their total correlated effect on the rate. To do this, we maximize or minimize the rate over χ 2 = 1 hypersurface. We find that at 68% confidence 10 −1411 < Γ V GeV 4 < 10 −533 . The range of decay rates allowed at 95% confidence is e −5660 < Γ V GeV 4 < 10 −386 .
Thus, the lifetime of the Standard Model universe is That is, to 68% confidence, 10 102 < τ SM years < 10 321 . To 95% confidence 10 65 < τ SM years < 10 1383 . To be more clear about what the lifetime means, we can ask a related question: what is the probability that we would have seen a bubble of decaying universe by now? Using the space-time volume of our past lightcone [15], (V T ) light-cone = 0.15 Since the bubbles expand at the speed of light, chances are if we saw such a bubble we would have been destroyed by it; thus it is reassuring to find the probability of this happening to be exponentially small. The phase diagrams in the m t /m h and m t /α s planes are shown in Fig. 2. In these diagrams, the boundary between metastability and instability is fixed by P = 1, where P is the probability that a bubble of true vacuum should have formed without our past lightcone, as in Eq. (6.28). The boundary between metastability and instability is determined by the gauge-invariant consistent procedure detailed Section 6.2 (and in [17,38]). Although the absolute stability boundary is close to the condition λ = 0 in Eq. (6.14), it is systematically higher and a better fit to the curve for λ = −0.0013.
Varying one parameter holding the others fixed, we find that the range of m pole Absolute stability is currently excluded at 2.48σ, which translates to a one-sided confidence of 99.3%. To exclude absolute stability to the one-sided confidence for 3σ, the top quark mass uncertainty must be reduced below 250 MeV. Similarly for α s for a 3σ uncertainty must be less than ∆α s < 0.00025. The dashed lines Fig. 3 indicate the scale at which new physics operators at the scale Λ NP can stabilize the SM, added as in Eq. (6.17). Recall that because tunneling is a nonperturbative phenomenon, higher-dimension operators do not decouple: new physics at an arbitrary high scale can destabilize the SM my opening up new tunneling directions [19,50,62,[85][86][87]. To stabalize the SM, they have to be strong enough to lift the potential from negative to positive. In Fig. 3 we see that the density of Λ NP curves increases near the absolute stability line. This happens because the absolute stability region is necessarily insensitive to the addition of a positive operator.

Mass Corrections
One remaining technical detail is how to handle the fact that the Higgs potential in the Standard Model is not exactly scale invariant, since there is a finite mass term for the Higgs field. We saw in Section 3 that with a scale-invariant classical potential, quantum corrections naturally pick out the scale µ where λ(µ) is minimal so that the action is dominated by bounces of a size R = 1 µ . One hopes that because the Higgs mass parameter m ∼ 10 2 GeV is much much smaller than µ ∼ 10 17 GeV, the corrections to the decay rate from the mass term will be completely negligible. Although normally classical effects, like the Higgs mass term, dominate over quantum effects, in this case the quantum scale violation can be dominant since it scales as an inverse power of (see Eq. (3.35)). Despite this convincing logic, producing a quantitative estimate of the effect on the decay rate of a finite mass term is surprisingly challenging.

A Bound on the m 2 Correction
Consider the potential V (φ) = 1 2 m 2 φ 2 + 1 4 λφ 4 with λ < 0 and m > 0. Trying to solve the Euclidean equations of motion for this potential, φ + 3 r φ − m 2 φ − λφ 3 = 0, one quickly discovers that the only solution is φ = 0. There are many ways to see this [57,59,60,62] such as with Derrick's theorem [99]. An intuitive way is to use Coleman's trick of thinking of the solution the Euclidean equations of motion as a ball rolling with friction down a hill shaped like −V (φ) starting at φ(0) and ending at φ = 0 at "time" r = ∞. For m = 0, the potential is scale invariant, so no matter where the ball starts it will get to φ = 0 only asymptotically at infinite time; different starting points correspond to different R for the bounce solutions in Eq. (3.1). Now, when we add 1 2 m 2 φ 2 to the potential, it creates a depression in −V (φ) near φ = 0. Since for any R the bounces just barely got to φ = 0 at infinite time, adding even an infinitesimal depression prevents solutions to the equations of motion from ever reaching φ = 0. Thus there are no bounces when m 2 > 0. Assuming m is small compared to µ , one might think we can write φ = φ b + m 2 ∆φ + · · · and evaluate the corrections to the action perturbatively. Trying this, one immediately finds This behavior is due to the non-normalizabilty of φ b . Thus Γ ∼ e −S = 0 confirming that even an infinitesimal m 2 seems to prevent vacuum decay.
To understand this unintuitive result, let us consider the alternative, more physical, treatment of tunneling described in [49,50]. There, a formula for the tunneling rate was derived inspired by the understanding of tunneling in non-relativistic quantum mechanics. In quantum field theory, the exponential factor determining the decay rate along a path parameterized by φ( x, τ ) is the integral where the energy functional is [75,100,101] Assuming m is small compared to µ , one might think we can write φ = φ b + m 2 ∆φ + · · · and evaluate the corrections to the action perturbatively. Trying this, one immediately finds This behavior is due to the non-normalizabilty of φ b . Thus Γ ∼ e −S = 0 confirming that even an infinitesimal m 2 seems to prevent vacuum decay.
To understand this unintuitive result, let us consider the alternative, more physical, treatment of tunneling described in [49,50]. There, a formula for the tunneling rate was derived inspired by the understanding of tunneling in non-relativistic quantum mechanics. In quantum field theory, the exponential factor determining the decay rate along a path parameterized by φ( x, τ ) is the integral where the energy functional is [75,100,101] In Eq. (7.2) τ is the Euclidean time and s is the proper time, determined by ( ds dτ ) 2 = 2U [φ]. Using s gives a formula exactly like the WKB exponent formula dx 2V (x) in quantum mechanics, but now with a contribution from gradient energy.
With this formulation let us now revisit the perturbative solution. If we try to calculate Γ φ along the m = 0 bounce path φ b (r = √ x 2 + τ 2 ), we find The integral over the first term gives S[φ b ] = − 8π 2 3λ , as for m = 0. The second term, however, shifts U (τ ) along this path up by a finite positive amount at each τ , and the resulting integral over τ is infinite.
One fault of using the path through the φ b bounces in a non-scale-invariant potential is that conservation of energy is violated. Since bubbles are produced at rest, we can see this If the Euclidean equations of motion are satisfied, energy is conserved and this cannot happen.
Let us consider instead a path through field space of the form These Gaussian bubbles were previously introduced and discussed in [50]. Their energy is With this value for φ 0 , the partial width for decaying to Gaussian bubbles is This is finite. Moreover, it has a nonzero maximum at R = 0. We conclude that the exponent is bounded from above. In other words, The lower bound comes from m = 0 and the upper bound comes from the Gaussian bubbles as R → 0. We conclude that the rate is finite.

Constrained Instantons
Now that we know the rate is finite, what is it? The Gaussian bubbles are in fact far from producing the optimal path through field space, as we will see. What we would like to do is directly minimize dτ U [φ(τ )] over field configurations with φ(0) = 0. Luckily, we do not have to find the absolute minimum; if we find any path which gives a finite rate close enough to the m = 0 case that we can neglect it for the Standard Model, we can conclude that ignoring m when m µ is justified. One way to find a finite-action path through field space is through the constrained instanton approach [57,60]. In brief, the idea is that instead of minimizing the action absolutely, we find the minimum along some surface. For example, we can look along the surface where d 4 xφ n = k for some k and some n. This constraint can be imposed through a Lagrange multiplier by writing the action as Taking n = 2 or n = 4 does not produce anything helpful since the new term is just like one of the old ones. Taking n 5 is also unhelpful, since for a normalizable solution we need φ → 0 at large distance, but then φ n term is subdominant to the λφ 4 term which produced the non-normalizable mode in the first place. Thus n = 3 is our only hope. (See [59] for a thorough discussion of constraints on the constraints).
For n = 3, the procedure for producing a normalizable well-behaved constrained solution that reduces to φ b at m = 0 is discussed in [59] and [58]. The solution and Langrange multiplier can be expanded perturbatively in m 2 : It is helpful to write φ 2 = φ 2,a + φ 2,b with To determine the constant in Eq. (7.13) and the O(m 2 ) value of the Lagrange multiplier, we note that perturbative solution is not normalizable. This non-normalizability is easy to understand: a solution perturbative in m can never describe the asymptotic behavior for r 1 m . no matter how small m is. Indeed, a normalizable solution should have φ → 0 as r → ∞ and therefore match on to φ K (r) = K 0 1 mr K 1 (mr) which solves ( + m 2 )φ = 0. At large r, φ K (r) ∼ K 0 π 2 1 (mr) 3/2 e −mr which is exponentially suppressed and gets contribution from all orders in m. We choose K 0 = 8 −λ m 2 R so that to order m 0 , φ K (r) matches on to φ b (r) at large r. This fixes the value for φ 2,b in Eq. (7.13) to be This is the unique solution allowing φ 2 at large r to match the m 2 terms of φ K (r). The Lagrange multiplier in Eq. (7.13) is then σ = −λφ 2,b + O(m 4 ) so that φ satisfies the (constrained) equations of motion to order m 2 . Higher order terms can systematically computed guaranteeing exponential suppression at large r [58,59]. Given these results, we now need to check that the decay rate along the constrainedinstanton path is finite. If we work strictly to order m 2 we find the action gets corrected by This is not surprising as φ b + φ 2 is not normalizble. The key to getting a finite action is to be careful with the boundary term. Without integrating by parts we can write The equations of motion on the constraint surface imply φ − m 2 φ − λφ 3 − 3σφ 2 = 0 so we drop this term. We also drop the total derivative term using that the exact φ vanishes exponentially at infinity (although not at any fixed order in m 2 ). The remainder can then be evaluated perturbatively:

Comments on Constrained Instantons
Before concluding, we add some comments on the constrained instanton approach. First, we made no claim that the constrained instanton produces the exact decay rate, as there may be lower-action configurations satisfying different constraints. For situations where m / 1 R it may be important to look for other solutions and Eq. (7.19) cannot be used in such contexts.
Second, we never try to integrate over the value for the constraint k. Normally one sets k at the outset and solves for the Lagrange multiplier σ as a function of k. Here we have found σ by requiring the solution have finite action; k is then fixed by σ. Although at order m 2 , k = ∞, k should be finite when the full solution is used, since in the full solution φ dies off exponentially at large distance. In any case, we do not need k to compute the tunneling rate, as we have seen.
Third, although the constrained instanton approach is useful to understand m > 0, it does not help resolve the divergent integral over instanton size R for m = 0. Since the true minimum of the action is known when m = 0, the divergence must be resolved from higher order perturbative effects, as we explained in Section 3.2. There has been some confusion about this point in the literature [11].
Fourth, we note that there is an apparent contradiction that the Euler-Lagrange equations have no solution with m 2 > 0, but we have proven there is a finite, non-zero minimum to the action. The resolution is that Euler-Lagrange equations are derived dropping a boundary term, but the behavior of the solutions at infinity are critical to finding a correct minimum. As we have seen, to any finite order in m 2 , the boundary terms cannot be dropped, so one would never come upon a perturbative solution like ours using the Euler-Lagrange equations alone. The importance of the boundary behavior is emphasized and discussed at length in [57][58][59] as a motivation for the constrained instanton approach.
Finally, we have done the whole analysis here, following [58,59] for the case m 2 > 0. The case with m 2 < 0 is also interesting. For m 2 < 0, the energy function U [φ b (τ )] gets shifted down and, for mR > 1, the tunneling rate is in fact infinite: there is no barrier to tunneling (as in quantum mechanics with a potential like V = −x 2 ). This result is also wrong. The argument is flawed, since U [φ b (0)] = 0 just like for the m 2 > 0 case, so the proposed tunneling path violates energy conservation is not allowed. Tunneling should speed up, for m 2 < 0, but only by an amount suppressed by factors of m 2 R 2 . By analytic continuation we can still use Eq. (7.19); for m 2 < 0 S[φ] now has a small imaginary part, but this produces a tiny effect on the final result, since Γ ∼ Im(ie −S[φ] ).

Conclusions
In this paper we have produced the first complete calculation of the lifetime of the Standard Model. Previous treatments were incomplete in a number of ways. First, there was a long-standing problem of how to perform instanton calculations when scale-invariance is spontaneously broken. The problem is that in a classically-scale invariant theory, the integral over instanton size R is divergent at next-to-leading order (NLO). We showed that in fact there are contributions which seem higher-order in but which in fact dominate over the NLO contribution after the integral over R is performed. Including all the relevant terms, to all loop-order, we are able to integrate over instanton size exactly giving a finite result.
The second problem we resolved is also related to instanton size. Since fluctuations associated with changing the size R are unsuppressed, one has to allow for large deviations in field space. Changing to collective coordinates allows the integral over all R to be done, however, it generates an infinite Jacobian. We showed that this infinite Jacobian is in fact compensated by an infinity in the functional determinant previously missed. To handle the infinity and the zero, we employ a judicious operator rescaling inspired by a conformal mapping to the 4-sphere. We find the spectrum of the rescaled operators exactly and give an analytic formula for the Jacobian (now finite) as well as the functional determinant with zero modes removed (also finite now).
The third problem we resolved has to do with fluctuations of vector bosons around the instanton background. When a global internal symmetry is spontaneously broken there are additional zero modes. In previous treatments the Jacobian for going to collective coordinates for these symmetries was found to be infinite. We show that this infinity was an artifact of working in R ξ gauge where the symmetry is actually explicitly broken by the gauge-fixing. Instead we work in Fermi gauges, and using the same technique as for the dilatation zero mode, show that the Jacobian for internal symmetries is finite.
The next new result in our paper is a complete analytic computation of the functional determinant around the instanton background for real and complex scalar fields, vector bosons and fermions. Moreover, we showed that the final result is gauge-invariant (of the parameter ξ in Fermi gauges and between Fermi and R ξ gauges). For the scalars, the insight which allowed for these exact results was to use the exact spectrum known from the operator rescaling and mapping to the 4-sphere [63][64][65]. For the vector bosons, we exploited a remarkable simplification of the fluctuation equations discovered in [39,40]. These authors found that the equations that couple the scalar and longitudinally polarized gauge bosons with the Goldstone bosons can be written in terms of a set of simplified equations using auxiliary fields. Although the treatment in [39,40] assumed a mass term for the scalar, so that their results do not exactly apply to the case of the Standard Model, our treatment very closely parallels theirs.
Combining all our results together we produced a complete prediction for the lifetime of our metastable vacuum in the Standard Model. We find the lifetime to be τ SM = 10 161 +160 The enormous uncertainty in this number is roughly equal parts uncertainty on the top quark mass, uncertainty on the value of the strong-coupling constant α s and theory uncertainty from threshold corrections, that is, from matching between observable pole masses and MS parameters at the electroweak scale. The uncertainty from error on the Higgs boson mass is small as is, thankfully, uncertainty associated with the unknown NNLO corrections to the decay rate.
Phase diagrams in the m t /m h plane and the m t /α s plane are shown in Fig. 2. This figure indicates that the SM seems to sit in a peculiarly narrow swath of metastability in the phase space of top quark mass, Higgs boson mass, and strong-coupling constant. An important fact to keep in mind when interpreting this tuning is that phase diagram assumes no gravity and no physics beyond the SM. In fact, any arbitrarily high-scale physics can destabilize the SM by opening up new tunneling directions [19,50,62,[85][86][87]. Moreover, near the absolute stability boundary, operators at an arbitrarily high scale can also stabilize the SM, as can be seen from Fig. 3. For the SM, which appears not to be on the stability boundary, the relevant scale of new physics is around 10 13 GeV.
Because of the importance of the top quark mass, the Higgs boson mass and α s in determining stability, it is interesting to look at their allowed ranges. We find that, varying each parameter separatly, the bounds for the SM to lie in the metastability window are If we hope to rule out absolute stability to 3σ confidence, assuming nothing else changes, we would need ∆m pole t < 250 MeV or ∆α s (m Z ) < 0.00025. Finally, we note that the predicted lifetime 10 161 years, while enormously long, has an exponent of roughly the same order of magnitude as the current lifetime of the universe, 10 9 years. Indeed, the long lifetime of the SM is due to the fact that the Higgs quartic coupling has a minimum value of λ = −0.0138. If the minimum of the coupling were smaller, say λ = −0.1, then the SM lifetime which scales like exp( 8π 2 3λ ) would be only 10 −20 seconds! Furthermore, since the lifetime is finite and the universe infinite, there is likely a bubble of true vacuum already out there, far far away. It is sobering to envision this bubble, with its wall of negative energy, barreling towards us at the speed of light. It seems the long-term future of our universe is not going to be slow freezing due to cosmic acceleration but an abrupt collision with one of these bubble walls.
functional determinant should be independent of the operator rescaling, this infinity must be compensated by something else. However previous investigations found a finite value for det . So something seems inconsistent.
To connect to previous work, let us perform the angular-momentum decomposition as in Section 4.3. This lets us write the functional determinant as with ∆ s in Eq. (4.22). For s = 0 there is one mode, the dilatation mode, which has zero eigenvalue, so R 0 = 0. For s = 1 there are four zero modes corresponding to translations. For s ≥ 2 all the eigenvalues are positive. Removing the zero modes from the numerator, Ref. [11] found R 0 ≈ −1 and R 1 ≈ 0.041. First, we look at s ≥ 2. Here there are no zero modes, so there are no issues with rescaling the operators for these values of s. That is, with R s≥2 in Eq. (4.33). For s = 0 and s = 1 there are zero modes. Since zero modes are still zero modes if the operator is rescaled, we know the explicit form of these modes. They are in Eq. (4.25) with n = 1 and s = 0 or s = 1 φ 10 = R R 2 − r 2 (R 2 + r 2 ) 2 , φ 11 = R 2 r 2(R 2 + r 2 ) 2 (A.4) It is easy to check that (∆ 0 − 3λφ 2 b )φ 10 = 0 and (∆ 1 − 3λφ 2 b )φ 11 = 0. Because of the zero modes, we need to compute Note that the zero modes become modes with eigenvalue of the shifted operators, so the shifted determinant will be proportional to as desired. For s = 1, the zero modes are translations and the Jacobian is finite (Eq. (2.15)). Thus we expect R 1 to be finite too. To compute R 1 we can first try the Gelfand-Yaglom method as in Section 4.4. The Gelfand-Yaglom method requires us to find a solution to ∆ s − 3λφ 2 b + φ 1ε = 0 (A.6) that scales like the free solution, φ 1 = r near r = 0 and r = ∞. Unfortunately, this does not work. At finite , φ 1 is oscillatory at large r while the free-theory solution φ 1 = r is not. Thus the two cannot approach each other and the Gelfand-Yaglom method does not seem to give a sensible answer.
To understand the failure of the Gelfand-Yaglom method we note that adding the term as in Eq. (A.5) is equivalent to adding a mass term 1 2 φ 2 to the potential. One would think that a small mass would be a small change in the theory, but it actually has a dramatic effect: it removes all bounce solutions to the equations of motion. Thus the limit → 0 is not smooth. One can deal with small masses using the constrained instanton approach described in Section 7, however, there is a simpler way to compute R 1 .
Since R 1 is supposed to be finite, we can rescale the operator as for s ≥ 2: We can now evaluate these determinants in the basis of s = 1 modes, φ n1 given in Eq. In this derivation we have used a property of Legendre polynomials, that drr 3 φ n1 φ m1 = R 2 6n(n + 3) δ nm (A.12) As → 0 the first term in Eq. (A.11) always dominants unless λ n = 1. Thus the n = 1 mode contributes εR 2 24 to the product and we can set = 0 for the other modes. We then find This factor of 1 24 = 0.041 matches the result from numerical calculations in [11]. Now we try the same approach for the s = 0 modes. In this case, the relevant mode is φ 10 in Eq. (A.8), corresponding to dilatations. Attempting the same calculation as for s = 1 we find R 0 = R 0 drr 3 φ 2 10 = ∞ · R 0 (A.14) Thus we conclude that the ratio of determinants in Eq. (A.1) is zero, even after the zero mode is removed. This result, while in disagreement with finite numerical extractions in [11,19], is expected from our treatment using rescaled operators in Section 3.
In conclusion, we find the infinite Jacobian for dilatations found in previous work to be compensated by an infinitity of the determinant after zero modes are removed. We therefore find no inconsistency in the functional determinant calculated with or without rescaling the operators.

B Divergent Graphs in Fermi Gauge
In Fermi gauge, the Lagrangian is given in Eq. (5.30). We treat all the mass terms as interactions. The Feynman rules are To proceed, we want to compute the divergent contributions with Feynman diagrams in d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected.
Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction The sum, after taking µ → µ √ To proceed, we want to compute the divergent contributions with Feynman diagrams in d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected.
Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction The sum, after taking µ → µ √ To proceed, we want to compute the divergent contributions with Feynman diagrams in d dimensions and subtract the corresponding contribution from the 4D result to sum over s. Unfortunately, performing the subtractions in Fermi gauge is difficult. In Fermi gauge, due to the gφ b A µ ∂ µ G term in Eq. (5.29) there is a Feynman rule picking up the momentum of virtual Goldstones. This extra loop momentum generates new UV divergences and makes the diagrams difficult. This is explained in more detail in Appendix B where we compute all the divergent parts (but not the finite parts). These divergences exactly correpond to those in Eq. (5.71) as expected.
Fortuntately, we can compute the regularized contribution in any gauge. Indeed we have checked numerically that our result for the finite s functional determinant is identical in R ξ gauge and Fermi gauges. In R ξ gauge, with ξ = 1 with Lagrangin in Eq. (5.27), the Feynman rules are (5.72) In addition, in R ξ gauge the ghosts do not decouple and have an interaction ← q = −g 2 φ 2 b (q), (5.73)