Computing the gauge-invariant bubble nucleation rate in finite temperature effective field theory

A gauge-invariant framework for computing bubble nucleation rates at finite temperature in the presence of radiative barriers was presented and advocated for model-building and phenomenological studies in an accompanying article arXiv:2112.05472. Here, we detail this computation using the Abelian Higgs Model as an illustrative example. Subsequently, we recast this approach in the dimensionally-reduced high-temperature effective field theory for nucleation. This allows for including several higher order thermal resummations and furthermore delineate clearly the approach's limits of validity. This approach provides for robust perturbative treatments of bubble nucleation during possible first-order cosmic phase transitions, with implications for electroweak baryogenesis and production of a stochastic gravitational wave background. Furthermore, it yields a sound comparison between results of perturbative and non-perturbative computations.


Introduction
A possible first-order electroweak phase transition (EWPT) in the early universe introduces two tantalizing possibilities: electroweak baryogenesis [2][3][4][5] and the production of a primordial gravitational wave background [6][7][8][9]. A first-order transition that proceeds through bubble nucleation is incompatible with the Standard Model (SM) of particle physics due to the large Higgs boson mass. Instead, electroweak symmetry breaking (EWSB) in the SM universe proceeds via a smooth crossover [10][11][12][13][14][15]. However, a modified SM scalar sectorwith potentially observable deviations of the SM Higgs boson properties or the discovery of new particles in future collider experiments -could imply that a first-order EWPT occurred at temperatures of T ∼ 100 GeV [16]. This possibility has invigorated numerous studies exploring prospects for a first-order electroweak phase transition in beyond Standard Model (BSM) scenarios (see  and also a more extensive list of references in [16]) with interest in particle physics model-building, phenomenology, and experiment. This research program heavily relies on calculating thermodynamic properties of the transition, such as its critical temperature and strength described by the released latent heat, as well as quantities that characterize the dynamics, such as the bubble nucleation rate. In fact, the existence of a first-order transition for a given choice of BSM parameters does not guarantee that it will have occurred. The nucleation rate, Γ, must be sufficiently large to ensure a transition out of the false vacuum. Obtaining theoretically robust computations of Γ is vital for testing the viability of a first-order EWPT.
Achieving a gauge-invariant, perturbative computation of the nucleation rate at finite temperature in the presence of a radiatively-induced barrier between the two phases of the transition has been a major challenge. While the presence of enhanced contributions from infrared bosonic modes casts doubt on the reliability of perturbative treatments, it is possible to alleviate other shortcomings of perturbation theory [46] through e.g. a consistent use of the renormalization group [47] and inclusion of higher-order effects [48]. Moreover, when the potential barrier between the two phases in the transition appears at tree-level, it is relatively straightforward to perform a gauge-invariant computation of Γ (see e.g. recent [46,47]). For tunneling at T = 0, gauge-invariant approaches for radiative barriers [49][50][51][52] are applied to the analysis of vacuum stability in the SM; see also Refs. [53,54]. Implementing gauge invariance at finite T when the barrier arises from gauge sector loops has persisted as a thornier problem.
A recent accompanying publication [1] has presented a practical, gauge-invariant framework for performing a perturbative computation of Γ(T > 0) for a radiative barrier. This framework was presented in the context of the leading-order high-temperature expansion familiar to the model-building and phenomenology communities, and generalizes the methods of Ref. [50] under clearly defined limits of validity: 1 (i) The nucleation temperature, T n , must be sufficiently close to the critical temperature, T c , such that the leading, thermal corrections to the effective thermal mass cancel against tree-level terms. Thus, one may adopt a power counting in the gauge coupling wherein the leading-order contributions to the radiative barrier are gauge-independent.
(ii) The gradient expansion is only well-behaved for the first few orders of perturbation theory. To go to higher orders requires other methods.
For (i) we assume thick-wall bubbles that require temperature ranges with large enough supercooling away from T c , but close enough to adopt a power counting for perturbative expansion that allows for a gauge invariant computation. We emphasize that gauge invariance is a necessary, but not necessarily a sufficient condition for a self-consistent and reliable perturbative description. Our framework has direct applications to strong transitions with thick-wall bubbles, as we describe in later sections. Importantly, the problem of radiatively-induced barriers -the subject of the present study -remains eminently relevant to the thermal history of EWSB in extended scalar sectors. For a recent general discussion, see Ref. [16]. The real triplet extension of the SM [55][56][57][58] provides a concrete example, wherein EWSB may occur in either one-or two-steps. In the one-step case as well as the first transition of the two-step scenario, the barrier is entirely radiatively generated. Both transitions may become first-order for suitable choices of the parameters. Moreover, for the one step case, thermal loop effects associated with the Higgs portal coupling are decisive, as the transition would be a smooth crossover in the absence of this interaction (e.g. the SM case). In both cases, the presence of the first-order transition provides the needed preconditions for successful electroweak baryogenesis.
From a theoretical standpoint, it is important to place the framework of [1] in a context allowing for a systematic treatment using the general effective field theory approach [59]. By complementing [1], we employ the powerful technique of high-temperature dimensional reduction [60,61], working in a three dimensional Euclidean effective field theory (3d EFT) as introduced in Refs. [62][63][64]. Then, extending [59], we construct an effective description for bubble nucleation within the 3d EFT with gauge fields. This allows us to show the attainability of a consistent description for nucleation and link the high-temperature nucleation to classical nucleation theory. Previous work has employed 3d EFT methods to analyse bubble nucleation in [65][66][67][68], yet these works have not demonstrated gauge independence of the nucleation rate. Although our setup is generic for gauge field theories with a radiative barrier, we focus on the Abelian Higgs model by following [49,50,69] to simplify the calculations and reasoning. A computation in this illustrative model captures relevant features for a gauge-independent calculation, and the methods described below can be generalized to realistic models of cosmic phase transitions, in particular with non-Abelian gauge fields. See, for example, the recent work of Ref. [67] on the SU(2) + Higgs theory that exhibits a radiative barrier.
To summarize the key features of our analysis, we determine the bubble nucleation rate (cf. Sec. 2) in powers of gauge coupling g (weak coupling constant) and the numerical coefficients a 0,1 . Two terms in the exponent are computable in the derivative expansion of the effective action at leading (LO) and next-to-leading orders (NLO), respectively. Higher order effects are inaccessible in the derivative expansion. 2 The prefactor has mass dimension four, and we leave its scaling in terms of the weak expansion parameter unspecified. After taking a logarithm, the rate reads The central argument here is that the first two terms in powers of g describe the leading behavior of ln Γ, which is relevant for obtaining the inverse duration of the phase transition. Thus, ln Γ is a key input parameter for determining the gravitational wave spectrum [8,9,46]. Our task is to compute the exponent terms (a 0,1 ) and demonstrate their gauge invariance.
In the framework of [1], this task is achieved at NLO in coupling expansion, but at leading order in the high-temperature expansion. Here, we extend this computation to resum all relevant NLO contributions from the hard thermal scale by utilizing dimensional reduction. Consequently, we show the cancellation of the renormalization scale related to thermal resummations, as well as among the two exponent terms, signalling a consistent perturbative treatment [47,64]. To be able to capture higher order thermal contributions, the EFT approach provides physical insight and consistency. It enables one to determine the order at which various contributions arise and to assess the limits of validity of the different expansions employed herein: the high-temperature, gradient, and coupling expansion. Earlier work [69] concluded that perturbative computations at high temperature introduce an artificial gauge dependence of the nucleation rate. This feature is related to the breakdown of the gradient expansion of the effective action in the symmetric phase, and seems to prohibit 2 Technically, many higher-order terms (at three-and higher-loops) related to contributions of parametrically heavier fields than the nucleating field are still accessible in the derivative expansion [67]. They appear at higher orders than the leading behavior of the prefactor A. a gauge-invariant treatment. In what follows, we show that this conclusion does not apply to the exponent terms in Eq. (1.1), provided that assumptions (i) and (ii) are satisfied.
This article focuses on a purely perturbative determination of the nucleation rate. However, due to infrared enhancement of the bosonic sector at high temperature, a fully comprehensive study of the phase transition thermodynamics requires non-perturbative simulations on the lattice [70]. For equilibrium properties of the transition, such lattice analyses appear e.g. in [10,58,65,[70][71][72][73][74][75]. For bubble nucleation, non-perturbative lattice studies are limited [65,66,68] with applications in [73]. The framework of [1] and this work at hand, provide a sound basis for comparing results of perturbative and of non-perturbative computations.
This article is composed as follows. Section 2 introduces the model and reviews the wellknown zero-temperature computation of [49,50] that utilizes the derivative expansion in the computation of the nucleation rate in perturbation theory. In addition, we discuss the extension to high temperature along the lines of the accompanying article [1]. Section 3 reformulates the same problem in 3d EFT language, using a general framework [59]. This formulation allows us to systematically organize thermal resummations and to better monitor intermediate gauge dependence pertinent to different scales in the perturbative computation. In Section 4, we summarize our computation and discuss its implications for other nucleation rate computations in the literature. Appendix A proves the gauge invariance of the dimensional reduction step by deriving high-temperature matching relations between the fundamental four-dimensional theory and three-dimensional effective theory. Appendix B explicates computational details in 3d EFT perturbation theory.

Nucleation, radiative barriers, and the derivative expansion
We work with a simple gauge field theory toy model as Refs. [49,50,69] to compute the bubble nucleation rate between different vacua. The Abelian Higgs model 3 can be defined by the Lagrangian density with B µ a U(1) gauge field (with gauge coupling g) and Φ a complex scalar. The covariant derivative for the complex Higgs reads and the hypercharge for the complex scalar Y φ = 1. Since our goal is to compute the Euclidean action, we already define the Lagrangian in the Euclidean, rather than Minkowski space. We expand the complex field in terms of real fields
where φ is a scalar background field and H and χ are propagating degrees of freedom. We apply general R ξ -gauge fixing [76]: with gauge fixing functional F (Φ, Φ * ). The latter implies the corresponding ghost Lagrangian where c, (c) are (anti)ghost fields. A priori bothφ and φ are unrelated but eventually identifiedφ = φ to eliminate the mixing between vector boson and Goldstone mode and remove mixed propagators. Gauge-fixing choices are comprehensively discussed in [76,77]. Before focusing on the central part of this article, the formulation of the thermal tunneling rate, we review its zero-temperature analog.

Zero temperature calculation: a review
Tunneling between two vacua in quantum field theory was first properly examined by Coleman and Callan [78,79] in analogy with the calculation of quantum mechanical tunneling rates. This analogy was re-examined by computing the tunneling decay rate directly from the Minkowski path integral by using a physical definition of the tunneling probability [80,81].
In the so-called bounce formalism of Coleman and Callan, the tunneling rate can be expressed as a path integral that is dominated by contributions from a bounce solution φ B (x), a field configuration that extremizes the Euclidean effective action of the theory, S eff . In perturbation theory we can expand where φ b extremizes the leading-order effective action and ∆φ = φ B − φ b collects higher order contributions that correspond to quantum corrections to the shape of the bounce. Formally, the tunneling rate (per unit volume) in a zero-temperature four-dimensional QFT can then be calculated through [82] where S 0 (φ) is the leading-order action, δ δφ S 0 (φ b ) = 0, and φ f.v. is the false vacuum from which the tunneling proceeds. The primed determinant, det ′ , excludes zero-modes related to translational-invariance. The O( ) term above encodes higher order corrections. However, to connect our calculation to the analysis of [50], we focus on the form given in terms of the effective action Γ = Im where V is the four-dimensional volume. Using the effective action to formulate the rate in this way has not been proven to be valid at all orders. However, it correctly reproduces the leading-order terms we discuss below (cf. related discussion of Sec. 6.1 in [81]). We continue to calculate the tunneling rate by using a derivative expansion of the effective action, where the ellipsis indicate terms involving additional powers of ∂ µ φ. For now, we assume that such an expansion is appropriate, 4 and that we can perturbatively expand the effective potential V eff (φ) and the kinetic field renormalization term Z(φ), where V (φ) denotes the leading-order potential and the + · · · higher order terms in the couplings. The leading-order bounce solution that extremizes leading S eff is then a radially symmetric solution of the following equation of motion and boundary conditions Next, we focus on a radiatively generated barrier, as considered in Refs. [49,50], and review their computation to establish the procedure for our finite-T generalization. Since the tree-level Lagrangian (2.1) contains no barrier, it admits no tunneling. A barrier arises via quantum corrections: integrating out the vector boson yields a barrier between two minima of the resulting effective action [83]. To examine this possibility, we first consider the background-field dependent squared masses of the fields: Here, B is the gauge boson and H the "Higgs" field, which are massive at the broken minimum; χ, (c) corresponds to the Goldstone (ghost) field which receives a gauge-dependent contribution to its mass in R ξ -gauge. Now consider the tree-level potential, V 0 (φ), and the one-loop correction, V 1 (φ): Here, the one-loop master function J 4 (x) is with master integral J d given in Eq. (B.30) and To induce a radiatively generated barrier, loop and tree-level effects need to be of similar size. As a consequence, the loop expansion breaks down, even if one may still retain a well-defined coupling expansion. Coleman and Weinberg [83] demonstrated that, indeed, a perturbative expansion in g remains valid by counting λ ∼ g 4 . In our case, where we are interested in tunneling between vacua, we also have a quadratic term (φ 2 ) with positive coefficient µ 2 . Counting λ ∼ g 4 requires µ 2 ∼ g 4 σ 2 for tunneling to be possible [50], where σ is a characteristic value of φ b (x = 0) ∼ σ, and in our case also the vacuum-expectation-value in the stable phase. The relevant power counting is To find the leading-order potential at O(g 4 ), we expand V 1 (φ) in powers of g and find . (2.20) This potential contains the tree-level potential (first two terms), the vector boson contribution to the one-loop potential (third term) and has two different minima that are separated by a barrier. Furthermore, it is gauge-independent. By inspecting the Goldstone and ghost terms in Eq. (2.16) we infer that the gauge-dependent terms are of O(g 6 ) according to the scaling relations in Eq. (2.19). One subtlety merits mentioning here. The cancellation of the nominally O(g 4 ) Goldstone and ghost contributions does not occur near the false vacuum. Nevertheless, the impact on Γ is suppressed by the small field values in the regions of non-cancellation and will be beyond the eventual accuracy goal of the computation, which is ln Γ = −a 0 g −4 − a 1 g −2 + ln A, in analogy to Eq. (1.2).
Since the leading-order effective potential (2.20) is gauge invariant, solutions to the corresponding leading-order equations of motion where ≡ ∂ µ ∂ µ , will also be gauge invariant. We expand the effective potential and the wavefunction renormalization in the coupling g . . , (2.24) which in turn are used for computing the effective action in the derivative expansion (2.8).
Expressions for the next-to-leading (NLO) corrections Z g 2 and V eff g 6 , where the latter includes both one-and two-loop contributions, have been computed in Refs. [50,84]. Here we merely need their counting in terms of g and not their explicit expressions. For illustrative purposes, we also include subdominant terms, Z g 4 and V eff g 8 additional to leading order (LO) and NLO terms presented in Ref. [50]. As we will see momentarily, these subdominant terms contribute at an order where the derivative expansion of the effective action breaks down.
Using a similar notation as in [49,50], the nucleation rate reads where the prefactor A (with mass dimension four) must be computed using the fluctuation determinants in Eq. (2.6) and results in a higher order effect than the LO and NLO exponent terms, regarding the logarithm of the rate. The exponent terms read We can determine the expected sizes of B 0,1 using the power-counting together with the characteristic size of the critical bubble, with radius R ∼ m −1 H ∼ g −2 σ −1 , which is determined by the leading-order potential. As a result (2.28) To understand the breakdown of derivative expansion, we can imagine calculating the next corrections in the same manner as above, i.e.
where the next order denoted by the ellipsis arises at O(g 0 ), and in fact an infinite number of higher order derivative terms in derivative expansion contribute at the same order -the derivative expansion does not converge. The breakdown occurs because loops containing propagating scalar fields first appear at this order. In general, a well-defined derivative expansion requires a separation of scales, leading to an expansion in powers of P/M with formal P ∼ ∂ and M being a mass scale M ≫ P . Taking a Fourier transform of the leadingorder bounce (2.22), we have which describes the characteristic "nucleation scale", P ∼ g 2 σ. The diagrams contributing to B 0 and B 1 have propagating vector bosons, so that M = m B . The resulting expansion parameter is P 2 /m 2 B ∼ g 2 . The loops at higher order, on the other hand, also include scalars, with the corresponding expansion parameter P 2 /m 2 H ∼ 1. Thus, a well-defined derivative expansion is applicable only when integrating out degrees of freedom that are heavy with respect to the nucleation scale, P (or, in another words the fluctuations with wavelengths much shorter than the nucleation length scale). In the present case, the heavy degrees of freedom are the physical vector bosons. Integrating them out yields the barrier, and the characteristic scale over which the bounce solution changes is small compared to the vector boson mass [49,50].
Based on these general observations, the derivative expansion nominally applies for B 0,1 . Since V eff g 4 is manifestly gauge-invariant, so is B 0 . The gauge invariance of B 1 remains to be demonstrated. Moreover, since the heavy scale m B ∼ gφ, and since the integrals in Eqs. (2.26) and (2.27) include regions of vanishingly small φ, one rightly worries whether contributions of order P 2 /m 2 B are, in fact, finite. This is manifested in the logarithmic φ-dependence of the NLO wavefunction renormalization, Z g 2 ∼ ln(φ/Λ) raising concerns about the finiteness of B 1 . However, upon closer examination, one finds that the contribution to B 1 is finite. To this end, consider the asymptotic behavior of the bounce φ b (r) at large r (small φ). In this region ∂V eff wherein the bounce equation and its solution read and c is an undetermined constant. Now, we divide the region of integration to two domains: (i) r ≤ R and (ii) r > R, with R being larger than the characteristic size of the bounce. Applying the asymptotic solution (2.33) to region (ii) the contribution from the possibly problematic terms ln(φ) (∂ µ φ b ) 2 to B 1 is proportional to which is finite.
In the following, we address the question of gauge invariance which can be understood with the help of the Nielsen identities [76,85]. For a derivation of the Nielsen identities in a derivative expansion, see the original result [50], and [69] for an amendment relevant for higher orders. The variation of the effective action with the gauge parameter can be expressed as which is the Nielsen identity with x ≡ d D x. The corresponding Nielsen functional is implied by the gauge transformation variations δ g Φ = igΦ, δ g Φ * = −igΦ * and the variation of the R ξ gauge fixing function (2.3) The above functional also admits a derivative expansion which, combined with the expansion of the effective action (2.23), results in the Nielsen identities for the effective potential and field renormalization factor In perturbation theory, the Nielsen coefficients C,D, andD are expanded as [50] C = C g 2 + C g 4 + . . . , (2.41) with an explicit derivation for C g 2 in [50,84]. In fact, neither the functionals D andD nor the correction C g 4 enter the test of gauge dependence of B 0 and B 1 due to their high scaling in g. Specifically, the derivatives of Z and V eff with respect to φ are all at least O(g 2 ), so that the first, third, and fourth terms on the right hand side of (2.40) are all at least O(g 4 ) whereas the leading non-trivial gauge-dependence of Z enters at O(g 2 ). Thus, for the quantities relevant to B 0,1 , the Nielsen identities become which have been verified explicitly in [50].
Using these formulas at each order in g, we can test the gauge dependence of B 0 and B 1 . The leading-order term B 0 is immediately gauge invariant as no gauge fixing parameter enters such that ξ ∂ ∂ξ B 0 = 0. The NLO term B 1 is also gauge invariant, which can be established with the help of the equations of motion [50,84]: We remark that in Eq. (2.21), we expanded in terms of m 2 G /m 2 c which is technically only allowed for φ that are not too small. But as we have seen in Eq. (2.34), even though formally Z diverges when φ → 0, the dominant contributions to the rate are still finite and gauge invariant. This completes our review of the zero-temperature computation.

Finite temperature calculation: the conventional approach
In the following, we generalize the zero-temperature analysis of [50] to finite temperature. We first investigate a conventional calculation of the effective action, S eff (T ), in accordance with earlier literature [69]. We, however, depart from the argument [69] that gauge invariance of finite temperature nucleation rate at NLO cannot be established due to a breakdown of the derivative expansion. This section complements the analysis in the accompanying article [1] by supplementing several technical details.
Nucleation rate at finite temperature was first discussed in [86,87]. In our analysis below, we define the thermal nucleation rate by where κ describes dynamical, real-time non-equilibrium phenomena related to thermal fluctuations from the meta-stable to the stable minimum, and Σ is a statistical part, that describes equilibrium, time-independent properties of nucleation. Above, we have only anticipated a form of Σ in terms of prefactor, and exponential, in which the leading orders are enhanced and computable from the effective action. Sec. 3.1 defines the statistical part more carefully within the EFT approach of [59]. We start with a brief review of concepts in the imaginary-time formalism of thermal field theory, required to compute the statistical part Σ. One may formulate the latter as a 4d Euclidean field theory with a compactified Euclidean time-dimension. This is the Matsubara formalism, in which the mode expansion for fields entails an integration over the threemomentum modes and a sum over Matsubara modes [88] viz.
where we denote P ≡ (ω n , p) for Euclidean four-momenta and the temperature-dependent bosonic Matsubara frequency is ω n = 2πnT . We use dimensional regularization in D = d + 1 = 4 − 2ǫ dimensions and the MS-scheme with renormalization scale Λ similar to zero temperature. It has been conventional to consider the high-T expansion of this formalism, wherein the bosonic one-loop function reads A high-temperature expansion in µ/T can be formally defined by assuming the scaling µ ∼ gT . The mass µ is a so-called soft mass scale of the theory, opposed to the hard scale of non-zero Matsubara modes that have parametrically larger mass ∼ πT . We now implement a commonly-followed approach for computing the T > 0 effective action: (i) integrate out non-zero Matsubara modes and (ii) implement the "daisy resummation" of zero-mode masses at leading order for the scalar field and temporal or time-like gauge field [89]. As a preview of subsequent sections, we note that this approach can be justified by 3d EFT methods in Sec. 3, and reproduces the correct leading behavior therein. However, within this conventional approach higher order corrections are not straightforwardly accessible, albeit their numerical importance due to slower convergence of the perturbative expansion at high-T [47,64]. We insert the leading contributions of Eq. (2.49) into the oneloop potential (2.16), discard field-independent terms, and obtain Here, the linear-in-T term requires resummation as implemented below. Before doing so, let us consider the sum of the tree-level, T = 0 one-loop, and V 1 (φ, T ) as given in (2.50) The presence of the linear-in-T term introduces a barrier between the symmetric and broken phases, implying the existence of a first-order transition at critical temperature T c . For the moment, however, we focus on the resulting term that is quadratic in φ, whose T -dependence governs the onset of spontaneous symmetry breaking: and where the + · · · denote the remaining non-quadratic terms. By examining the behavior of µ 2 eff , we observe that for very large temperatures, it will be positive and large, as the positive-definite ∝ T 2 terms dominate the negative µ 2 -term. For such large temperatures the µ 2 eff φ 2 term dominates the effective potential which implies that only the symmetric phase is attainable at large temperatures. For very small temperatures we instead have that µ 2 eff < 0: the symmetric phase is unstable. And in between these temperatures there is a temperature T 0 where µ 2 eff = 0. In the absence of the other terms in the potential, T 0 would define the critical temperature for a second order transition from the symmetric to the broken phase. Thus, one must have that for T near T 0 , µ 2 ≈ − 4λ + 3g 2 T 2 12 . For temperatures different from, but close to T 0 , the cancellation between the µ 2 and g 2 T 2 components of µ 2 eff is not exact, but there will exist a region for which µ 2 eff /T 2 ≪ O(4λ + 3g 2 ). It is natural to parametrize the degree of µ 2 eff suppression with additional powers of g: Here, we will assume that T c lies within a temperature range for which N = 1 applies. Under this assumption, one may define a consistent power counting in g. Near the phase transition, all terms in the potential should be roughly of the same order of magnitude, a feature one may implement by taking It is possible that for a range of temperatures near T 0 the suppression of µ 2 eff is stronger: i.e. the effective mass is ultrasoft, cf. Eq. (3.12). In this case perturbation theory breaks down and the system is non-perturbative. Henceforth, we focus on scenarios for which the relations in (2.54) hold for a range of temperatures near T 0 that includes T c . Generally, this counting λ ∼ g 3 -required for a radiatively generated barrier at finite temperature -was introduced by Arnold and Espinosa [89], and further studied in Refs. [69,90]. For a discussion on the scaling of µ 2 eff ∼ g 3 T 2 also cf. [91]. Including the O(T φ 3 ) term in V eff changes the nature of the transition from second to first order. Before analyzing the implications of (2.54), we recall that in this context a consistent treatment of thermal loops requires to perform a "daisy resummation". In practice, the latter amounts to replacing the field-dependent masses in (2.50) by the corresponding thermal masses. Gauge invariance implies that the masses of the spatial gauge bosons, and ghost fields, remain unchanged. One then has are the squares of the temporal gauge boson Debye mass, Higgs boson, and Goldstone boson, respectively and where The temperature-dependent part of m 2 L corresponds to the leading Debye screening: hard thermal excitations of the plasma screen the temporal gauge field which renders it massive even in the unbroken phase.
The contributions from the ghost and Goldstone bosons to the linear-in-T term introduce an explicit ξ-dependence. However, the power counting (2.54) implies that this gauge dependence appears at higher order in g than the gauge-independent contribution from the transverse and temporal gauge fields: The resulting LO effective potential then becomes Proceeding with the tunneling rate calculation, the leading-order bounce is solved from Here ≡ ∂ i ∂ i is the 3d Laplacian operator. As in the T = 0 case, one may in principle solve for corrections to the bounce solution, ∆φ, by including the higher order terms in (2.65).
In general, these corrections enter the ln Γ beyond the two leading orders of interest here.
Exceptions may occur, such as in the thin-wall regime [59]; see the end of Sec. 3.4 for a detailed discussion. For the action in the derivative expansion, we need expansions of the effective potential and field renormalization factor Expressions for V eff g 4 and Z g are given in the accompanying article [1] and we present Z g also below in Eq. (2.71). We compute both terms in detail within the 3d EFT in Appendix B. V eff g 4 (Z g ) contain contributions from transverse and longitudinal gauge bosons and ghosts at two-loop (one-loop) level. In addition, V eff g 4 includes the leading difference of Goldstone and ghost terms at one-loop, Eq. (2.62). Both V eff g 9/2 and Z g 3/2 arise from Higgs loops at one-loop level, and are not required. In analogy to the notation of [49,50], we write (the statistical part of) the nucleation rate as where β ≡ 1/T . As at zero-temperature, we do not compute the prefactor A. The gdependence of B 0,1 follows from the power counting of Eq. (2.54) and the characteristic bubble size R ∼ µ −1 eff ∼ g −3/2 T −1 . The latter is determined by the leading-order potential: The next-order exponent term arises at O(g 0 ), and in analogy to zero temperature (cf. discussion around Eq. (2.30)), calculating this order in the derivative expansion would require an infinite amount of terms. Thus the derivative expansion breaks down also at finite temperature, though again the leading two terms are attainable (cf. also Ref. [59]). Before discussing gauge invariance, let us first ensure that the contribution from B 1 is finite. First note that Here, the second term corresponds to the contribution from the gauge field temporal mode. We compute these contributions within the 3d EFT approach in Appendix B (cf. Eq. (B.66)), and have here converted to parameters of the fundamental theory at leading order. Importantly, at high temperature the leading correction to the field renormalization Z does not explicitly depend on the gauge fixing parameter unlike at zero temperature.
Here, we contrast to an existing computation in Ref. [69]: Our expression in Eq. (2.71) agrees with the result for the broken phase in Eq. (B.15) therein. However, the authors argue that one must use their result in Eq. (6.1) outside the broken phase (which corresponds to our Eq. (B.65)), which would introduce additional gauge dependence. In our power counting, this issue does not arise and we discuss this in more detail after Eq. (3.25) The presence of the 1/φ term in (2.71) renders Z g (φ) more singular than the T = 0, Z g 2 (φ) wavefunction correction. One may thus wonder whether the 1 2 Z g (∂ µ φ b ) 2 contribution to B 1 is finite. We proceed as before using the asymptotic behavior of the bounce solution: and study the possibly problematic contribution by dividing the radial integration into two regions (see discussion around Eq. (2.34)): which is finite. We now demonstrate the gauge invariance of B 0,1 . As before, B 0 is trivially gauge invariant as a gauge fixing parameter is absent at this order. The Nielsen identities also hold at finite temperature [69], and their coefficients expand as where the coefficients C, D,D are derived in Appendix B at leading order. However, like at zero temperature, we merely need the leading-order expression for C and the identities which have been verified to hold in Ref. [1], and which we further validate explicitly within the 3d EFT approach in Appendix B. In particular, both sides of Eq. (2.78) vanish identically since at leading order the correction to Z is ξ-independent (cf. Eq. (2.71)) and C is φindependent (cf. Eq. (B.77)). At high temperature, the proof that ∂ ∂ξ B 1 = 0 compares to its zero-temperature analog in Eq. (2.45): The third line moved C g outside the integrand due to its φ-independence. Note that the sequence of steps looks manifestly different from what occurs in (2.45) since the kinetic contribution to B 1 is explicitly gauge invariant as implied by the φ-independence of C g . Both cases rely on the vanishing of the surface integral, which follows from the asymptotic behavior of the bounce solution in Eq (2.73). This completes the proof of gauge invariance of the exponent B 1 . In fact, with the help of the Nielsen identities we reached this conclusion by merely knowing the powercounting for the next-to-leading order V eff g 4 . Its explicit expression is, however, relevant for numerical applications. To this end, the next section will employ the technique of high-temperature dimensional reduction and use 3d effective field theory. This allows us to implement transparently the required thermal resummations and to organize the two-loop level computation systematically.
Finally, also here one can question whether the derivative expansion is well-behaved since Eq. (2.62) is an expansion in powers of m 2 G /m 2 c , which diverges in the limit φ → 0. But the situation is similar to zero temperature. Inconsistencies introduced this way in the derivative expansion enter only at higher orders. For a detailed discussion of non-local terms that are ignored in the limit φ → 0, see [59].

High temperature effective theory
We now place the foregoing discussion in the context of the dimensionally-reduced, hightemperature effective field theory (3d EFT). Thus, we can (A): define the thermal nucleation rate, by matching to classical nucleation theory, as was done in [59], (B): systematically incorporate thermal resummations, and access important higher order corrections. This allows us to eliminate otherwise problematic renormalization-scale dependence in the perturbative expansion [46,47], (C): assess the limits of validity of the Γ-computation, while remaining within the context of perturbation theory.
We begin with some general remarks. The characteristics of the nucleating bubbles are set by the long wavelength behavior of the theory. In this context, the infrared physics of the high-temperature plasma is related to the static modes of the theory. Effectively these modes are described by a three-dimensional theory where heavy modes in the temporal direction are integrated out -this is the idea of dimensional reduction. Recalling that the thermal plasma gives rise to a rigorous scale hierarchy, at every distinct scale it is possible to construct such a dimensionally reduced EFT. Initially established in the context of non-Abelian gauge theories [60,61], the formalism [62,63,92] is used widely in hot QCD (cf. [93] for a review), and is also becoming increasingly popular for electroweak theories beyond the Standard Model [46-48, 57, 58, 71-75, 94-109]. For a recent tutorial that applies dimensional reduction to the singlet scalar field theory, see [110]. The Lagrangian density of the high-temperature 3d EFT for the Abelian Higgs model [111][112][113][114], has a structure similar to Eq. (2.1). Its fields and couplings are replaced by 3d quantities and are denoted by subscript "3" i being the field strength tensor for the spatial U(1) gauge field B 3,i (with gauge coupling g 3 ) and where Φ 3 is the 3d complex scalar. The 3d EFT character of the vector boson is labelled by the first subscript while the second index is a spatial Lorentz index i = 1, 2, 3. The covariant derivative for the complex Higgs reads In addition, there is a thermal remnant of the gauge field temporal component: the temporal scalar B 0 . Its 0-subscript merely labels its origin from the temporal component of the 4d gauge boson.
The corresponding temporal sector of the effective Lagrangian reads with m 2 D the Debye mass and κ 3 the B 0 self-interaction coupling. For an Abelian gauge field, its temporal remnant B 0 is a singlet that merely couples to the Higgs via the portal coupling h 3 . Couplings κ 3 and h 3 originate in analogy to the Debye mass: screening of the hard scale induces thermal corrections to interactions and not just the mass. In the conventional approach only the mass is resummed while dimensional reduction also accounts for the resummation of interactions. Couplings to the spatial gauge field B 3,i are absent because the Abelian gauge field does not self-interact in the fundamental 4d theory.
We use generic R ξ -gauge defined in analogy by the gauge fixing Lagrangian Eq.
The 3d gauge fixing parameter is denoted as ξ 3 . Here Φ 3 is the scalar field andφ 3 an external, generic background field, that is in general separate from the field expectation value φ 3 ≡ Φ 3 . In the end, we identifyφ 3 = φ 3 which eliminates the mixing between the Goldstone mode and the gauge field. The relevant Faddeev-Popov ghost Lagrangian [69] reads after varying the gauge-fixing function F 3 (ψ i ) of Eq. (3.3) with respect to its fields where the partial derivative acts in direction of the arrow and θ parameterizes infinitesimal gauge transformations. By construction [69] L 3d,GF and L 3d,FP have the same relative sign compared to L 3d .
We omit higher dimensional operators from the 3d EFT which is justified at next-toleading order dimensional reduction [62]. While the dimensional reduction for the Abelian Higgs model is known [64], we independently reproduce results and explicitly ensure gauge invariance of the reduction step along the way. We relegate details of the reduction to Appendix A and details of computations within 3d perturbation theory to Appendix B.
One obtains the parameters of the 3d theory by requiring equality between 3d and 4d Green's functions. Implementing this requirement leads to a set of matching relations that where we denote Green's functions of the parent 4d theory by G and field renormalization factors by Z. By virtue of matching, only the hard scale contributes to G and Z above via non-zero Matsubara modes, as detailed in Appendix A. Expanding the above expressions at the desired order in couplings gives rise to the 3d EFT parameters that depend on the original model parameters and temperature. We emphasize that at NLO these 3d parameters do not depend on the gauge fixing parameter associated to dimensional reduction. The perturbative computation introduces an individual gauge dependence on both G's and Z's which cancels exactly in the above matching relations, as we detail in Appendix A. However, this intermediate gauge dependence in the construction of the 3d EFT should not be confused with the gauge fixing within the 3d EFT perturbation theory (cf. Appendix B).

Overview of effective field theory setup for bubble nucleation
This section formulates the computation of the bubble nucleation rate within the 3d EFT along the lines of [59]. Thermal fluctuations from the meta-stable to the stable minimum correspond to the nucleation of bubbles. Though this is a non-equilibrium process, the nucleation rate can be factored into (cf. Eq. (2.46)) a dynamical part, κ, that captures the non-equilibrium phenomena and a statistical part, Σ, calculable within the 3d EFT ) corresponds to the nucleation scale effective action obtained by integrating out heavier excitations than the nucleating degree of freedom. The critical bubble, φ b , is a stationary configuration of S nucl . Despite similarities in the formulas between zero temperature vacuum decay in Eq. (2.6) and the nucleation rate in Eq. (3.11), there is a different physical, effectively classical, picture to the nucleation [59], wherein the nucleation rate formula follows from Langer's nucleation theory [115][116][117]. Evaluating the determinant around the critical bubble corresponds to the contribution of differently shaped nucleating classical bubbles. A primed determinant means that the translational-invariance zero-modes are excluded. The statistical part is normalized such that κ is the exponential growth rate of the nucleating bubbles [116]. The remainder of this article focuses on the statistical part Σ without further discussing the dynamical prefactor κ.
The spatial extent of the critical bubble sets a length scale, the nucleation scale Λ −1 nucl . The nucleation scale can be identified with the mass of the nucleating d.o.f., Λ nucl ∼ µ eff (cf. Eq. (2.54)), away from the thin-wall limit. Given that there is a hierarchy between the nucleation scale and higher intermediate (soft) 3d scale, one can use the effective field theory framework to organize the calculation.
The relevant scales within the Abelian Higgs model are 5 πT thermal scale Step 1 ≫ gT intermediate scale Step where both thermal and intermediate scales are higher scales to be integrated out for the nucleation scale effective description. 6 The computation separates into two parts as presented in Tabs. 1 and 2.
Step 1 in Tab. 1 is the usual dimensional reduction, as described at the beginning of this section and in Appendix A. This step takes care of the highest energy scale of the theory, the thermal scale (πT ), when constructing the 3d EFT with Lagrangian in Eq. (3.1) describing the length scales of (gT ) −1 .
Step 2 in Tab. 2 matches the intermediate (soft) scale to the nucleation scale, by integrating out heavier degrees of freedom than the nucleating field. Thus, we can create a local description for the length scale, µ −1 nucl , of nucleating bubbles where µ nucl ∼ (g 3/2 T / √ π). Here, we included a factor of π coming from the one-loop order of the heavy contributions to λ ∼ g 3 /π. The effects from the higher scales enter the local nucleation scale effective action, S nucl , that can be used to find an approximation for the critical bubble.
The nucleation scale and lower scales enter the fluctuation determinants in the prefactor. We will concern ourselves with the two leading orders of the exponent, and leave the determination of the next order, i.e. the nucleation scale contributions, as future work. The critical bubble background has strong effects on these scales: The gradient expansion diverges as variations of the bubble background cannot be treated as small external momenta within the loop integrals, and the fluctuations around the critical bubble background contain an unstable negative eigenmode and zero-modes (cf. fluctuation determinant in Eq. (3.11)), which need careful treatment. Contrary to the higher scale contributions, the effects of these scales 5 The assumed scaling for the nucleation scale does not necessarily apply in the full parameter space, and our discussion is limited to those regions where it is valid. 6 Dimensional reduction literature interchangeably refers to the thermal scale as hard or superheavy and the intermediate scale as soft or heavy.

Scale Validity Dimension Lagrangian Fields Parameters
Hard Step 1: Integrate out n = 0 Matsubara modes Table 1: Conventional dimensional reduction of (d+ 1)-dimensional Abelian Higgs model into effective d-dimensional theories based on the scale hierarchy at high temperature [64]. The effective couplings are functions of the couplings of their parent theories and temperature and are determined by a matching procedure. The first step integrates out all hard nonzero modes. The second step integrates out the temporal scalar B 0 with soft Debye masses m D . Lattice studies of thermodynamics in this model [111-113, 118, 119] were performed at the ultrasoft scale, where only ultrasoft spatial gauge fields B r (with corresponding fieldstrength tensors F rs ) remain along with a light Higgs that undergoes the phase transition. In our perturbative computation for bubble nucleation, Φ 3 is not assumed to be ultrasoft and hence we do not perform the second reduction to the ultrasoft scale, but instead construct nucleation EFT; see Tab. 2.

Scale Validity Dimension Action Fields
Intermediate on the critical bubble are not enhanced by a scale hierarchy. The leading-order effect from these scales is encoded in the fluctuation determinants.
Lastly, we note that there are scale-shifting fields. They belong to the intermediate scale on the main body of the critical bubble and to the nucleation-and lower scales on the bubble tail. These scale-shifter fields include all the three dimensional fields, except B 0 and the nucleating d.o.f. φ 3 . Contributions related to scale-shifters entail subtleties, but they only appear at higher orders [59].

Matching from intermediate scale to nucleation scale
The matching to the nucleation scale is the second step in Eq. (3.12); see Tab. 2. It proceeds with matching the 1PI actions for the nucleating field [59]. This matching procedure differs from the dimensional reduction from thermal to intermediate scale. The latter can be performed in the symmetric, unbroken phase perturbation theory. The difference arises since the background of the nucleating field, φ 3 , affects the masses of the soft scale fields at leading order. The nucleation scale spatial variations of the background can still be treated as small external momenta in the loop integrals, and can hence be treated with a gradient expansion.
The matching to the nucleation scale EFT reproduces the two first orders from the conventional computation of the nucleation rate with the full effective action; see the description in Sec. 2.2. This is because the scale hierarchy enhances the two first orders coming from the intermediate scale. By inspecting Eq. (1.2), we see that terms enhanced by negative powers of g result from the intermediate scale (and implicitly from the thermal scale), compared to the contributions from the nucleation scale, encoded by the (ln A)-term that is of O(1) up to possible logarithms of g. Thus, one can interpret the two first orders as creating an effective description for the nucleating field. Now we match to the nucleation scale by matching the 1PI action for the φ 3 -field in analogy to a scalar field theory [59] and SU(2) + Higgs theory [67]. 7 The matching helps to understand better the physical picture coming from scale separations. Formally, we equate actions at the intermediate and nucleation scales which is illustrated diagrammatically  [120] for an alternative derivation) which foreshadows which terms contribute to the LO (NLO) action B 3d 0 (B 3d 1 ). This discussion also applies for the field renormalization factor Z which only obtains contributions from the intermediate scale in the matching. This detail ensures that a derivative expansion for the action applies.
To understand how Eq. (3.14) arises, let focus on contributions (a)-(c) of the nucleating field in Eq. (3.13). The EFT construction requires the IR behavior on both sides of the equation to be equal. In other words, the Higgs field contributions collected by (a) on both sides of the equation are equal (with minor subtlety related to diagrams (b) and (c) that are discussed below). Concretely one can achieve this with dimensional regularization and a strict perturbative expansion [63,92], where the mass of the nucleating field is treated as a perturbation which renders its propagator massless. As a consequence, all loop contributions with only Higgs fields vanish as they give rise to scale-free integrals. Therefore, the only contribution remaining on the l.h.s. is the "tree-level" leading order nucleation scale effective action. Eventually, the resulting action S nucl contains non-polynomial and cubic terms, which are induced since the bubble background contributes to the leading-order in the masses of the other 3d fields; see Eqs. (B.2)-(B.6). Furthermore, in strict perturbation theory the mixed bubble diagram (c) vanishes.
The case of mixed sunset diagram (b) is slightly more subtle. Therein, the momentum of the Higgs propagator can be either at the intermediate scale, or at the light nucleating scale. In the EFT description, this division can be depicted as The first term after the arrow corresponds to a case where the Higgs propagator is at the intermediate scale and contributes to the matching in Eq. (3.14). From the perspective of the nucleation scale EFT this is a local contribution to S nucl , which is depicted as diagrams shrinking to a point. In the second term, the Higgs propagator is at the nucleation scale, only the heavy loop shrinks, and in the EFT this is a resummed one-loop contribution. This IR contribution is formally already included in (a) and in fact vanishes in strict perturbation theory, as described above. Similarly other resummations are automatically implemented in Therefore, it should not contribute to the matching which is indeed the case since it vanishes by construction in strict perturbation theory.
We have computed all different diagrams for Eq. (3.14) in Appendix B and collect here the result where different terms correspond to formal expansions 8

(3.21)
At leading order where the last term is the one-loop contribution of spatial gauge fields and the temporal scalar B 0 , which provides a barrier between minima of the potential. We emphasize that the LO potential is gauge invariant.
Here, we observe the benefit of the higher order thermal resummations. The first term with µ 2 3 includes O(g 4 ) resummations (which arise at two-loop level, cf. (A.24)), whereas µ 2 eff in Eq. (2.63) is only resummed at one-loop level, which is correct at O(g 3 ). A similar discussion applies to λ 3 and other 3d EFT parameters. Since 3d EFT parameters are dimensionful, the proper power counting to organize the perturbation theory does not directly follow from Sec. 2.2. We postpone the discussion on the power counting to the next section, but highlight here that one should not simply expand in powers of g, which would compromise the benefit .

(3.24)
Here we complement the comparison to Ref. [69] below Eq. (2.71). The field renormalization factor contributes to B 3d 1 as where we used the characteristic length ∼ µ −1 eff in analogy to Eq. (2.70). This contribution originates within a region of the critical bubble, which has a field value close enough to the broken phase, i.e. within the characteristic radius µ −1 eff . Within this radius, the field value has the power counting φ 2 3 ∼ T (cf. Eq. (2.54)). Outside the characteristic radius of µ −1 eff , the field value of the critical bubble becomes parametrically smaller, φ 2 3 ∼ gT ; see Fig. 1. In this region, Eqs. (B.65) and (B.66) no longer agree to O(g). However, contributions from this region are suppressed by the field value (cf. Eq. (2.74)), and they contribute to either O((ln g) 2 ) or O(1) [59], which are beyond NLO accuracy for B 3d 1 . Furthermore, the derivative expansion breaks down in this region. Since the field value is φ 2 3 ∼ gT , no scale hierarchy protects the derivative expansion.
Before discussing gauge invariance of the nucleation EFT action, we summarize the effect of the matching. The statistical part of the nucleation rate reads formally The EFT description is physically intuitive as one obtains a consistent statistical description for the length scale of the nucleating bubbles by integrating out the shorter length scales. In the nucleation EFT, it is possible to find the leading-order critical bubble before computing the one-loop contributions at the nucleation scale, that are described by functional determinants in Eq. (3.11) and account for the critical bubble background. Furthermore, two-loop contributions of the nucleation scale could also be pursued. We have schematically depicted them above in Eq. (3.26), but leave their determination to future work, together with the computation of one-loop determinants. In such two-loop computation also one-particle reducible (1PR) contributions appear when expanding around the saddle point of the nucleation scale action, S nucl . The 1PR dumbbell diagram in Eq. (3.26) follows from radiative corrections to the critical bubble from the nucleation scale.
For a more detailed discussion on 3d EFT for bubble nucleation, see [59] as well as [67] for technical details on computing higher order contributions.

Gauge invariance of the nucleation scale effective action
This section demonstrates the gauge independence of the nucleation scale effective action in Eq. (3.17). We start by discussing the power counting within the 3d EFT. While all 3d EFT parameters are dimensionful quantities, we introduce dimensionless quantities (cf. [121]) by extracting powers of g 3 from 27) and the dimensionful couplings (3.28) Thus, one can monitor the relative sizes of the different terms in the leading-order potential in a dimensionless form concluding they are of equal size when (note that ρ ∼ 1) Therefore, in 3d EFT the proper expansion parameter is x, the dimensionless ratio of the scalar quartic coupling and the gauge coupling squared [64]. We emphasize that the connection between the power counting here and the earlier assumption of µ 2 eff ∼ g 3 T 2 in Sec. 2.2 is encoded in the scaling of the dimensionless ratio y. In the full parent theory, the scaling for µ 2 eff applies for a temperature regime in the vicinity of T c and this can be translated to the 3d EFT via dimensional reduction matching. The 3d EFT, however, can also be studied as an independent entity and for such analyses the power counting of Eq. (3.30) describes a first-order phase transition in the presence of barrier at leading order.
The leading-order potential in Eq. (3.29) is gauge invariant because this power counting enables an expansion in powers of m 2 χ,3 /m 2 c,3 ∼ x, in analogy to Sec. 2.2 and explicit gauge dependence arises at NLO at O(x −2 ). The mass squared eigenvalues scale as Comparing to the previous section and using Eqs.
The resulting scaling for S nucl is B 3d 0 ∼ x −3/2 , B 3d 1 ∼ x −1/2 . The Nielsen functional C (cf. Eq. (B.76)) expands and scales as with Nielsen identity Appendix B explicitly verifies that the above Nielsen identity holds. Consequently, the proof of gauge invariance ∂S nucl ∂ξ = 0 is analogous to the conventional calculation in Sec. 2.2, and we do not repeat the steps explicitly.

Accuracy of the results
We already foreshadowed the parametric accuracy of our computation in Sec. 1 around Eq. (1.2). This section further discusses the accuracy and some limitations of our computation. Following Refs. [46,47], we demonstrate perturbative accuracy of our computation by inspecting leftover renormalization-scale dependence in our results, since this indicates the convergence of the perturbative expansion. The implicit, leading renormalization group running of parameters in LO expressions cancels against explicit logarithms of NLO expressions. A leftover running of parameters inside NLO terms indicates the size of missing NNLO corrections, and NLO running of LO terms.
Due to the EFT setup of our computation, two different renormalization scales emerge which are related to the different scales that have been integrated out in the EFT construction. The 4d scale Λ is related to the hard thermal scale, and appears (explicitly) inside dimensional reduction matching relations and (implicitly) in the running of 4d parameters. Appendix A discusses the cancellation of scale dependence between LO and NLO expressions, that proves renormalization group improvement related to the scale Λ. On the other hand, we also have the 3d renormalization scale Λ 3d related to integrating out the intermediate scale and running of the 3d mass parameter (cf. (B.52)). 9 We can gain insight in the perturbative accuracy within the 3d EFT by the following analysis: We express the leading contributions to ln Γ using the effective action ∆S We can then analyse the 3d renormalization-scale dependence of this expression, keeping in mind that the leading-order bounce solution φ b will depend implicitly on the renormalization scale through its (again implicit) dependence on the parameter µ 2 3 (y in the nondimensionalized theory, cf. Eq. (3.28)). But as it turns out this dependence is irrelevant at leading order, In the first line, the term related to the y-dependence of φ b cancels since φ b extremizes ∆S LO . The second line uses the definitions of B 3d 0 and B 3d 1 . To understand the remaining 9 The 3d EFT is super-renormalizable: running of the mass parameter arises at two-loop order and is the exact running -it does not receive any higher order corrections [64]. Furthermore, 3d coupling constants do not run at all and are RG-invariant.
renormalization-scale dependence, we consider the beta function (cf. dimensionful version in Eq. (B.52)): Interestingly, the three terms in the beta function are of different orders (β x 0 y , β x 1 y , β x 2 y ) (terms of higher order in x correspond to terms of higher order in λ 3 in Eq. (B.52)). The first term β x 0 y will cancel the Λ 3d -dependence of the corresponding sunset-diagrams in B 3d 1 (cf. Appendix B and Eq. (B.57)), but the two remaining uncancelled terms correspond to sunsets that appear at higher orders. This indicates that where in the second line we used that y ∼ x −1 and B 3d 0 ∼ x −3/2 . In summary, we demonstrated RG-improvement at the order we work, i.e. S nucl ∼ x − 3 2 + x − 1 2 , with a leftover 3d scale dependence at O( √ x) which exceeds the accuracy of our computation.
There is an interesting parallel to the case of a tree-level barrier that is the focus of [47]. In that case the leading contribution goes as g −1 [59], and this is the only contribution calculable in the derivative expansion. The absence of any sunset diagrams to cancel the running leaves all of the leading running uncancelled. But because this perturbation theory is more well-behaved, the end result is that Λ 3d d dΛ 3d ln Γ = O(g), which is formally better than the radiative barrier case.
Next, we discuss limitations of our calculation when approximating the equation of motion with the leading-order potential and field renormalization in Eq. (2.64). In this context it is useful to consider the phenomenological free energy of the bubble in the thin-wall limit where σ is the surface tension, ∆p = −∆V eff is the pressure difference between the phases, and A b and V b are the surface area and the volume of the bubble, respectively. For a study of the surface tension in electroweak theory, within the 3d EFT, see [122]. The equivalent for finding the critical bubble as extremizing the action (which leads to equation of motion in Eq. (2.64)) is to extremize the free energy with respect to the bubble radius. This results in , (3.42) for the radius of the critical bubble, where we expanded the pressure in powers of g. The surface tension, σ, should be expanded in similar manner. To simplify the discussion, we merely keep its leading behavior σ ∼ g 3 , that comes from V g 3 . Inserting Eq. (3.42) into (3.41) yields where the term linear in p g 4 vanishes due to the extremality of the critical bubble and the first correction appears at O(g 2 ). At first sight one could think that this indicates that NLO corrections to the equation of motion in Eq. (2.64) are indeed negligible. However, p g 3 goes through zero in the thin-wall limit and therefore the ratio p g 3 is not necessarily suppressed and the expansion in Eq. (3.42) breaks down. This indicates that our computation is invalid in the vicinity of thin wall limit. At the end of this section, we revisit this exercise and highlight that the thin-wall limit is only relevant for relatively weak transitions.
The formal consistency of our perturbative computation relies on another subtlety. In our prescription, the nucleation rate is obtainable only for T < T c,LO , since the solution for the leading-order bounce φ b is non-existent at higher temperatures. However, it is conceivable that the NLO contributes to the effective potential such that T c,NLO > T c,LO . This does not compromise the computation of equilibrium properties of the transition such as T c or the latent heat. It will, however, pose an apparent problem for computing the nucleation rate. For temperatures in the range T c,LO < T < T c,NLO (cf. Fig. 2) one expects bubble nucleation to occur at NLO description, since the symmetric phase is metastable and the minima are separated by a barrier. Crucially, now the LO bounce solution is non-existent. Hence, the prescription described in earlier sections is not immediately applicable. Nevertheless, this should not be interpreted as a failure of our perturbative setup. The reason for this is that in practice for strong transitions, the relevant temperature range for bubble nucleation is far below the critical temperature, as we explain in more detail below.
On the other hand, the aforementioned issue can be resolved already at NLO based on the freedom to choose the 3d RG-scale. At LO, the effective potential can be tuned by the 3d RG scale Λ 3d since the 3d mass parameter is running in terms of this scale. This way, we can fix a value for which T c,LO = T c,NLO . The corresponding renormalization point is unique since the combined effective potential at LO and NLO, as well as the nucleation rate in our prescription, are RG-invariant. In effect, we can use this observation to resum NLO contributions to the effective potential at LO: such resummation propagates to the LO equation of motion for the field, i.e. to the bounce solution (φ b ) at LO. This exercise provides a handle to understand why the "physical leading order" is not only the LO (which is RG scale dependent), but rather LO and NLO combined (which is RG-invariant).
More formally, the effect of higher order corrections on the critical bubble can be estimated as (to reproduce the exercise, see [76]): where we formally expanded S eff = S eff LO + ∆S eff and φ B = φ b + ∆φ. Here, ∆S eff incorporates higher order corrections. The resulting estimation for the correction to the critical bubble, ∆φ(y), assumes a 3-dimensional x-integration and depicts the leading result as a propagator with a tadpole insertion. While this estimate is computable from LO and NLO actions, its effect to the nucleation rate is of higher order than required in our power counting. By fixing the 3d RG scale such that the LO and NLO critical temperature are equal, the dominant NLO contribution from ∆S eff is effectively resummed into the LO equation of motion for the bounce where S eff LO ≡ B 3d 0 . Consequently the tadpole expansion could partially be resummed into φ b . In the thin-wall bubble regime, the described procedure and RG scale resummation is still compromised close to the critical temperature. However, we can argue that this is not an issue for strong transitions that are interesting for cosmologically relevant applications, for which S eff ∼ 100 [46]. In the thin-wall regime, this condition translates to σ 3 /(∆V eff ) 2 ∼ 6; see Eq. (3.43) and herein S eff ∼ F thin-wall . For the thin-wall approximation to hold, ∆V eff must be sufficiently small, which in turn -to satisfy the previous numerical relation -requires the surface tension σ to be small and conversely the transition to be weak 10 [68]. In physically interesting applications, such as the electroweak phase transition in BSM theories, transitions are strong, and hence interesting temperatures for nucleation lie in a region with supercooling sufficiently far away from T c . The emergent bubbles form with rather thick than thin wall [46,73]. Finally, if the system supercools significantly from the critical temperature before nucleating, the effective mass can become parametrically smaller than g 3 2 T , i.e. ultrasoft (g 2 T ) 2 . In this case, the system becomes non-perturbative in the symmetric phase and the nucleation rate is unattainable with perturbative methods.

Discussion
In what follows, we first summarize the key components of the procedure for achieving a gauge invariant, perturbative computation of the nucleation rate in the presence of a radiative barrier at high temperature. We subsequently discuss the implications for earlier work on nucleation.

Summary
While the EFT procedure for implementing gauge invariance follows the general framework presented in Ref. [59], it is worth highlighting the steps necessary for consistency with the Nielsen identities: Step (1): Implement dimensional reduction to obtain the high-temperature 3d EFT. This step is shown to be gauge invariant in Appendix A.
Step (2): Establish a power counting that yields a gauge invariant, leading-order potential with a radiative barrier arising from parametrically heavier fields. Compute next-to-leading order corrections from the heavy fields to the effective potential and field renormalization factor of the nucleating field, and use them to construct the nucleation scale effective action S nucl = −(B 3d 0 + B 3d 1 ) (cf. Eqs. (3.18) and (3.19)).
Step (3): Use the leading-order bounce equation to obtain the gauge invariant, critical bubble configuration φ b .
Step (4): Compute the leading exponential contribution to the nucleation rate given by ∆S nucl evaluated at φ b (cf. Eq. (3.11)).
Step (1) assumes a high-temperature scale hierarchy, which is valid for the electroweak phase transition and many other thermal phase transitions. If dimensional reduction is performed at LO using effective one-loop masses, tree-level couplings, and a 3d effective potential at one-loop, its accuracy matches the one-loop daisy resummed thermal potential. However, it was recently shown in [46,47] that the determination of equilibrium thermodynamics in this approximation contains large theoretical uncertainties, and in order to reduce them, dimensional reduction at NLO is essential. Furthermore, without a 3d EFT setup that systematizes thermal resummations, it is more cumbersome to compute the two-loop level contributions necessary in step (2).
In step (2), fields such as gauge fields and other scalar fields are parametrically heavier than the nucleating field. This gives rise to an additional scale hierarchy within the 3d EFT, ensuring that it is now justified to integrate out these other fields. Formally this corresponds to an EFT matching between a heavy intermediate scale and a lighter nucleation scale. 11 We emphasize that power counting in step (2) works when T n is sufficiently close to T c .
Step (4) computes the leading exponent term. As we have shown, this computation is gauge-independent. This step can be performed by numerical integration or by using tunneling potentials [84,129].

Implications and outlook
The companion article [1] presented the gauge invariant procedure for thermal bubble nucleation. In the article at hand, we embedded this procedure in the general EFT approach of [59]. This makes it possible to include NLO thermal resummations, as well as to justify how the computation is organized according to the chain of scale hierarchies in Eq. (3.12). In this regard, our work extends the examples considered in [59] to a gauge theory. We have hence demonstrated and established gauge invariance of the nucleation rate computed in this procedure. Furthermore, we showed the theoretical importance of the NLO exponent B 3d 1 since it cancels the leading RG running of B 3d 0 , and therefore illustrates a consistent perturbative expansion.
Our findings on gauge invariance differ from the conclusion of [69]. Therein it is argued that the breakdown of the gradient expansion in the symmetric phase results in an irreducible gauge dependence in a perturbative computation of Γ(T > 0). In fact, this breakdown indeed induces uncertainty to the nucleation rate at NNLO and beyond, and it must be considered when addressing such higher orders of the computation. The leading orders, however, remain gauge-independent as our articles demonstrate. While Ref. [69] mainly focuses on the resummation of infrared effects at finite temperature, dimensional reduction to the 3d EFT is not utilized drawing their approach closer to [1] than this article at hand.
In the future, a non-perturbative computation of the nucleation rate could be performed in analogy to Refs. [65,66,68]. Our perturbative computation provides a sound, gaugeinvariant comparison for such non-perturbative analyses, and allows to further investigate the validity of the perturbative approach. Another interesting future direction would be to apply recent developments [67] to the Abelian Higgs Model discussed in our work and focus on the computation of nucleation determinants without derivative expansion, which allows one to pursue higher order corrections to the nucleation rate.

Acknowledgements
We wish to thank Andreas Ekstedt, Oliver Gould and Anders Thomsen for enlightening discussions. Specially we thank Suntharan Arunasalam for his contributions at the early stage of this project. In addition, we thank the anonymous referee for pointing out the resummation scheme based on the freedom to choose the 3d RG scale, in Sec

A. Gauge independence of the dimensional reduction
This Appendix explicates results for the dimensional reduction matching in Eqs. (3.5)-(3.10) and in particular proves their gauge independence. In this part of the calculation, we use general covariant gauge, or Fermi gauge, defined by the gauge fixing term where the Fermi gauge fixing parameter, ξ ′ , is primed to discriminate it from the R ξ -gauge fixing parameter ξ, in Eq. (2.3). In the case of U(1) gauge symmetry, the ghost field decouples in 4d perturbation theory for any Green's function computation and gives no contribution to the matched EFT parameters. Before focusing on the technical computation, we amend an incorrect statement made in [135]. The latter stated that gauge invariance in the EFT construction is maintained by keeping only Wilson coefficients of O(T 2 ) terms. This corresponds to one-loop in thermal masses and tree-level in couplings which is the LO dimensional reduction and indeed gauge invariant. However, also NLO dimensional reduction -two-loop in thermal masses and one-loop in couplings -is gauge independent: all 3d parameters are constructed from correlators G (Wilson coefficients) and field renormalization factors Z. While individually gauge-dependent, their combination in the matching parameters yield gauge-independent results at NLO dimensional reduction. 12 This important detail is mentioned in [62] but perhaps not stressed enough since it is lacking a technical demonstration. Also an earlier Ref. [64] states that the NLO 3d EFT parameters are gauge invariant. Ref. [136] demonstrates this at one-loop order, but not at full NLO. These works are also referred in [137,138], which describe a gauge-invariant analysis within the 3d perturbation theory, but without the reduction step from 4d to 3d EFT. Recently gauge invariance of the dimensional reduction in the Standard Model was demonstrated [46]. By performing the computation in Fermi gauge, the gauge fixing parameter was shown to explicitly cancel at NLO. Our computation in this section follows the same strategy.
In the four-dimensional computation with dimensional regularization, we use the following conventions for renormalization. These are the definitions of bare quantities in terms of renormalized parameters and counterterms: Here Λ is the 4d renormalization scale in dimensional regularization and counterterms are defined to cancel ultraviolet (UV) divergences of the correlation functions. The renormalization group equations (RGE) and corresponding β-functions for parameters can be obtained by requiring that the bare quantities are independent of the 4d RG scale. At one-loop level they read Λ d dΛ The β-functions and counterterms are zero-temperature objects since finite temperature does not alter the UV structure of the theory. In particular, by separating correlation functions into their soft IR and hard UV parts, UV divergences in hard pieces are cancelled by counterterms in the matching procedure. While suppressing explicit expressions for counterterms, we note that they are closely related to the renormalization group equations above and the field renormalization factors.
In the general context of low-energy effective field theories, Ref. [139] reviews the rationale for dimensional reduction. It discusses the required resummations to remove the hightemperature infrared divergences by matching the correlation functions at the higher scale and lower scale EFT. Ref. [110] presents a practical tutorial for the matching procedure for a real scalar field. 13 Below, we present a formal recipe for this matching. For a generic field ψ, we denote n-point correlation functions by where n > 2. We distinguish the 2-point function Π and expand in soft external momenta K = (0, k) with |k| = k ∼ gT : Here, G denotes the correlator at zero external momenta and Π ′ is the quadratic-momenta correction that contributes to the field renormalization factor Z By matching the effective actions in both theories, the leading (quadratic) kinetic terms yield the relation between 3d and 4d fields where we denote the scalar background fields by ϕ and illustrate the separation into soft (k ∼ gT ) and hard (K ∼ πT ) modes. For simplicity, we omit the field subscript from Z and Π ′ for a moment. By construction of the 3d EFT, contributions Π ′ 3d = Π ′ soft cancelthis is a requirement that the theories are mutually valid in the IR. Therefore, only the hard modes contribute to the last line in Eq. (A.14). By equating the quartic terms of the effective actions, we get 13 DRalgo [140], an automated package for dimensional reduction for generic models was put forward recently.
where we omitted the field subscript from Γ for a moment and also for illustration separated the tree-level part of the correlator from loop corrections. Again by virtue of the EFT construction, terms Γ soft = Γ 3d cancel. After inserting Eq. (A.14) for the field normalization into Eq. (A.15) one can solve for the 3d quartic coupling λ 3 : where we reinstated the corresponding field subscripts, and further marked at which loop order the NLO contributions arise. Terms contributing at next-to-next-to-leading order (NNLO) are neglected. Here, the first line shows the form presented in Sec. 3, and the second and third lines illustrate the Taylor expansion and the composition of the NLO contribution. In this NLO piece, the gauge fixing parameter will cancel between contributions from G and Z, as we will show explicitly below. We have also highlighted that at LO the coupling is an implicit function of the RG scale Λ. Below we show that its running is cancelled by explicit NLO logarithmic terms.
The remaining parameters are matched analogously. However, for the mass parameters µ 2 3 and m 2 D leading contributions arise at tree-level and one-loop, and therefore, the NLO result contains two-loop 2-point diagrams with vanishing external momentum. Also, correlators with gauge field external legs cannot be generated by an effective potential with scalar background field. These correlators require a background field for the gauge field [141] or a direct computation of correlation functions [62].
To compute all correlation functions in Eqs. (3.5)-(3.10), we employ in-house FORM [142] software developed and demonstrated in [46,110]. We employ dimensional regularization and Taylor-expand in soft scales, namely in the external momentum K = (0, k) ∼ gT and the scalar masses µ 2 ∼ (gT ) 2 . We denote spatially transverse projectors and define We consequently obtain the correlators To achieve NLO accuracy, we compute 4-point correlators at one-loop, to zeroth order in the scalar mass, and with zero external momenta.
All above expressions are ξ ′ -independent which completes the proof of gauge independence of dimensional reduction at NLO. Additionally, these matching relations are independent of the renormalization scale Λ. The implicit running of LO terms in Eqs. (A.7)-(A.9) cancels the explicit Λ-logarithms in L b of NLO terms. In fact, the cancellation of ξ ′ in the above matching relations is related to the following observation: all ξ ′ -dependence at NLO is associated with L b pieces and these ξ ′ -dependent logarithmic terms have to cancel identically since they cannot be cancelled by the running of LO terms since β-functions are gauge invariant. The matching for the scalar mass parameter in Eq.(A.24) is subtle as we can replace [62,64] which describes the exact running in 3d (on the r.h.s. we set the hypercharge Y φ = 1); see Eq. (B.52). The contribution in terms of 3d parameters is formally of higher order in power counting but it is natural to include it due to the super-renormalizable nature of 3d EFT: this running receives no higher order corrections. The above matching relations reproduce previous results [64,111] and improve the Debye mass m 2 D to two-loop level. At this point, we would like to address the power counting. Higher order NNLO terms are parametrically of the form Xg 4 , Xλ 2 , Xλg 2 , Xµ 2 λ, Xµ 2 g 2 , where X can be λ, g 2 , µ 2 . 14 This means that our NLO matching relations are missing terms which are at most O(g 6 ). The main interest of this work is a radiatively generated phase transition, for which the relevant power counting is λ ∼ g 3 . This means that terms of λ 2 ∼ g 6 can be discarded in the matching relations Eq. (A.24) and (A.26) for µ 2 3 and λ 3 . After this omission, all matching relations are accurate up to O(g 5 ) and terms g 2 λ and µ 2 λ are suppressed compared to O(g 4 ) terms. This completes the dimensional reduction step.

B. Three-dimensional perturbation theory in R ξ -gauge
This Appendix details the computations within the 3d EFT of Sec. 3. Since we exclusively work in the Euclidean 3d EFT, we use a separate notation and subscript three-dimensional quantities to indicate that they belong to the EFT.

Mass eigenstates
We split the scalar fields as Here, we also shift the gauge fixing background field and treatH 3 as an external auxiliary quantum field that appears only in external legs for different correlation functions in the derivative expansion below. Background fields φ 3 andφ 3 are real and the mass squared eigenvalues read for the Higgs, Goldstone, ghost, gauge field and temporal scalar mass eigenstates, respectively. Here, we keepφ distinct from φ which should only be understood as a notation to separate which terms arise from the gauge fixing Lagrangian (3.4). In our computation below, we do not account for the mixing of gauge field and Goldstone, i.e. we identifyφ → φ.

Propagators
The propagators of scalar fields (denoted generically by ψ = H 3 , χ 3 , B 0 ) and ghosts have the form whereas the gauge field propagator reads where D ij (p, m 1 , m 2 ) is the shorthand notation used below and the transverse projector is defined as the spatial part of Eq. (A.17)

Vertex Feynman rules
The Feynman rules correspond to coefficients for cubic vertices and quartic vertices As an example C HHH equals −(3!) times coefficient of the H 3 3 term in the Lagrangian, and similarly for other vertices. The negative sign arises due to the Euclidean metric and the factorial factor takes care of field combinatorics in our diagrammatic computations below. External momentum dependence and Lorentz structure are implicit in these vertex coefficients. This corresponds to the Lorentz structure δ ij for vertices of two gauge fields which is absorbed in the appearing integral structure of each diagram. Similarly, for a cubic vertex with a single gauge field we have the Feynman rules where an imaginary unit appears due to the Euclidean metric. This sign convention assumes inflowing momenta p and k at the vertex. The computations below take into account the momenta dependence of these vertices in the appearing integral structures, that we list next.

Loop integrals
Before the diagrammatic computation, we list and compute all integrals encountered. As an example we consider a contribution to the 2-point correlator H 3 (k)H 3 (−k) where internal fields are gauge fields. This pure gauge field diagram is calculated at vanishing external momentum, at one-loop level and contributes to the field renormalization factor Z below.
Using the Feynman rules, we obtain where the leading numerical factor is the symmetry factor of the diagram. First of all, Lorentz indices (inside the gauge field propagator D ij ) are contracted by software such as FORM [142] giving rise to scalarized integrals that can be computed using algebraic manipulations. One example is employing d-dimensional rotational symmetry by (p · k) 2 → 1 d p 2 k 2 inside the integrands to reduce to a set of known master integrals. In turn, these integrals are then computed by techniques such as Feynman and Schwinger parametrizations. Since we assume the external momentum to be soft, we Taylor expand the purely spatial integrals in external momenta k. This suffices to extract the quadratic coefficient Z by recursively applying the identity 15 before integration. Consequently, all different terms reduce to the one-loop master integral where γ is the Euler-Mascheroni constant, for which we need a special case Eventually, for the integral in Eq. (B.21), we find in d = 3 − 2ǫ dimensions. Scalar-scalar and scalar-gauge one-loop integrals, that arise in the computation of the field renormalization below, are given by The first few orders of (B.22) yield: where d = 3 − 2ǫ dimensions and the limit ǫ → 0 was taken. The above integrals are finite in dimensional regularization and no 1/ǫ poles emerge. A similar computation of the same integrals in D = 4−2ǫ dimensions produces the field renormalization Z presented recently [84] at T = 0. The vector-scalar integrals VS yield the three structures (B.26)-(B.29) that depend on whether the external leg is H 3 orH 3 . At one-loop level, the effective potential has the master integral At two-loop level, the effective potential contains factorizing one-loop bubbles I 3 1 (B.24) and the sunset integral Once vector bosons are present, more sunset diagrams arise and we need: The computation ofD from Eq. (2.38) requires the following integrals In a similar manner, we could list all integrals appearing in the computation of D from Eq. (2.38). However, due to their large number, we decide to suppress them here.

The effective potential
The effective potential [64] at tree-level reads At one-loop level, the effective potential contains only the master integral J 3 (B.32) and reads At two-loop level, the effective potential formally comprises (cf. [64] as well as [138,143]) where the latter directly corresponds to the β-function for the 3d mass parameter The 3d EFT is super-renormalizable and therefore the β-function describes the exact running without receiving higher order corrections [64]. One can verify that the running of the 3d mass parameter within the tree-level potential in Eq. (B.40) cancels the explicit Λ 3dlogarithms inside the two-loop contribution and provides renormalization group improvement. Accordingly, all such logarithmic terms are ξ 3 -independent, as is the running. The explicit result of the two-loop effective potential is lengthy and its display less illuminating. Instead, we focus on the condition for a first-order phase transition due to a gauge-loop induced barrier at leading order. Accordingly, the one-loop vector boson and temporal scalar contributions have to match the tree-level terms in magnitude Note that 3d couplings have dimension of mass and the 3d field has mass dimension 1/2. Thereof, we identify the leading-order effective potential that is formally of O(g 3 3 ) V 3d eff,LO = The next-to-leading correction arises at O(g 4 3 ) and comprises several contributions. Here we formally count powers of g 3 despite it having non-zero mass dimension. The more careful treatment in Sec. 3 considers dimensionless ratios of 3d parameters, Eq. (3.28), instead. At one-loop level, the difference between contributions of the Goldstone and ghost field is a NLO contribution where m 2 G,3 ≡ µ 2 3 +λ 3 φ 2 3 . Another contribution at O(g 4 3 ) arises from two-loop diagrams. These can be accounted for by setting masses m H,3 → 0 and m χ,3 → m c,3 , which are their leading contributions inside the two-loop computation. Additionally, the Hχχ-sunset diagram is suppressed by λ 2 3 and dropped (by λ 3 → 0). As a result, one obtains 1 (4π) 2 g 4 3 φ 2 3 − 1 + ξ 3 + ln (B.57) Above, the ξ 3 -dependent terms arise partly form the one-loop part (m 2 G,3 ) and the two-loop part from the (VS) topology (−g 3 3 φ 3 /(2π)) and similarly for the B 0 term in the second line.

H 3H3
where the second line holds for a constant field expectation value φ 3 = φ * 3 . We employ the gauge transformation variations δ g 3 Φ 3 = ig 3 Φ 3 , δ g 3 Φ * 3 = −ig 3 Φ * 3 and the variation of the R ξ gauge fixing function (2.3) in three dimensions (cf. (2.37)) with a line prescription as in Fig. 3 and the additional bold leg being the external source J(x) (cf. [69]). The integral reduces to the one-loop master of Eq. (B.23) by writing the rational expression as a sum of terms with minimal denominators. At leading order this becomes Notably, in the 3d EFT, the leading-order Nielsen coefficient C g 3 is independent of the scalar background field φ 3 . This is consistent with our earlier observation that the leading-order Z g 3 is independent of the gauge fixing parameter ξ 3 , so that the Nielsen identity is satisfied. In addition, the second Nielsen identity The computation of D follows similar steps and is diagrammatically illustrated in Fig. 13 of [69]. Due to its length, we do not showcase it further and comment that the same Feynman rules are applied as forD which gives rise to additional integrals that require evaluation. We merely show the final result The functional form of these expressions is same as in [69] albeit here all expressions are those of the 3d EFT. Importantly, at leading order both D andD are of O(g −1 3 ). Hence their exact expressions are irrelevant for our discussion in Sec. 3.3.