Bubble Nucleation to All Orders

This paper extends classical results by Langer and Kramers and combines them with modern methods from high-temperature field theory. Assuming Langevin dynamics, the end-product is an all-orders description of bubble-nucleation at high temperatures. Specifically, it's shown that equilibrium and non-equilibrium effects factorize to all orders, and that the nucleation rate splits into a statistical and a dynamical prefactor. The derivation clarifies, and incorporates, higher-order corrections from zero-modes. The rate is also shown to be real to all orders in perturbation theory. The methods are applied to several models. As such, Feynman rules are given; the relevant power-counting is introduced; RG invariance is shown; the connection with the effective action is discussed, and an explicit construction of propagators in an inhomogeneous background is given. The formalism applies to both phase and Sphaleron transitions. While mainly focused on field theory, the methods are applicable to finite-dimensional systems. Finally, as this paper assumes an effective Langevin description, all results only hold within this framework.


Introduction
The observation of gravitational waves is a game-changer, granting us new eyes to gaze at the cosmos [ -]. There is now a realistic chance of glimpsing phase transitions that occurred a few nanoseconds after the Big Bang: if the transition is first order. Contrary to continuous phase transitions, a first-order transition can leave tracks in the form of a stochastic gravitational-wave background. This is because these transitions proceed through nucleating bubbles [ , , -]. The interior of these bubbles is permeated by a true vacuum; the outside situated in a metastable vacuum. Once these bubbles expand they can generate gravitational waves through collisions and turbulence in the primordial plasma [ -]. Furthermore, these bubbles might facilitate Electroweak Baryogenesis-thus potentially explaining the observed Baryon asymmetry [ -].
Although, it is unclear if a first-order phase transition occurred in our cosmological history, not to mention whether gravitational-waves from such a transition can be detected by upcoming experiments [ -].
To answer these questions both model building [ -] and theoretical calculations [ -] are required. Theoretical calculations, in particular, need to be improved to reduce uncertainties in the description of rate-of-production, and successive growth, of nucleating bubbles [ , -].
In the context of perturbative calculations, one issue is that calculations converge slowly at high temperatures, so higher-order corrections can be large. As such, calculating these higher-order corrections provide much-needed cross-checks. This is especially relevant for the nucleation rate, where radiative corrections can be enhanced for large bubbles.
To that end, this paper shows how to calculate the bubble-nucleation rate to higher orders.
The main result of this paper is that the nucleation rate factorizes as Here A dyn accounts for non-equilibrium effects and A stat captures equilibrium effects. These factors are known as the dynamical and the statistical prefactor respectively, and they can be calculated order-by-order in perturbation theory. Equation . is derived in Section , and applied to a real-scalar model in Section . Crucially all derivations assume a classical Hamiltonian system. See [ -, -] for the details and applicability of this effective description. It is important to stress that formulas similar to Equation . have been suggested before [ , , -]. However, it is unclear from past literature how higher-order corrections should, even in principle, be included.
In this paper, the nucleation rate is derived from first principles, and the factorization in Equation . is proved to all orders. Moreover, Feynman rules for calculating the statistical and dynamical prefactors are derived. It is also shown how to treat zero-modes at higher orders. In addition, explicit calculations show that the rate is renormalization-scale invariant to two-loops.

. Non-equilibrium dynamics
Thermal escape is a non-equilibrium process. As such, damping and related effects are important [ , , ]. However, these effects are hidden in Equation . . In addition, the saddle-point approximation is expected to break down for small damping [ ]. This motivates going beyond leading-order results [ , ].

. Thermal escape and tunneling
To see how Equation . should be modified, it is informative to discuss the physics behind thermal escape.
The formulas for quantum tunneling and thermal escape look similar However, the physics is quite different. For tunneling, the rate is¹ The effective action is evaluated on a solution of δS eff [φ b ] = 0. In essence the tunneling rate comes from a saddle-point approximation around the bounce solution. The functional determinant (one-loop contribution to S eff ) arises from fluctuations around the bounce, and higher-order terms in S eff consist of vacuum diagrams in the bounce background [ , ]. Crucially the effective action is imaginary due to a negative eigenvalue of δ 2 S B [φ b ].
The thermal-escape formula is similar to the tunneling decay rate. Indeed, except for the prefactor, the determinant and exponent could be the first terms of e −S eff /T , where S eff is the effective action in three dimensions. This has led to the conjecture [ , ] The leading-order dynamical prefactor is A dyn = λ 2π , and λ is normally taken as |κ|, where κ is the only negative eigenvalue of Yet thermal escape is a classical process.
To illustrate the physics it is useful to consider a classical particle moving in three dimensions. We assume that the potential of this particle has one metastable minimum at x = x S , and a lower-energy minimum at x = x TV . Furthermore, let us assume that these two minima are separated by a barrier. Lowest at x = x B . Without loss of generality we take x B to lie in the x 1 direction.
If this was a tunneling process, the particle would start at the metastable minimum and go through the barrier with some probability. The most likely path of escape, the bounce, results in a negative eigenvalue, which gives an imaginary energy. The rate is identified But the picture is different for thermal escape, where the particle travels over the barrier. That is, the probability for the particle to have an energy E is proportional to e −E/T . This means that the particle is most likely to escape where the barrier is lowest (V B = V (x)| x= x B ), with probability Γ ∼ e −V B /T . Although, the particle can also escape close-by x B by moving slightly in the x 2 or x 3 directions.² A better approximation is Γ ∼ x 1 =x B d x 2 d x 3 e −V (x)/T . Crucially only variations in orthogonal directions to x B are included in this integral. The rate is real because there is no integration in the concave x 1 direction.
What about non-equilibrium physics? When the particle jumps over the barrier, it first needs to get sufficient energy to escape, and then move to the other side. If the negative curvature (κ) in the concave direction is large, the probability to get moving over the edge is naturally big. Likewise, large damping (η) slows down the particle. This is the physical reason for the leading-order dynamical prefactor: A dyn = |κ| 2πη . Intuitively the Boltzmann factor describes the probability to jump up to, and the dynamical factor describes the flow from, the barrier. To leading order this flow is in the concave direction, while higher-order corrections make the flow veer down the barrier from orthogonal directions.
Note that it does not matter that V is concave in the x 1 direction because we should only integrate in orthogonal directions to x 1 when calculating the rate. Thus no imaginary terms can appear in the saddle-point approximation.³ These considerations motivate the main result of this paper where the effective action omits contributions from negative eigenvalues. This means that the effective action is real to all orders. Equation . is derived in the next section for an overdamped system; an equivalent derivation for general damping is given in Section . The complete nucleation rate, in field theory, is given in Section . .

The nucleation rate
In this section we derive a complete formula for the nucleation rate. To that end, consider a system with n degrees of freedom, indexed by i. Translating the results to field theory is straightforward, and is done in Section .
For simplicity this section assumes an overdamped system. While not generic, an overdamped system makes the physics transparent. The steps are similar for the generaldamping case, differing only by longer intermediate formulas, which are given in Section . As mentioned, the calculation of the nucleation rate is classical. This follows because field theory is, to first approximation, classical at high temperatures [ -]. This is indicated by the Boltzmann factor At high temperatures T ω, the occupation number is large: n(ω) ≈ T /ω 1. Equivalently, the system is classical if the thermal-wavelength is small: However, although nucleation occurs at low energies, large energy modes with ω ∼ T still contribute through loops. These high-energy modes can be integrated out, and give effective couplings; effective damping parameters; and effective thermal noise [ ].
In this section the damping is taken as a free parameter. See [ -] for calculations, and discussions, of the damping for specific models.

. Classical nucleation theory
This section reviews classical nucleation theory. We stress that the results in this first subsection are known [ , , , , ].
The Einstein summation convention is used for the indices i, j, k, l, where the vertical position of these indices is unimportant.
Consider the Fokker-Planck equation for a system with large damping (also known as Here V (x) is a generic potential, ∂ i is the derivative with respect to x i , and η is the damping.⁴ The function P(x, t; x 0 , t 0 ) is the probability to find the particle at a position x, given that it was at x 0 at time t 0 . Then, if P is initially localized to a metastable state we can ask what is the probability-flow from this state.
Determining P is in general hard, yet to calculate the rate it suffices to determine P close to the barrier. The rate can then be approximated by the probability flow, defined by J i , across this barrier [ , , ]: The calculations become tractable if one introduces artificial sources and sinks so that P is static [ , , , ]. In that case the problem is reduced to solving the static Fokker-Planck equation. After the equation is solved, and we have our current, the sinks and sources can be removed. This is known as the flux-over-population method. Although, it should be mentioned that this method fails for underdamped systems [ ].
It is assumed that all components [x i ] have the same damping η.
To use the flux-over-population method, consider first the metastable state, taken to be at x = x S . The (static) equilibrium solution is P 0 = e −V (x)/T , which is the normal Boltzmann distribution and follows from Equation . . Further, if the potential is of the close to the metastable point, the denominator in  Equation . can be calculated with a saddle-point approximation.
Next, the idea is to determine P close to the barrier-generally a saddle point. This requires two approximations. First, we assume that the surface integral in Equation . is dominated by the saddle point. Second, we assume that the surface can be approximated as a surface normal to the saddle point.
To find P close to the saddle point it is useful to follow Langer and make the ansatz P = ζ( x)e −V (x)/T . The function ζ( x) describes deviations from thermal equilibrium and behaves as where x TV denotes the true (lower-energy) vacuum state. Using this ansatz in Equation . one finds To solve this equation, assume that the saddle point is at x i B . At this point the potential is to first approximation quadratic: This equation can be solved with the ansatz ζ( x) = ζ(u), where u is a linear combination of x i . Namely, u = U i x i . Then where U i must be an eigenvector of ω i j for consistency. That is For simplicity we choose the normalization U 2 = 1. The solution to Equation .
is [ , ] The boundary conditions for ζ(u) force κ < 0, and so U i is the eigenvector corresponding to a negative eigenvalue. As expected, deviations from equilibrium only depend on the negative-eigenvalue direction.
Given ζ(u), the current is where P 0 = e −V (x)/T . To find the rate we integrate the current on a surface normal to u = 0. This gives the rate⁵ where to leading order wheredet means that all negative and zero eigenvalues are excluded, see Section . for the details. Crucially u = 0 excludes the negative-eigenvalue direction, and the determinant is manifestly real. Normalizing with the denominator in Equation . gives the known result [ , , ] We refer to this as the leading-order rate from now on. Consider now corrections to Equation . , of which there are two kinds. First, corrections to the dynamical prefactor come from including more terms for V in Equation  . . Second, the statistical prefactor receives corrections from additional terms in V via Equation .
. To leading order we can write the rate as: As the next two sections show, the statistical prefactor consists of normal (exponentiated) vacuum diagrams; the dynamical prefactor, in field language, consists of operator insertions.
Henceforth we leave the metastable normalization implicit. This normalization factor is straightforward to calculate via the effective potential [ , ].
. Corrections to the statistical prefactor Consider the statistical prefactor. As mentioned, the negative eigenvalue does not cause any issues because everything is calculated at u = 0.
To proceed, consider the integral in Equation . : Integrating over k enforces u = 0, and we leave this integration for last. Barring some subtleties with the negative eigenvalue, the integrals give Technically U i ω −1 i j U j < 0, so the integral over k formally diverges. Being more careful one can let k → −ik and do a Wick-rotation. This gives We identify Z[U] as a generating function: correlators can be calculated via We see that Z[U] is a generating function with a non-vanishing current. In terms of Feynman diagrams this means that in addition to usual diagrams, there are now diagrams with external-current insertions.
To see the form of these external-current insertions, consider the potential The Feynman rules for this potential are given in Figure . = Figure : Feynman rules for the statistical prefactor.
In addition, if there are 2n external-current insertions, the diagram picks up a factor (−1) n κ n (2n − 1)!!. This last rule comes from integrating over k. Also, diagrams with an odd-number of current insertions vanish since ∞ −∞ d kk 2n+1 e −k 2 A = 0. At next-to-leading order (NLO) the diagrams are given in Figure . Using the Feynman rules in Figure   the second 1 Note that the propagator can be written where V i a are eigenvectors of ω i j with eigenvalue λ a . After using this form of ω −1 i j , the last two diagrams remove κ from the first.
In summary, the statistical prefactor is whereS eff denotes all vacuum diagrams omitting zero and negative eigenmodes.

. Corrections to the dynamical prefactor
Consider now the dynamical prefactor. Again, we assume a potential Including higher-order terms in V does two things. First, the Boltzmann factor changes. This is encoded in the statistical prefactor. Second, the rate of flow across the saddle-point changes. This flow is controlled by deviations from thermal equilibrium, which is encoded in the dynamical prefactor.
Because Equation . describes the deviation from equilibrium, we need to find ζ to higher orders. Explicitly, we need to solve T ∂ 2 i ζ − ∂ i V ∂ i ζ = 0 with the inclusion of λ 3 and λ 4 terms. The leading-order solution, denoted henceforth as ζ 0 , is given in Equation . .
The boundary conditions force ζ( x) to vanish as u → ∞, which is enforced directly by the ansatz It is assumed that ζ 1 = O λ 3 , λ 4 . The equation for ζ 1 can be found by combing Equations . and . : ( . ) We use the notation x i x j x k x l . All homogenous solutions of Equation . grow exponentially with u, and so do not satisfy the boundary conditions. This leaves the particular solution.
To find the particular solution, note that if we make a polynomial ansatz [in x i ], neither of the terms on the left-hand side of Equation . increase the polynomial's degree. Whence the ansatz go through the equations. It is also useful to work in the eigenbasis of ω : This basis is useful because V i a is orthogonal to U i and since ω i j x j reduces to a sum of eigenvalues. The idea is to make a polynomial ansatz in u and v a , where a denote all eigenvalues except for the negative one.
For simplicity, consider first the two-dimensional case (x i ) = (x 1 , x 2 ) and ignore the x 4 term for the moment. That is, we assume the potential There are two eigenvectors: Since U x 2 is a degree 2 polynomial, the relevant ansatz is After using this ansatz in Equation . , and collecting terms, one finds Note that most of the terms in ζ 1 are irrelevant. This is because the rate is given by For the above example only the C 3 term contributes since if the derivative (∂ i ) hits anything else, that term is proportional to u or U i V i . Both which vanish. This holds in general: we only have to care about terms containing one factor of u.
The n-dimensional case is similar. Including the quartic term in Equation . , the ansatz is a polynomial of degree 3: Following the discussion above, only the B 1 , C a , and D ab terms contribute to the rate. One finds The complete solution for ζ 1 is given in Appendix B, and the result for potentials with x 5 and x 6 terms is also given in Appendix B. Note from Equation . that ζ 1 can be written in terms of Green's functions of the form As before a excludes the negative eigenvalue.
Before showing how ζ 1 affects the rate, note that ζ 1 can considered as a set of operator insertions. To see this, pick one term from Equation . , for example D ab uv a v b . When calculating the probability current, the factor of u disappears, and the term is proportional This operator can also receive loop corrections like in normal field theory.
Additional corrections to ζ can be found by using ζ = ζ 0 + ζ 1 ζ 0 + ζ 2 ζ 0 + . . .. Where it is assumed that ζ n = O (λ n ). In general, ζ n is a polynomial of degree 4n − 1, and is completely determined by ζ n−1 . Namely, for n > 0 Naively one would think that ζ 2 is sub-leading compared to ζ 1 . However, from Equation . we see that the c a term does not contribute directly to the rate. Instead this term contributes via the one-loop tadpole shown in Figure   This tadpole insertion scales as λ 2 3 , and similar terms exist in ζ 2 . As such, all λ 2 3 terms in ζ 2 are relevant.
We here only note that these terms are of the form The explicit expressions are given in Equation C. In summary, the dynamical prefactor consists of operator insertions, which multiply the full statistical prefactor.

. Comparison with known results
Sections . and . derived the nucleation-rate for systems with n degrees of freedom. There are however pre-existing results for n = 1. Take a potential The NLO rate is [ , ]: where as before the normalization from the metastable state is left implicit. This formula includes all corrections to the dynamical and statistical prefactors. It is encouraging that the statistical prefactor in Equation . together with the dynamical prefactor in Equation .

Generic damping coe cients
Consider now the generic-damping case. There are a few new features as compared to the overdamped system. First, the velocity, or rather the conjugate momenta, contribute to the Boltzmann factor. Second, the damping is not solely a multiplicative factor as in Equation .
. Third, the negative-eigenvalue direction contributes to the Boltzmann integral. As before, the starting point is the Fokker-Planck equation [ , , ]: here π i are conjugate momenta associated with x i , and η is the damping.
To derive the rate, we start by assuming Next, we make the ansatz P For consistency so U i x is also an eigenvector of ω i j . Putting everything together one finds [ , ] It will be shown below that λ needs to be positive, so perforce κ < 0 and we must choose λ = 1 2 η 2 − 4κ − η . Next, using Equations . and . gives [ , ] γζ The solution satisfying the boundary conditions is Given ζ(u), it is straightforward to calculate the rate. Explicitly [ , ] where E = 1 2 π 2 i + V (x). Doing the Gaussian integrals one finds up to a normalization As before there are two types of contributions from higher orders. First, corrections to ζ, and second, higher-order corrections to E in the integral u=0 d n x d n ve −E/T . Before calculating these corrections, note that the rate is still real. To see this, denote the negative-eigenvalue direction by x u ≡ U i x x i . The problematic term in the energy is If the rate is calculated at We see that the Gaussian integral in Equation . converges, and that the rate is real.

. Physical interpretation and applicability
The rate is calculated at u = 0. Let us now see what this means. For clarity we use That is, x u and π u are the position and conjugate momentum in the negative-eigenvalue direction. Denote the energy in this direction by E u ≡ 1 2 π 2 u + 1 2 κx 2 u . Consider now the motion of a particle starting at the saddle point. To do so, recall the static Fokker-Planck equation: Namely, deviations from thermal equilibrium follow Equation . , which is known as the adjoint Fokker-Planck equation. Now, the Fokker-Planck equation corresponds to a particle obeyinġ the last term represents thermal noise and vanishes on average. Likewise, the backward Fokker-Planck equation follows froṁ Imagine now releasing a particle close to x u = π u = 0. If the particle evolves according to Equation .
, an unstable solution is [ ] On this solution the energy is As expected d d t E = −ηv 2 < 0. Next, consider instead a particle following Equation .
. The unstable solution is now with associated energy This is precisely the solution implied by u = 0: π u = |κ| λ x u . Note that the energy in Equation . suppresses fluctuations away from the saddle point via the Boltzmann factor. This is why the Boltzmann-integral is localized at x u = 0 for η → ∞.
Equation . also highlights a problem for small damping. Namely, for η = 0 the energy in Equation . identically vanishes. Effectively the negative-eigenmode turns into a zero-mode. While the rate is finite in the η → 0 limit to leading order, this does not generalize to higher orders.
There are other indications that the formalism breaks down for small damping [ ], and one should be cautious when using existing formulas [ , ] for an underdamped system.

. The statistical prefactor
Consider now corrections to the statistical prefactor. As in Section . we start with the integral where We see that Z[U x , U π ] is a generating function, and the rule is The integral converges because with γ given in Equation .
. As before, consider the potential The Feynman rules for this potential are given in Figure . In addition to these rules, if there are 2n external-current insertions, the diagram picks up a factor (−1) n γ −n (2n−1)!!; diagrams with an odd number of current insertions vanish.
Similarly, for integrals over the conjugate momenta, the external current gives a factor Feynman rules for general damping.
x , and the 2n rule is the same. The derivation of the dynamical prefactor is given in Appendix A.

Field theory
As yet we have not considered field theory, however, the field-theory limit is straightforward [ ]. From Section we known that these fields live in three dimensions. Furthermore, as discussed in Section . , high-temperature effects can be captured by using effective parameters in this theory [ , ]. Yet for the purposes of this paper we do not consider an explicit effective theory [ , ], and three-dimensional parameters are taken as free. We will however absorb all temperature-dependence by rescaling couplings and fields. In addition, we omit details associated with the assumed, effective, Langevin description. Uncertainties associated with this description are important [ -, -], but lie beyond the scope of this paper.
Consider a real-scalar model with three-dimensional action The barrier position x B is replaced with the bounce φ b (| x|) [ , ]; the action evaluated on the bounce is analogous to the barrier height: In Section we assumed that In field theory this potential is replaced with (to fourth order) where δ y is shorthand for δ δφ( y) . We see that the dictionary is where ∂ 2 ≡ ∇ 2 . Furthermore, the solution to the Fokker-Planck equation is [ ] where for an overdamped system ζ(u) is given by . Feynman rules for the statistical prefactor Consider the potential Expanding around the bounce gives⁶ Linear terms are omitted since δS LO Feynman rules for a real-scalar theory.
Using the result in Section . , we obtain the Feynman rules given in Figure . The propagator is defined by ; all zero eigenvalues are projected out as explained in Section . . Also, this propagator depends on the bounce, and needs to be found numerically. See Section . for the details.
The Feynman rules work as usual, with the addition that if there are 2n external-current insertions, the diagram picks up an additional factor (−1) n κ n (2n−1)!!, and diagrams with an odd number of current insertions vanish.
As an example, consider the diagram in Figure   Using the Feynman rules given in Figure we find The dynamical prefactor is calculated via operator insertions as explained in Section . . These operators are given in Figure . The first vertex acts as a tadpole, the second as a 2-point insertion, and the third is a number. Note that the second vertex can be reduced to a sum of propagators once it contracts with other lines. For example, because V a (x) are eigenfunctions of the propagator, we have Then the leading correction involving the 2-point insertion is given in Figure    After using the Feynman rules we find where the shifted propagator satisfies .

Zero-modes
We now turn to zero-modes. In field theory there are three zero-modes associated with the bounce [ , ]. These zero-modes are readily identified: We can remove these zero-modes by using collective coordinates [ , ] φ where a excludes zero-modes. Essentially the zero-modes turn into a volume factor. In particular, by using collective coordinates the zero-modes never appear in propagators.
Using collective coordinates also gives a Jacobean. To leading order this Jacobean is (S B ) 3/2 . However, as emphasized in [ ], there can be additional terms. To see the effect of these terms, one can extend the result of [ ] to a three-dimensional field theory. The result is The first term in J Z M gives the familiar prefactor (S B ) 3/2 , which cancels when calculating the statistical prefactor [ , ]. The other terms are operator insertions that need to be taken into account order-by-order. Note that these zero-mode operators talk with the dynamical prefactor.

. Connection with the e ective action
It is sometimes useful to view the rate as an effective action. Take for example the statistical prefactor This effective action is the sum of all vacuum diagrams in the bounce background. In getting to the right-hand side, we have assumed that φ b is a generic background field. To obtain the rate one should fix φ b so that δS eff [φ b ] = 0.
We could of course calculate the integral in Equation . numerically without any reference to the effective action. However, the effective action can be useful. On the practical level the effective-action approach removes tadpole diagrams. For example, omitting external-current insertions, the statistical prefactor in Equation . is given by the diagrams in Figure   The third diagram is not part of the effective action. Instead, corrections to the leadingorder bounce incorporate all tadpoles. Explicitly, if the tree-level action is S LO and the The second term in A 3 can be symmetrized over x, y and z.

one-loop correction is
Including φ NLO in the effective action gives Diagrammatically S NLO [φ LO ] comes from the double-bubble and sunset diagrams, and is the dumbbell diagram. The actual expression for φ NLO is given by This follows from the definition φ NLO : Note that both the propagator, and implicitly the effective action, are defined with zeromodes projected out.
. Renormalization-scale invariance of the rate Calculations should be independent of the renormalization scale. To show that the rate is scale-independent is however non-trivial. Note that renormalization-scale dependence first shows up first at two loops in three dimensions. Moreover, the two-loop beta function is the complete beta function, and only mass parameters depend on the scale in three dimensions. For the real-scalar model considered in Equation . , the beta function is Before considering the rate, it is instructive to first calculate the effective potential. That is, we treat φ b → φ as a constant background field.
The relevant diagram is given in Figure   Integrations are left implicit.
Using normal Feynman rules, the contribution to the effective potential is where the last step only works if the background field is constant. The propagator in the background field is defined by The solution, in d = 3 − 2ε dimensions, is [ ] here K 1/2−ε (x) is the modified Bessel function. Note that for ε = 0 the propagator is Note that all ε poles, and associated scale-dependence, come from the |x| → 0 region. While contributions from the |x| → ∞ region are finite. The idea is to introduce a radial cut-off R, and for |x| < R one should use the d = 3 − 2ε expression for the propagator, while for |x| ≥ R one can directly take the ε → 0 limit [ ]. Using properties of the Bessel function one finds |x|<R d 3 x∆(x) 3 = 4ε log(µ 3 R) + (2 + 4γ)ε + 1 64π 2 ε .
( . ) Note that the log R term cancels when including the region |x| ≥ R [ ].
Consider the scale-dependent piece: which cancels with the running of the tree-level mass: With the effective-potential case done, most steps carry over to the nucleation rate. First, note that the scale dependence of the leading-order bounce action is⁹ d d log µ 3 φ b terms cancel since the action is evaluated on-shell.

The sunset diagram gives
We omit terms linear in φ b as they correspond to tadpoles. The region of interest is x ≈ 0, and in this limit the propagator satisfies¹⁰ Thus to leading order this propagator is the same as for the effective-potential case [ ]. We find The integral with respect to y is the same as in Equation .
, and gives We see that the effective action is renormalization-scale invariant. Next, consider the dynamical prefactor. Up to a factor of 2π, the renormalization-scale dependence is and the relevant -loop diagram is given in Figure   Using the Feynman rules from Section . this diagram is¹¹ We find The same result follows after explicitly calculating the propagator as shown in Section . .
The minus sign comes from the U n rule. Note also the sign of κ as compared to |κ|.
The dynamical prefactor is manifestly renormalization-scale invariant once we multiply the above diagram with the leading-order prefactor |κ| 2π . Although this section only treated the real-scalar model, everything carries through once vector-bosons are included. The same methods can also be used to check that the Sphaleron rate is renormalization-scale invariant.

. Radiative barriers and the dynamical prefactor
The previous section dealt with scale-dependence arising from two-loop diagrams. Consider now instead a model where a heavy field is integrated out. For example a vector boson. The procedure of integrating out particles, and consistently calculating the rate, has been studied in [ , , , , ].
This scenario is relevant as even if a barrier is absent at tree-level, it can be generated from loops [ , ]. As an example, consider a SU(2) gauge theory with a scalar in the fundamental representation; the leading-order potential is The next order includes a scale-dependent contribution [ ] Normally the scale dependence in V NLO cancels with the beta function of the mass parameter: This is manifestly the case for the statistical prefactor, but it is more complicated for the dynamical prefactor. Indeed, naively the negative eigenvalue satisfies which implies Let us check that this procedure is correct by using the methods in Section . . Treating V NLO as a perturbation, the scale-dependent piece in ζ 1 is (see Equation . ) 256π 2 e 4 u log µ 3 .
( . ) As shown above, the scale-dependent piece in κ is d d log µ 3 κ = 51 256π 2 e 4 . This implies Using the contribution from δζ 1 , we see that the dynamical prefactor is renormalizationscale invariant.

. Factorization of the rate
In the previous sections we have seen that the rate factorizes as The statistical prefactor is Here S eff [φ b ] omits all contributions from negative and zero eigenmodes. The zero modes also generate operator insertions; these are included in J Z M as discussed in Section . . For concrete calculations, the Feynman rules given in Section . can be used. The rate should be normalized with e − d 3 x V eff (φ S ) , where V eff (φ S ) denotes the effective potential in the metastable state.¹² The dynamical prefactor is defined as The virtue of these definitions is that the renormalization-scale dependence for A dyn and A stat are simple; both A stat and A dyn are manifestly real to all orders; A dyn and A stat can be calculated independently; there is a clear physical interpretation with A stat as the Boltzmann suppression, and A dyn as the rate of probability-flow across the saddle-point.
To leading order A dyn can be interpreted as the growth-rate of the bubble.

Considerations for concrete calculations
This section shows how to calculate propagators in an inhomogeneous bounce background; what diagrams appear at two loops; and how to isolate divergences in dimen-This last point is important since A stat is not finite without the normalization. In practice V e (φ S ) should be included diagram-by-diagram as discussed in Section . . sional regularization.

. Propagators
The propagator satisfies the equation Because the bounce is spherically symmetric, it is useful to expand the propagator in spherical harmonics [ ] with the short-hand r ≡ | x|, r ≡ | y| and cos α = x· y Given the bounce, Equation .
can be solved numerically for each l.
There are two caveats. First, for practical reasons one can not solve Equation . for infinitely many l. Also, the propagator diverges in the x → y limit.
Both of these problems have a common solution [ , ]. The idea is to solve Equation . for l ∈ (1, L) where L is typically between and . The remaining part of the sum can then be done analytically by using the WKB approximation.
There are two cases. First, assume that cos α < 1. Solving Equation .
(see Appendix D) for large l one finds to leading order The δ l term can be evaluated numerically for given L.
For the second case we assume cos α = 1, and find 1 4π Note that the leading δ L term is finite in the |x| → | y| limit. Higher-order terms are given in Appendix D. Zero modes only contribute when l = 1. To remove the zero-modes, start with Equation Now use the representation where with eigenvalue λ a . From Section . we know that the normalized zero-modes are [ , ] This leads us to consider the equation where we can remove the zero eigenvalues by defining This modified propagator, also known as a generalized Green's function, satisfies Again, this equation can be solved by expanding ∆ (x, y) in spherical harmonics; as mentioned, the ∂ µ φ b (x)∂ µ φ b ( y) term only contributes to l = 1.
In summary: Use Equations . and . to find the propagator for l = 0, . . . L. Choose L according to the required precision.
If α < δ, use Equation . to do the sum from l = L + 1 to ∞.
When calculating integrals of the form d 3 x d 3 y∆(x, y), use the result from point to analytically calculate the contribution from the x → y region.
All ε poles, with d = 3 − 2ε, come from the x → y region.
Integrate over the remaining phase-space numerically.
Finally, external-current insertions can be omitted after replacing . The statistical prefactor Consider now possible diagrams appearing at the two-loop level. As before we use the real-scalar model defined in Section with a large damping. Diagrams contributing to the statistical prefactor can be divided into three classes: vacuum diagrams, external-current insertions, and zero-mode insertions. Note that at NLO-equivalent to two loops for equilibrium observables-there are contributions from lower-loop diagrams with operator-insertions.
The Feynman rules are the same as in Section . , and κ denotes the negative eigenvalue, while U(x) is the corresponding eigenvector normalized as d 3 x U(x) 2 = 1. All propagators are defined with zero-modes removed.

. . Vacuum diagrams
There are three vacuum diagrams at two loops, which are given in Figure . The bubble and sunset diagrams are only finite after normalizing with the effective potential evaluated in the metastable state. As such, for actual calculations one should always add the corresponding diagram for the effective potential with a relative minus sign. For example, the double-bubble diagram in Figure   The divergence comes from the |x| → ∞ region. For numerical stability one should include the effective-potential term directly: where ∆ S (x, x) is the propagator for a constant φ = φ S . Since all diagrams exponentiate we can directly include them in the effective action. This gives an extra minus sign for each diagram. To estimate the size of each diagram we use the power-counting The  and the dumbbell diagram gives Finally, the sunset diagram scales as .

. External-current insertions
Here diagrams are labeled by the order they appear, left to right and top to bottom. All diagrams scale as λ m . The diagrams are given in Figure . D In addition, the notation A(x, y) ≡ (λφ b (x) − g)(λφ b ( y) − g) is used. In the above, D 5 represents the square of diagram , D 6 the square of diagram , and D 7 the product of diagrams and . Note that all diagrams are finite in the | x| → ∞ limit.

. . Zero-mode insertions
As noted in Section . , when removing zero modes we have to include the operators given in Equation .
. We denote these operators with black circles to differentiate them from other operators. In addition, external-current insertions are omitted for brevity.
The relevant diagrams are shown in Figure . In terms of propagators these diagrams Figure : Zero-mode diagrams contributing to the statistical prefactor. give Note that D 1 ∼ D 2 ∼ λ m , while D 3 ∼ λ 2 m 2 . This means that only D 1 and D 2 are relevant at NLO.

. Dynamical prefactor
The dynamical prefactor is given in Section . . There are two types of contributions. Those coming from a quartic vertex and those coming from two cubic vertices.
The quartic contribution is given in Equation . and gives The double-cubic contribution has one part coming from Equation . , and another from C. . The contribution from Equation . results in two diagrams, which are given  in Figure . The first diagram is

y)∆( y, y)A(x, y),
and the second is .

. Zero-mode insertions
As mentioned, zero-mode operators also contribute to the dynamical prefactor. There are two one-loop diagrams (ignoring external-current insertions), these are given in Figure  . The first diagram gives The second diagram scales as λ 2 m 2 , and so is not relevant at this order.

. Summary
In this paper we have shown that the nucleation rate factorizes into an equilibrium and a non-equilibrium contribution-a statistical and a dynamical prefactor. This factorization holds to all orders assuming an effective Langevin description. The formula derived in this paper lays the groundwork for calculating the nucleation rate beyond leading order. This is critical for getting a hold on uncertainties, and in particular, for estimating when the formalism breaks down.
To illustrate the calculations, this paper derived special Feynman rules. In addition, explicit calculations confirm that both the dynamical and the statistical prefactors are renormalization-scale invariant to two-loops.

. Future work
This paper does not show that the dynamical prefactor is gauge invariant. Even though the statistical prefactor is gauge invariant [ , , , ], it is still important to verify that this holds for the full rate.
Another avenue is to apply the formalism of this paper to Sphaleron transitions. This is of great importance for Baryogenesis scenarios [ ]. In particular, the literature lacks a robust calculation of the Sphaleron rate with a radiative barrier. Moreover, current calculations of the functional determinant are limited to the Standard Model [ , ].
Lastly, this paper is based on classical field theory. This is motivated at high temperatures: integrating out high-energy modes results in a classical field theory with effective damping, couplings, and thermal noise. Explicit calculations of the damping are however sparse [ -, ]. There is also the question whether a Langevin description is applicable at higher orders. Answering these questions is essential for robust predictions.

. Conclusion
It is the hope that the results of this paper will aid the theoretical understanding of phase transitions, and thereby reduce uncertainties for gravitational-wave predictions. Yet there are several open questions, and many problems remain. Solving these problems requires applying existing methods to new models; pushing high-temperature field theory to higher orders; and doing vital cross-checks with lattice computations.

A The dynamical prefactor for general damping
Consider the potential The goal is to solve Equation . and find ζ. That is, we need to solve The leading-order solution, refereed henceforth as ζ 0 , is given in Equation .
. Following Section . , make the ansatz ζ = ζ 0 + ζ 1 ζ 0 (u). Plugging this ansatz into Equation A. one finds As in Section . the idea is to use a polynomial ansatz for ζ 1 , however, this ansatz is more involved than in the overdamped case. To see how things work, focus on the U π x 3 term.
The idea is to use the eigenbasis ω: x . The ansatz is Solving Equation A. , and keeping only terms that contribute to the rate, gives Recall that λ = 1 2 η 2 − 4κ − η . Other terms are given by similar formulas.

A. Comparison with regular perturbation theory
To check the results, consider ω i j → ω i j + εω i j 1 where ε is a small number. In this case we already know the answer because the correction to κ follows from usual perturbation theory: where U i x ω i j = κU j x . As before, make the ansatz ζ = ζ 0 + εζ 1 ζ 0 where ζ 0 is given in Equation . . After some algebra one finds Following the previous section, the ansatz is ) where V i a are eigenvectors of ω i j with positive eigenvalues λ a . After collecting terms we get We use the notation 〈U V 〉 ≡ U i ω i j 1 V j . Just as in the overdamped case, the rate involves integrating over a surface normal to u. Since eigenvectors are orthogonal, all A a and B a terms drop out.
Adding ζ 1 changes the rate by With ζ 1 taken care of, this still leaves contributions from the Boltzmann factor. This contribution is Adding these two contributions together is equivalent to taking |κ| → |κ| − U This result agrees with regular perturbation theory.

B The dynamical prefactor for dimension operators
Recall Equation . T where the ansatz was For brevity we use v a bc ≡ v a v b v c and V abc ≡ V a V b V c ; all repeated indices are summed over.

Consider now a potential
In the language of SM-EFT the above potential would correspond to including a dimension operator in the tree-level Lagrangian, see for example [ ].
This new potential necessitates including the terms One finds

Remaining terms are given by Equation B. with the modification
The above factors, once included in the rate, can be written in terms of propagators. See Section for the details.