Theoretical uncertainties for cosmological first-order phase transitions

We critically examine the magnitude of theoretical uncertainties in perturbative calculations of first-order phase transitions, using the Standard Model effective field theory as our guide. In the usual daisy-resummed approach, we find large uncertainties due to renormalisation scale dependence, which amount to two to three orders-of-magnitude uncertainty in the signal-to-noise ratio for gravitational wave experiments, such as LISA. Alternatively, utilising dimensional reduction in a more sophisticated perturbative approach drastically reduces this scale dependence. Further, this approach resolves other thorny problems with daisy resummation: it is gauge invariant which is explicitly demonstrated for the Standard Model, and avoids an uncontrolled derivative expansion in the bubble nucleation rate.


Introduction
A first-order phase transition in the early universe gives rise to a stochastic gravitational wave background (SGWB) which may be observable today. 1 As a consequence, upcoming gravitational wave experiments open a new window into particle physics phenomenology in the early universe. The frequency window of space-based interferometer experiments, such as LISA, may probe the nature of the electroweak phase transition which typically produce a gravitational wave background peaking in the mHz range [5,6]. A first-order electroweak phase transition is motivated in particular by the baryon asymmetry of the universe (BAU), as it provides the necessary departure from equilibrium [7,8]. Although state of the art calculations of the Standard Model (SM) indicate a crossover transition [9][10][11][12][13], it is straightforward to extend it by new scalar fields [14][15][16][17][18][19][20][21][22][23][24] or effective operators [25][26][27] to catalyse a strong firstorder electroweak phase transition (EWPT). 2 Aside from the electroweak transition, hidden sectors can have dark transitions [38][39][40][41][42][43], which could provide a unique window on a dark sector, which only interacts gravitationally.
In order to accurately predict the SGWB in any of these models, one requires a robust mapping between the observables and Lagrangian parameters. This problem can be partitioned in two [6,39]: the mapping between Lagrangian parameters and thermal parameters; and the mapping between thermal parameters and the SGWB, in particular its peak frequency and amplitude. Only if both these mappings are well understood, can the SGWB serve as complimentary to other probes, such as collider experiments [44][45][46] and direct detection experiments [47,48]. In some cases, the reach of a detector such as LISA may even improve collider constraints on effective operators compared to collider upgrades [27]. However, to reach such precision, the theoretical uncertainties associated with the SGWB predictions need to be under control. With these motivations in mind, this paper will study the theoretical uncertainties for different methods of calculating the thermal parameters of a phase transition.
As has long been recognised, at high temperatures the long-wavelength modes of bosons become strongly coupled [49]. This thwarts the usual perturbative expansion. (At least) two different approaches to resumming perturbation theory have been developed to ameliorate this problem at high-T : daisy resummation and dimensional reduction. We compare and contrast these two approaches, estimating the numerical magnitude of theoretical uncertainties.
It was realised early on [50] that the so called daisy diagrams cause the largest infrared contributions, and should be resummed. Concrete resummation methods were developed and utilised to two-loop order. Dubbed the Parwani [51] and the Arnold-Espinosa [52] resummations, these approaches differ in details though are methodologically similar. The Parwani resummation method allows for a smooth transition to the correct low-T behaviour, but generates unphysical linear terms in the potential which shifts the symmetric minimum away from the origin [53,54]. The Arnold-Espinosa resummation method avoids unphysical linear terms, though by only screening the IR bosonic modes, it fails at sufficiently low temperatures. Various attempts to go beyond these leading contributions, and to resum so called super-daisy diagrams were made (see for example Ref. [55] and references therein), though we do not consider these approaches here.
The idea of dimensional reduction is clearest in the Matsubara formalism. Therein, the equilibrium properties of 4-dimensional QFTs at nonzero temperature, T , are described by fields living on R 3 × S 1 β with the radius of the circle equal to β = 1/T . Phenomena on length scales much longer than 1/T do not see the compact direction, and hence should be describable by a purely 3-dimensional effective field theory. As a concrete method for resummation of perturbation theory, this idea dates back to the 1980s [56][57][58][59]. The 1990s revived this approach [60][61][62][63], developed a simpler method to derive the effective coupling constants of the 3d theory and laid down a generic recipe. An important further development was the use of lattice Monte-Carlo simulations to study the 3d effective theory [9,11,64], which led to the discovery that the electroweak transition in the Standard Model is a crossover [10]. In this work we do not discuss lattice simulations in detail, being mostly interested in theoretical uncertainties of perturbation theory.
A first-order EWPT requires physics beyond the SM (BSM), in particular new Higgs interactions. If the particles responsible for these new interactions are significantly heavier than the electroweak scale, it should be possible to integrate them out at T = 0. The resulting effective theory contains all possible operators of the SM, with higher-dimensional operators suppressed by the heavy scale of these new interactions. This is the SM effective field theory (SMEFT) [65,66]. When truncated at dimension-six and keeping to 3rd generation fermions, it contains an additional 60 independent operators. However, as we are only interested in the electroweak phase transition, which takes place due to symmetry breaking in the Higgs sector, we will restrict ourselves to considerations of the single effective operator where φ is the SM Higgs doublet. 3 Such a dimension-six operator may imply a potential barrier between symmetric and broken phases even before considering loop corrections. In particular, a first-order phase transition can be triggered by a relatively small dimension-six coefficient and negative quartic coupling. However, the additional operator O 6 is merely one of many in a complete basis for the SMEFT. Derivative couplings with kinetic and gauge covariant terms (cf. Ref. [65,66,72,73]) make themselves felt at higher energies and can additionally affect the EWPT. Since O 6 is arguably the dominant higher-dimensional operator in composite Higgs scenarios [26], and expected to dominate in scenarios with extended scalar sectors and large portal couplings, we refrain from considering the full SMEFT basis. This does not change the qualitative scope of our analysis.
The phase transition in this effective field theory was studied previously, in Refs. [25][26][27]74] (see also Ref. [75]). These works adopted the daisy-resummed approach with and without a high-temperature approximation. It was found that as the cut-off scale M is lowered, the EWPT strengthens, and both the critical temperature and nucleation temperature decrease. At sufficiently small M 650 GeV, the transition is strong enough to be observable at nearfuture gravitational wave experiments [27]. For M 550 GeV, the bubble nucleation rate never exceeds the Hubble rate and the phase transition never completes.
The strength of the transition is dictated by the single tunable parameter, M . This makes a convenient model to study theoretical uncertainties, the focus of this article. It is also significant to note that as one varies M , the thermodynamics of the phase transition in this SMEFT qualitatively reproduces that of single-step transitions in scalar extensions of the SM, such as the two-Higgs doublet model (see e.g. Figs. (5) and (7) in Ref. [6]). However, we will not be concerned with the validity of the SMEFT, or the minimal truncation that we consider, for describing any particular extension of the SM.
To date, many studies have examined the reliability of perturbation theory for firstorder phase transitions, recent examples include Refs. [4,20,24,54,55,[76][77][78][79][80][81][82][83][84]. Previously the daisy-resummed and dimensionally-reduced approaches were compared in Refs. [24,83,85], which also include comparisons to lattice simulations. In this paper, we go beyond such previous studies by carrying out a comprehensive study of a wide range of different theoretical uncertainties relevant for these two approaches, focusing on the implications for the gravitational wave spectrum.
By way of example, Fig. 1 shows results for the generated gravitational wave spectrum calculated in the daisy-resummed approach with parameters matched at the Z-pole. The calculations are performed for two different renormalisation scales, µ 4 = T /2 and 2πT . 4 Because any dependence on the renormalisation scale is unphysical, the width of the band in Fig. 1 estimates a corresponding theoretical uncertainty. The theoretical uncertainty is multiple orders of magnitude such that the sensitivity of LISA to this parameter point is completely ambiguous.
Later, we will show that when one includes the running of such parameters in the daisyresummed approach, the scale dependence of the gravitational wave peak amplitude reduces by about an order of magnitude. This compares well both analytically and numerically to the same calculation performed using dimensional reduction at one-loop level. However, there is a systematic difference due to the breakdown of the gradient expansion in the daisy-resummed approach in the calculation of the nucleation rate.
In dimensional reduction, the inclusion of next-to-leading order terms is comparatively amenable, and somewhat standard. It has been argued that these terms are essential for higher order corrections to thermodynamic observables to be perturbatively small [63]. However, the proof really is in the pudding: by explicit calculation, we find that, with the inclusion of these next-to-leading order terms, the theoretical uncertainties in the dimensionally-reduced approach are numerically much smaller than in the daisy-resummed approach.
The outline of this paper is as follows. Section 2 outlines and compares the daisyresummed and dimensionally-reduced approaches at a theoretical level. Daisy resummation is introduced in Sec. 2.1, and the recipe used to calculate thermodynamic quantities is given in Sec. 2.1.2. Following this, Sec. 2.2 shows how higher order resummations are incorporated by dimensional reduction, giving an explicit comparison of the effective potentials as series expansions in the couplings. In Sec. 2.2.2 we give the recipe used to calculate thermodynamic quantities in the dimensionally-reduced approach, including an overview of the -expansion.
In Sec. 3, we show explicit numerical comparisons of different theoretical uncertainties. In particular, we study the importance of scale dependence in Sec. 3.1, gauge dependence in Sec. 3.2, the high-temperature approximation in Sec. 3.3, higher loop orders in Sec. 3.4 and corrections to bubble nucleation in Sec. 3.5. In Sec. 3.6 we also gather existing nonperturbative estimates of the effect of the breakdown of perturbation theory.
Finally, Section 4 summarises our findings and discusses their consequences for the predicted gravitational wave signal. We have endeavoured to make the main of the document intellectually self-contained, though some topics and many detailed results have been transferred to the appendices. Appendix A presents the purely 4d parts of our calculations for the SMEFT. Appendix B provides a hands-on introduction to dimensional reduction, followed by the details of our calculation within the 3d approach for the SMEFT. Appendix C discusses various approximations to the nucleation prefactor in-depth.

Thermodynamics of the phase transition
In studying the thermodynamics of the phase transition, we calculate four thermodynamic parameters that play an important role in the SGWB: the critical temperature T c , the percolation temperature, T p , the inverse duration of the phase transition, β/H p , and the strength of the phase transition, α. In the following, they are defined independently of the methods to calculate them.
The critical temperature for a first-order phase transition defines the temperature at which the free energy of both phases are equal. For homogeneous phases, the free energy density is equal to the effective potential.
To define the percolation temperature, requires first discussing the rate per-unit-volume at which bubbles of the broken phase are nucleated. For first-order phase transitions on cosmological timescales, the bubble nucleation rate takes an exponential, or semiclassical, form, Γ = Ae −Sc . (2.1) Here A is the nucleation prefactor and S c is the Euclidean action of the critical bubble, or bounce. For a thermal transition this is the Boltzmann suppression S c = E c /T [86,87], which we assume throughout. To calculate the nucleation rate from first principles, one can start from the semiclassical result [88,89], Here Z[φ c ] is the contribution to the partition function from the region around the critical bubble, φ c , suitably analytically continued [88]. And Z[0] is the contribution from the region around the high-temperature phase, or false vacuum. The factor V is the volume of space and ω c is a real-time frequency which gives the inverse decay time of the critical bubble. The percolation temperature, T p , is taken to be the temperature for which a fraction 1 − 1/e ≈ 0.63 of the universe has transitioned to the broken phase, following the conventions of Refs. [6,90]. In terms of the action of the critical bubble, this condition can be written as where β is the inverse time scale of the transition For strongly supercooled transitions, the second derivative of the tunnelling action can modify the relation between the tunneling action and the inverse time scale of the transition [91,92], however we will not consider this dynamical effect further, as our focus is on quantum field theoretic uncertainties.
Note that we do not calculate the wall velocity, v w , as it requires a real-time calculation which is beyond the scope of this article. We discard the last term in Eq. (2.3), noting that corrections due to the wall velocity will have small numerical impact on our results as long as v w = O(1) which becomes more likely for stronger phase transitions and even for moderately weak phase transitions, especially if there are no BSM light particles to generate additional friction.
Finally, the appropriate measure of the strength of the transition, α, consistent with the conventions of Refs. [93,94], is determined by the difference in the trace anomaly, Θ, between the two phases, d∆V d ln T , (2.5) evaluated at the percolation temperature. Here V denotes the effective potential of the theory and ∆ denotes the difference between the broken and symmetric phases. The numerator, ρ rad is the radiation density of the high-temperature phase, equal to 3/4 of the enthalpy of that phase, We take the effective number of relativistic degrees of freedom g * = 106.75, equal to its high-T Standard Model value throughout (without e.g. right-handed neutrinos). A recent work [95] reappraised the correct definition of α, and proposed a generalisation of Eq. (2.5). In this article we adopt Eq. (2.5) throughout, justified because we focus on quantum field theoretical uncertainties which will be present regardless of the precise definition of α. Concrete calculations will be carried out in the simplest SMEFT truncation. Its Lagrangian includes the gauge, fermion and Yukawa parts of the Standard Model, and extends the Higgs sector by the single additional operator (1.1), L SMEFT = L gauge + L fermion + L Yukawa + L Higgs , (2.7) These parts, the covariant derivatives D µ , the gauge fields with corresponding field strength tensors, the associated gauge couplings, and ghosts follow the conventions of Ref. [96]. In the following, we will also use c 6 ≡ 1/M 2 for the coefficient of the higher dimensional operator, as it is more convenient to work with c 6 when carrying out Feynman diagrammatic calculations, but M , being related to the energy scale of new physics, aids intuition. For experimentally measured physical parameters we will use the central values presented in Table 1 throughout the paper, taken from Ref. [97]. Our perturbative calculations use the MS-scheme for renormalisation, with the 4-dimensional renormalisation scale denoted byμ. We match experimental results to MS-parameters at 1-loop order, matching pole masses using the full

MS-Parameter Observables
Central Value 1-loop self-energies. This includes momentum-dependent terms additional to those from evaluating the second derivative of the 1-loop effective potential at the minimum.

The 4d approach: Daisy-resummation
The effective potential (the free energy density) encodes the equilibrium properties of a phase transition, such as its character (or order), critical temperature and latent heat. While it is possible to compute the effective potential in perturbation theory, the perturbative expansion at high temperatures suffers from problems at low energies. This is Linde's infamous Infrared Problem [49]. Namely, at high temperatures infrared bosonic modes become highly occupied, enhancing the effective loop expansion parameter for modes with energy E ≪ T , where n B is Bose-Einstein distribution and m is mass of the bosonic mode. At sufficiently high temperatures comparable to m/g 2 , the infrared bosonic modes become strongly coupled. Furthermore, infrared divergences appear at finite loop order: at four-loop order for the effective potential [49]. This means that although the electroweak theory is weakly coupled at zero temperature, massless bosonic modes are nonperturbative at high temperatures and should be treated with appropriate (lattice) techniques. However, it is still possible -and economical -to use perturbation theory as a first approximation in studies of phase transitions. 5 For this reason, this section pedagogically describes a recipe for the purely perturbative analysis of cosmological phase transitions. In particular, it describes how to consistently perform resummations to mitigate the Infrared Problem. We also comment on how to find a nonperturbative solution, after resummations are performed perturbatively in an infrared safe manner.

Resummation at leading order
We use dimensional regularisation in D = d + 1 = 4 − 2ǫ dimensions and the MS-scheme with renormalisation scaleμ. We define the notation P ≡ (ω n , p) for Euclidean four-momenta where the bosonic Matsubara frequency is ω n = 2πnT and This last definition is the Bose(Fermi)-distribution with E p = p 2 + m 2 . In addition, we parametrise the perturbative expansion in terms of the weak gauge coupling, g, and assume the usual power counting for the other coupling constants [63] so that the loop expansion and the expansion in powers of g 2 are equivalent at zero temperature. Due to the nonrenormalisability of the c 6 term, that relation contains an explicit energy scale, denoted by Λ, which should be typical of the low energy SMEFT. At high temperatures we will assume Λ ∼ T .
As an illustrative starting point, let us consider the one-loop correction to the two-point correlator at high temperature. This contributes to the 1-loop thermal mass of the Higgs field zero-mode. For a scalar field with MS-mass parameter m 2 , this correction is of the form (dropping overall symmetry and coupling constant factors) . (2.14) Where I 4 1 (m) is the UV divergent zero-temperature piece and I T 1 (m) is the UV finite but IR sensitive finite-T piece. It is this last temperature-dependent term that leads to problems in the infrared. In order to see this, it is more useful to write this integral in the form , (2.15) separating the soft zero-mode from the hard non-zero Matsubara modes. In fact, at high temperature m/T ∼ g (note that this choice merely parametrises the high-T limit), the mass of the zero-mode scales as ∼ gT while all non-zero modes exceed this, with masses of ∼ πTthis signals a scale hierarchy. Therefore, in the high-T limit, non-zero mode excitations of the thermal plasma effectively screen the zero mode. The zero mode acquires an effective thermal mass m 2 T = m 2 + #g 2 T 2 , where the numerical coefficients, denoted generically by #, depend on the group structure and representation of the fields in question. Physically this thermal mass arises as a screening mass due to the heat bath. Since m 2 can be negative, thermal corrections can trigger a phase transition around where they cancel the zero-temperature contribution, i.e. when the effective mass of the zero-mode becomes ultrasoft (we will expand upon this point later).
Similarly, the gauge field zero-mode is also screened by non-zero modes. Since it is massless in the symmetric phase (or gauge eigenstate basis), its thermal mass is solely dictated by the hard modes and reads m D = #gT . This thermally induced mass is called the gauge field Debye mass, in analogy to Debye screening of the electric plasma. In Appendix B.1 we show the calculation of these thermal masses in detail.
The infrared problem manifests itself when considering higher loop contributions, the so-called daisy diagrams where we have omitted an overall combinatorial factor and replaced λ by g 2 according to its assumed scaling. In these diagrams the hard mode contributions (double dashed lines) screen soft zero-modes (single dashed lines). When the inner loop of Eq. (2.16) is a zeromode, i.e. has a soft momenta P = (0, p) and all N outer loops or petals have hard momenta Q with non-vanishing Matsubara frequencies, this contribution is of order O(g 3 ) for any N . Furthermore, it is IR-divergent for N ≥ 2 in the limit of vanishing mass. For scalar fields and zero-components of gauge fields this IR-problem can be partially cured. The recipe is called daisy resummation. One must first calculate the thermal corrections from the non-zero modes to find the corrected mass of the zero-mode. Then, in computing the contributions of the zero-mode, its mass is upgraded to the thermally corrected mass m T . The following subsection describes this prescription in more detail. Now, we turn to the effective potential, which at one-loop is of the following form, in terms of the background field φ where the one-loop piece is composed of the master sum-integral 6 where m 2 is a φ-dependent mass eigenvalue, and for the full effective potential all mass eigenvalues are summed over with proper coefficients for scalar, gauge and fermion fields. Additionally, in renormalised perturbation theory, there is a term with counterterms that we have omitted for simplicity. For the complete effective potential, see Appendix A.4. Customarily, the sum-integral (2.18) is split into a zero-temperature (Coleman-Weinberg) piece and a temperature-dependent piece (thermal function) ≡T J soft (m) . (2.20) Next, daisy-resummation can be performed by replacing the masses of the zero-modes by thermal screening masses 7 where Π T is the one-loop thermal contribution to the screening mass. By writing we end up with the Arnold-Espinosa type [52] -or ring-improved -resummed effective potential where V T ≃ J T ,b/f and V daisy ≃ J daisy . Note, in order to reach this familiar form where the zero-temperature pieces are separated from thermal pieces, we had to subtract the original soft contribution from the resummed one, in order to avoid double counting. Instead of this form of resummation -that is encoded in the daisy term -we could simply write the resummed effective potential as Technically, this resummation can be achieved by adding and subtracting thermal masses for soft modes in the Lagrangian, such that terms with plus sign contribute to the mass and terms with minus sign are treated as counterterm-like interactions [98,99]. This reorganises the perturbative expansion while the original Lagrangian stays untouched.
with V hard ≃ J hard . This form is equal to the Arnold-Espinosa form in the case where only mass parameters have been resummed. However, in the special case of the SMEFT, the new six-leg vertex introduces qualitatively new features to resummation. In fact, the leading so-called flower contributions of the dimension-six coupling c 6 in SMEFT (cf. Eq. (2.14) for notation) appear at 2-loop order for the mass parameter µ 2 h and at 1-loop order for the scalar selfcoupling λ. Their effect incorporates thermal screening by resumming not only the mass parameter but also the self-coupling. Appendix A.4 fully derives the one-loop effective potential in the SMEFT with leading order daisy resummations.
Finally, let us comment on gauge invariance. As Secs. 2.2.2 and 3.2 explain in detail, in perturbation theory a gauge-invariant treatment requires an -expansion, in which the effective potential is expanded around its tree-level minimum. However, the Arnold-Espinosa type resummed effective potential -with an inclusion of thermal corrections in the soft parts -reorganises the perturbative expansion and departs from the strict -expansion, since the thermal correction Π T is of order O( ) [76]. In Ref. [76], a prescription to cure this problem to ensure gauge-invariance has been proposed, and we will comment on this proposition in Sec. 2.2.1. We therefore do not implement the -expansion in our 4-dimensional approach for computing thermodynamic parameters. As such, our 4d analysis retains an unphysical gauge-dependence which leads to a theoretical uncertainty we calculate in Sec. 3.2. For a recent introductory review of daisy resummation, see Ref. [4].

Daisy-resummed recipe for thermodynamics
Here we outline the calculation of the thermodynamic parameters, which are calculated from the effective potential in 4-dimensional perturbation theory to one loop, including both its scale and gauge dependence. In way of summary, a brief recipe of the approach follows: • Fix the zero-temperature MS-parameters by matching to physical observables at the input scale, here m Z , and then run them to a scale characterising the phase transition, e.g.μ ∼ T . Optimiseμ according to the principle of minimal sensitivity [100].
• Calculate the effective potential of the 4d theory by summing the tree-level potential, the zero-temperature Coleman-Weinberg piece and the finite-temperature piece, with daisy resummation.
• Numerically find the minima of the real part of the effective potential to determine the phase structure and pattern of phase transitions.
• Solve the bounce equation with the potential given by the real part of the effective potential. From this, solve Eq. (2.4) to find the percolation temperature T p and the inverse duration of the transition, or β/H p .
The first step in the daisy-resummed recipe consists of zero-temperature physics, and hence is the same as for the 3d approach. For the SMEFT, the details are given explicitly in Appendix A.1. The phases, distinguished by different Higgs vacuum expectation values, are found by numerically minimising the real part of the effective potential, with respect to the background field φ. The imaginary part of the effective potential can be related to the growth rate of long-wavelength modes about a constant background field [101]. We treat the presence of this nonzero imaginary part as a source of systematic uncertainty in our daisy-resummed calculation. Following standard practice in the literature [26], we content ourselves with checking that the imaginary part is much smaller than the real part of the effective potential at its minima. The nonzero imaginary part of the effective potential, which has to be removed by hand in Eq. (2.27), gives a hint that something is not right. The interpretation of this imaginary part as a decay rate [101] does little to allay this suggestion, as this decay rate is not exponentially suppressed and hence is generically a much faster process than bubble nucleation. Further, this decay rate may be nonzero at the broken minimum solved numerically using Eq. (2.27), suggesting that the broken phase itself will decay into another phase with nonhomogeneous Higgs vacuum expectation value (vev). Both this problem and the problem of gauge dependence are circumvented in the -expansion, which leads to a real effective potential, but this method unfortunately is incompatible with daisy resummation.
To calculate the rate of bubble nucleation, and in particular, the effective tunneling action, we assume O(3) symmetry and solve the bounce equation with boundary conditions, This approach essentially follows Ref. [87]. Equation (2.28) is typically solved using the shooting method, here we employ AnyBubble and BubbleProfiler [102,103]. 8 Evaluating the Euclidean action on this solution then yields S c (T,μ), from which the thermal parameters can 8 Very small differences resulting from these different methods are at the percent level and as they are not quantum field theoretic uncertainties. We do not present these differences here. Where inconsistent, we take the geometric mean of results.

The 3d approach: Dimensional-reduction
Dimensional reduction is a general framework for studying the thermodynamic properties of quantum field theories at high temperatures. It applies widely, and has been particularly fruitful in application to non-Abelian gauge theories. While its use in hot QCD is standard and by now approaches impressive orders in perturbation theory [106,107] (cf. Refs. [108,109] for reviews), its success within electroweak theories and studies of EWPT is far less exploited -even though it proved essential in understanding the phase transition of the Standard Model [10,13,63] and various bulk thermodynamic properties therein [110,111]. Despite featuring in early studies of supersymmetric extensions of the SM [112][113][114][115] and of the two-Higgs doublet model (2HDM) [116], only more recently has the use of dimensional reduction in cosmology been reinvigorated, in studies of the SM with extended scalar sectors, such as the real singlet extension (xSM) [83,96], the real triplet extension (ΣSM) [117,118], and the 2HDM [24,119,120]. High-temperature dimensional reduction (DR) is based on a hierarchical separation of energy scales. In accordance with the effective expansion parameter (2.10), the underlying scales render the theory perturbative at the hard scale (p ∼ πT ), barely perturbative at the soft scale (p ∼ gT ), and non-perturbative at the ultrasoft scale (p ∼ g 2 T ). Here p = |p| denotes a momentum scale of particles in the heat bath. Note, that related literature [63] interchangeably refers to the hard scale as superheavy, the soft scale as heavy, and the ultrasoft scale as light. This hierarchy classifies degrees of freedom when constructing an effective field theory (EFT) for its ultrasoft sector; see Table 2. In the Matsubara formalism of thermal field theory the hard scale screens the purely spatial (static) zero-modes which live at the soft scale. At sufficiently high temperature the infinite tower of non-zero modes is integrated out in a conventional EFT sense. This includes all bosonic non-zero modes and all fermionic modes. Their effect and temperature dependence is encoded solely in the parameters of the resulting EFT of lower dimension. Due to the heat bath breaking Lorentz invariance for temporal gauge fields, the 3-dimensional EFT contains temporal remnants of gauge fields that are adjoint Lorentz scalars (A 0 , B 0 , C 0 ). They get screened at the scale of their respective Debye masses m D , m ′ D , m ′′ D ∼ O(gT ). Furthermore, since the spatial gauge bosons are only Debye screened at the next natural order O(g 2 T ), an additional scale separation emerges between 9 Note that this expression for the nucleation prefactor is a rough guess based on the results of Refs. [104,105]. The expression, however, does not reproduce the temperature dependence of the prefactor derived in Refs. [104,105], nor is applicable beyond the parameter point for the SM with light Higgs studied in Ref. [105]. Further, Appendix C shows that Refs. [104,105] contain a significant error in their result for the (statistical part of the) prefactor. Regardless, as we argue in Sec. 3.5, even the definition of the prefactor is problematic in the daisy-resummed approach, so we adopt this estimate nevertheless.
Validity Dimension Lagrangian Fields Parameters h , λ, c 6 , g, g ′ , g s , g Y   Integrate out n = 0 modes and fermions End: d-dimensional Pure Gauge Table 2. Dimensional reduction of (d + 1)-dimensional SMEFT into effective d-dimensional theories based on the scale hierarchy at high temperature. 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 non-zero modes. The second step integrates out the temporal adjoint scalars At the ultrasoft scale, only ultrasoft spatial gauge fields A r , B r , C r (with corresponding field-strength tensors G rs , F rs , H rs ) remain along with a light Higgs that undergoes the phase transition.
the soft scale of adjoint temporal scalars and the ultrasoft scale. The effective theory of the ultrasoft scale is then non-perturbative, since g 2 n B ∼ O(1). Note that, massive bosonic scalar fields may assume all three scales depending on their zero-temperature mass.
The separation of scales defines the high-temperature regime and generically holds for phase transitions involving scalars in weakly coupled theories. With decreasing temperature the zero Matsubara modes of the scalars signal the absolute instability of the high-temperature phase, below some temperature T 0 . A scalar thermal mass, of the general form m 2 T = m 2 + #g 2 T 2 , goes through zero at the temperature This is generically in the high-temperature regime, at least regarding the scalar field undergoing the transition, because the temperature is larger than the vacuum mass parameter by a factor of 1/g. Note also that T 0 < T p < T c , as the thermal mass is necessarily positive at both T p and T c , so implying that T p and T c are generically in the high-temperature regime. This further suggests that bubble nucleation will almost always take place via a purely spatial, O(3) and not O(4) symmetric, instanton [87,121,122]. In contrast to the scalar undergoing the phase transition, for the temporal gauge bosons the lack of a (negative) vacuum mass implies that they are always of the soft scale, and hence are integrated out in constructing the EFT of the ultrasoft scale. The Higgs zero Matsubara mode is treated as ultrasoft throughout our analysis of the SMEFT. At temperatures relevant for the dynamics of the phase transition, between the percolation and critical temperatures, the thermal mass of the scalar zero mode is positive. But following Eq. (2.32) it is at most of order gT , and hence is either of the soft or ultrasoft scale. An ultrasoft Higgs mass is certainly correct in the vicinity of T 0 where the vacuum and thermal mass contributions exactly cancel, but should also hold near T p and T c due to a remaining partial cancellation of vacuum and thermal mass contributions. Regardless, we do not expect any significant discrepancies between treating the Higgs as soft versus ultrasoft due to the small numerical effects of the temporal gauge fields.
The philosophy of dimensional reduction is to treat perturbative modes perturbatively and nonperturbative modes nonperturbatively. Fermions and bosonic non-zero Matsubara modes are perturbative, and are treated perturbatively when integrated out in the construction of the EFTs. The bosons of the soft scale are also perturbative, and are treated similarly. Since only the ultrasoft scale is nonperturbative this scale is then normally treated with non-perturbative lattice studies. Existing lattice studies utilise the super-renormalisability of the EFT to perform an exact mapping between bare lattice parameters and MS-parameters [123,124]. However, in the EFT we consider for the SMEFT, the presence of the marginal, sextic Higgs field operator O 6 means that the EFT is merely renormalisable and not super-renormalisable, aggravating the matching of lattice parameters to known physics. Nevertheless, recent lattice computations in scalar-extended BSM models [24,118] have indicated that, for relatively strong transitions in weakly coupled theories, two-loop perturbation theory within the ultrasoft EFT describes the phase transition with reasonable accuracy; see also Sec. 3.6. There are a few reasons for this perhaps surprisingly good agreement between lattice and perturbation theory. On the one hand, by constructing this effective theory for the ultrasoft modes, dimensional reduction makes it easier to hone in on these important modes and to treat them to higher loop order than is otherwise possible. On the other hand, at least in the case of a strong transition, the transition depends most strongly on the scalar sector, which is, in a concrete sense, less nonperturbative than the spatial gauge bosons, for which there are true IR divergences in the symmetric phase at finite loop order. Further, the IR divergences of the spatial gauge bosons only arise at higher loop order, for example, at four-loop order for the free-energy. So when the first few terms of the loop expansion converge well, one can expect the nonperturbative effects to be relatively small.
In practice dimensional reduction is performed along the modern EFT recipe. One first identifies the most general Lagrangian that respects the symmetries of the full theory, and then matches static Green's functions to determine the parameters of the EFT in terms of temperature and parameters of the parent theory. For a fuller explanation of dimensional reduction, we refer to our Appendix B, which accounts step-by-step of how to construct such effective theories in phenomenologically relevant models. Therein, Appendix B.1 present a breakdown of the calculation of the SU(2) Debye mass, which we hope suitably introduces the nitty-gritty of dimensional reduction. Appendix B.2 presents our explicit results for the dimensional reduction of the SMEFT at full NLO.

Resummations at higher orders and gauge invariance
In dimensional reduction, higher order resummations are systematically incorporated order-byorder in powers of the couplings. This is achieved by careful power counting, necessary because thermal screening breaks the alignment between the loop and coupling expansions. 10 By contrast, in the 4d approach, resummation is carried out in a more ad hoc way, by identifying and resumming infrared sensitive parts at the level of Feynman diagrams. As has long been recognised [52], at higher orders it is necessary to resum new classes of diagrams beyond just the daisy diagrams.
One-loop daisy resummation, as presented above, generates the effective potential accurately up to O(g 3 ). However -as argued in Sec. 2.2 in Ref. [63] -one must go beyond this and achieve O(g 4 ) accuracy in order to obtain perturbatively small fractional corrections to many infrared observables. Further, the RG running of the leading order effective potential starts at O(g 4 ), so one must reach this order to control the RG scale dependence. This requires two-loop contributions, both to the effective potential and to the resummed thermal masses. It also requires additional resummations at one-loop order: both to the couplings and to the field itself, the latter due to the momentum dependence of thermal screening at O(g 4 ). Dimensional reduction provides a systematic means to keep track of these disparate resummations, and is extendable to still higher orders.
The effective potential provides a convenient means to show how the differences between the 3d and 4d approaches manifest in concrete calculations. Schematically there is a relation of the form which holds up to O(g 3 ). Note that at leading order in powers of g 2 , the 3d and 4d fields are related as φ 3d = φ 4d / √ T . At higher orders, momentum dependent thermal screening modifies this relation, as captured in the 3d matching relations. Here, for simplicity, we compare the 3d effective potential at the soft scale, leaving discussion of the effects of integrating out the soft scale to later in this section.
To understand in more detail where the two approaches differ, we break down Eq. (2.33), giving From the construction of the dimensionally-reduced EFT, one can deduce the following approximate equality for the hard contributions This follows since the effective potential is the generator of one-particle irreducible (1PI) correlation functions and 3d parameters are defined by matching the 1PI correlation functions to the 4d theory. Utilising φ ∼ T for the dimensionful background field we arrive at the following schematic power counting, As we indicate, this equation is free from nonanalytic dependence on g 2 because it involves only hard modes. Our 4d approach correctly captures only the leading order term in this expansion. Both this leading term and the O(g 4 ) term are captured in the NLO matching relations of the 3d approach. Appendix B.2 presents the NLO matching relations for the SMEFT. Daisy resummation is engineered to correctly describe the leading effects of the soft scale. This results in the following approximate equality for the remaining soft parts In the soft sector, the presence of infrared modes leads to a nonanalytic dependence on g 2 , As we have explicitly verified in the SMEFT, the 4d approach correctly reproduces the O(g 3 ) term. In the 3d approach, by including two-loop corrections to the effective potential, we capture also the O(g 4 ) term. Appendix B.4 yields an expression for the 3d two-loop effective potential in the SMEFT. To compute the full O(g 5 ) term requires a three-loop computation [109,125], whereas the O(g 6 ) term is nonperturbative [49]. The comparisons made in this section utilise the 3d effective potential of the soft sector. However, to simplify the thermodynamic calculations in our 3d approach we integrate out the soft temporal bosons and instead utilise the 3d effective potential of the ultrasoft sector; see the end of Appendix B.2.2. This additional step incorporates both the O(g 3 ) and O(g 4 ) effects of thermal fluctuations of the soft sector into the parameters of the ultrasoft theory. A difference does arise, though, regarding the dependence on the Higgs vev of the masses of the temporal gauge bosons. For the SMEFT with the 3d EFT truncated at (φ † φ) 3 3d , this difference arises at O(g 3 (φ † φ) 4 3d /T ) for the 3d potential. Although formally of O(g 3 ) for a transition with φ ∼ T , this discrepancy is accompanied by a sizable numerical suppression, O(10 −6 ), due to combinatorial and loop factors. Thus it is expected to have only a very small numerical effect, though it could become more significant for very strong transitions.
Beyond aiding higher order computations, an additional benefit of the 3d approach, is that one can achieve exact order-by-order gauge invariance by applying the -expansion inside the 3d effective theory, c.f. Sec. 2.2.2. In this expansion, the value of the effective potential is computed as an expansion around the minimum of V 3d tree . This possibility depends upon the gauge invariance of the 3d matching relations. In Appendix B.2 we show this explicitly for the dimensional reduction of the SMEFT: choosing a general covariant gauge, the ξ i -dependencies cancel duly in the matching relations up to O(g 4 ) (and higher order terms can be discarded).
Reference [76] -in order to maintain gauge invariance -proposed an alternative resummation approach. In this the soft part of the LHS of Eq. (2.34) is evaluated at the (temperature dependent) minimum of V 3d tree 11 , but the remaining tree-level and hard parts are evaluated in an expansion around the (temperature independent) minimum of V 4d tree . This approach differs from the approaches presented here already at leading O(g 2 ) order, since the minima of the 4d and 3d tree-level potentials differ at leading O(g 0 ) order. Had the authors used the 3d minimum also in V 4d tree + V 4d hard in addition to the soft, resummed part, their potential would have matched those presented here up to O(g 3 ).

Dimensionally-reduced recipe for thermodynamics
Here we outline the calculation of the thermodynamic parameters in the dimensionally-reduced approach. In way of summary, a brief recipe of the approach follows: • Fix the zero-temperature MS-parameters by matching to physical observables at the input scale, here m Z , and then run them to a thermal, matching scaleμ ∼ πT . Optimisē µ according to the principle of minimal sensitivity [100].
• Carry out dimensional reduction at this thermal scale, by matching the static, infrared correlators of the 4d theory to an effective 3d theory. Further match to a reduced effective theory at the infrared energy scale g 2 T /π. This amounts to integrating out all modes with energies ∼ πT and ∼ gT .
• Calculate the effective potential of the effective 3d theory, and find its minima. If possible, maintain a strict -expansion.
• By taking derivatives of the 3d effective potential with respect to the parameters of the theory, calculate the gauge-invariant condensates, such as φ † φ . From this one can find the critical temperature T c , as well as the strength of the transition.
• Calculate the bubble nucleation rate in the 3d effective theory. If possible, maintain a strict -expansion. From this, solve Eq. (2.4) to find the percolation temperature T p and the inverse duration of the transition, or β/H p .
The first step in the dimensionally-resummed recipe is the same as for the daisy-resummed approach. For the SMEFT, the details are given explicitly in Appendix A.1.
In the second step of our recipe, the hard and soft modes are integrated out and the 3d effective theory for the ultrasoft modes is constructed. An explanation of this procedure at a synoptic level has been given above, at the beginning of Sec. 2.2. An example application is worked through in Appendices B.1 and B.2. In application to the SMEFT, we utilise the high-temperature approximation, partially for simplicity and partially because we expect this approximation to be valid, following the argument given around Eq. (2.32). Note, however, that this approximation is not an inherent limitation of dimensional reduction [54,96,126].
Once we have arrived at the ultrasoft 3d effective theory, the advantages of dimensional reduction manifest themselves. The 3d effective theory is simpler than the full 4d theory. Not only has the theory a reduced field content: all fermions, plus any bosons with masses of order or greater that ∼ gT have been integrated out. More importantly, all sum-integrals have been evaluated and one is effectively studying the zero-temperature vacua of a 2 + 1 dimensional theory. The temperature merely enters the parameters of the effective theory. As a consequence, perturbative calculations within 3d can be performed as a vanilla loop expansion, so 3d loop orders are not mixed. This allows one to perform strict -expansions, thereby maintaining order-by-order gauge invariance (see Sec. 3.2), as well as avoiding doublecounting in computing the bubble nucleation rate (see Sec. 3.5). It also makes it more feasible to go to higher loop orders than would otherwise be practical.
The third step of the dimensionally-reduced recipe calculates the effective potential in the 3d effective theory. For our calculations in the SMEFT, we carry this out to 2-loop accuracy. The calculation utilises previous work of Refs. [60,127,128], and is given explicitly in Appendix B.4.
For the SMEFT, due to the presence of the c 6,3 (φ † φ) 3 term in the Lagrangian of the 3d EFT, there is a first-order phase transition at tree-level. 12 From the perspective of the 3d effective theory, the transition takes place as the effective parameters change with temperature. At least for transitions that are not too weakly first-order, 13 higher loop orders will not change the order of the transition, so a strict -expansion should converge well. On the other hand, without the c 6,3 term some one-loop contributions must compete with tree-level contributions to give a first-order phase transition, signalling a breakdown of the -expansion [127,128]. In this case, spurious imaginary parts arise and "little constructive information" can be gained from the -expansion [127] (see also Ref. [128]). The difference is the presence or absence of a tree-level barrier. In the case of a strong first-order transition, one can instead recover a consistent expansion parameter in terms of ratios of couplings (see for example Refs. [52,[129][130][131][132][133][134]) such as λ 3 /g 2 3 for the SM. However, for very weak first-order, second-order or crossover transitions, even this option is no longer possible and one must resort to lattice simulations [9,11,24,118,[131][132][133].
To perform the -expansion in the 3d effective theory, one expands all quantities in powers of , the loop-counting parameter of the 3d EFT. In particular, the 3d effective potential and 12 Note that this is contrary to the full 4d theory, for which the tree-level potential is temperatureindependent. 13 Such that the 3d loop-expansion parameter is perturbative (see Appendix B.3), and loop corrections are small compared with the distance to the second-order point in the phase diagram. scalar vev are expanded as Equations, such as V ′ 3 (v 3 ) = 0 for finding the loop-corrected vev, are then solved order-byorder in . Note that, following the notation in vacuum, we denote by the loop counting parameter, though Planck's constant scales out of the 3d effective theory.
In calculating thermodynamic quantities in the -expansion, one first carries out all the necessary computations at tree-level. Once completed, including higher-order contributions is a simple algebraic exercise in matching powers of and solving corresponding linear equations. The tree-level vev solves, where multiple solutions signify a coexistence of phases. Expanding the effective potential around the tree-level vev, one finds where all the potential terms on the right hand side are evaluated on v 3(0) . Solving for the broken minimum order-by-order in , the solution to The two-loop expression for the broken minimum, Eq. (2.42), contains infrared (IR) divergences in V ′′ 3(1) (v 3(0) ) and V ′ 3(2) (v 3(0) ) due to the vanishing Goldstone mass in the tree-level broken minimum. These Goldstone IR divergences are a feature of the Landau gauge, and do not occur for positive gauge fixing parameters. However, if we regularise the loop integrals by taking the Goldstone mass to zero from above, the IR divergences in Eq. (2.42) precisely cancel. This is equivalent to taking the gauge parameters to zero from above in taking the Landau gauge limit. For a more detailed discussion of this point see Refs. [127,135], and, for an alternative approach see Ref. [136].
At the critical temperature two phases coexist with equal free energy density, which for homogeneous phases equals the effective potential. Thus, to solve for the critical temperature, we have the additional equation to solve, Here, in light of our intention to apply this formalism to the SMEFT, we have assumed that one of the two phases lies at the origin. Following Eq. (2.5), we use ∆ generally to refer to the broken phase value minus the symmetric phase value of some quantity.
As before, Eq. (2.43), is solved order-by-order in . As the 3d EFT does not directly see the parameter T , it makes sense to solve this equation instead for the mass of the scalar undergoing the transition, which is what we do for the SMEFT. For now, let us denote by m 3,c the value of the scalar mass at which Eq. (2.43) holds. The equation determining criticality then takes the form, where the terms on the right hand side are functions of the other parameters in the 3d EFT. Hence this equation defines a surface in this space of the parameters of the 3d EFT. As the temperature changes, a line is traced out in the space of parameters of the EFT, a line which pierces the critical surface at T c . Written more explicitly, the O( ) and O( 2 ) corrections to the critical mass take the form, Here the ∆V 3(i) are evaluated at the tree-level minima, and at the tree-level critical mass. One question then arises: do we also expand T c in powers of ? This has been discussed in Refs. [76,135]. Given that here is the loop counting parameter of the 3d EFT, we have chosen not to expand T c in , instead solving Eq. (2.44) numerically. Both options yield gauge invariant results. Due to the presence of a tree-level barrier between phases in the SMEFT, we expect any difference between these two approaches to be very small, as the difference is formally of higher order in the EFT loop expansion.
The measure of the strength of the transition, α, defined in Eq. (2.5), is calculated via the trace anomaly. Once we have determined ∆V 3 to some order in , following Eq. (2.41) above, the trace anomaly follows from differentiation with respect to temperature. The 3d EFT does not depend explicitly on temperature, so all temperature dependence must arise through the effective couplings. Thus, we find that, where the sum over {κ i } runs over the parameters of the EFT, and the factor of 1/T on the left hand side follows from the basic, dimensional relation between 3d and 4d physics. Note that as ∆V 3 is calculated order-by-order in , this expression is completely independent of the gauge fixing in the EFT [137,138], and, if the effective couplings are themselves gauge invariant, then the whole expression is gauge invariant. Furthermore, the -expansion also ensures that ∆V 3 is manifestly real, so there is no need to ad hoc take the real part, as in the daisy-resummed approach, Eq. (2.27). However imaginary parts can arise in the -expansion in the absence of a tree-level barrier in V 3 , at least when T c is also expanded in powers of [135]. We note that at 2-loop order, a dependence on the 3d renormalisation scaleμ 3 arises. This dependence can be diminished by solving the appropriate renormalisation group equation (RGE) (see e.g. Ref. [60]), though the choiceμ 3 = O(g 3 v (0) ) is sufficient to avoid large logarithms. In Sec. 3.1 we treat thisμ 3 dependence as a source of theoretical uncertainty, and vary it over some appropriate range, to estimate its magnitude. As the dependence onμ 3 only arises at two-loop order, its effect is expected to be small.
The remaining thermodynamic quantities, T p and β/H p , can be determined in terms of the bubble nucleation rate. Just as for the other quantities, this can be calculated in anexpansion within the 3d EFT, and there are significant benefits from doing so. The calculation begins with the semiclassical expression for the bubble nucleation rate, Eq. (2.2), written in terms of the partition function in 4d. The relation between the partition functions of the full theory and the 3d EFT is where f 1 is the coefficient of the unity operator in the dimensional reduction, and (the logarithm of) the equation holds up to some order in powers of the coupling constant (see e.g. Refs. [61,62]). A semiclassical evaluation of the partition function expands around a background configuration. For a background configuration that varies, at most, on the long length scales of the ultrasoft theory, Eq. (2.47) follows directly from the matching of correlation functions carried out in DR. As the coefficient of the unit operator, f 1 , is independent of the field configuration, Eq. (2.2) takes the same form in the 3d EFT as in the full 4d theory, contributes to the partition function of the EFT from the region around the critical bubble, or bounce, and Z 3 [0] contributes from the region around the symmetric phase, or false vacuum. The only factor which cannot be computed purely within the 3d EFT is ω c . In principle it requires a real-time calculation, though it has been estimated in the literature (see Appendix C.1). Since the bubble nucleation rate is an intrinsically real-time quantity, it is quite surprising that all but ω c can be calculated within a timeless EFT. Starting from Eq. (2.48), the calculation of the bubble nucleation rate is now an unambiguous application of the original, vacuum bounce formalism [139,140], except in three Euclidean dimensions. The temperature only enters the couplings of the EFT, so that one only needs to calculate the vacuum tunnelling rate as a function of the couplings. Thus, one can derive from first principles the following tree-level bounce equation for the critical bubble, φ c , of the 3d EFT,  49) contains the tree-level potential in the 3d EFT, and should be contrasted with the 4d approach, in which tree-level and one-loop contributions are mixed in the calculation of the tunnelling action. In the 3d approach, the -expansion of the nucleation rate reads One-loop fluctuations around the critical bubble make up the nucleation prefactor, A, and do not enter S c . The calculation of the nucleation prefactor, the O( ) contribution to the nucleation rate, is, in general, a significantly more difficult task than calculating the tunnelling action, S c . Nevertheless, in the 3d approach, this difficult task is significantly easier than in the 4d approach. This is mostly because loop orders are not mixed in the calculation, and hence a vanilla semiclassical analysis applies out-of-the-box [139,[144][145][146]. The difficult task simplifies further because the temperature has already been eliminated from the calculation, and because the field content is reduced. As a consequence, we are able to reasonably estimate the nucleation prefactor, which we carry out in Appendix C.
Once the rate of bubble nucleation has been calculated, one can solve Eq. (2.3) for the percolation temperature, and evaluate Eq. (2.4) at this temperature to find β/H. Finally, one also evaluates Eqs. (2.46) and (2.5) at the percolation temperature to find α. Nonperturbative corrections unknown unknown Sec. 3.6 Table 3. Sources of theoretical uncertainty and relative importance quantified by the parameter ∆Ω GW /Ω GW defined in Eq. (3.1) over the range M = {580 − 700} GeV in the SMEFT. Although we do not have reliable estimates for the uncertainties of the 4d approach due to higher loop orders and nucleation corrections, they are expected to be much larger than the corresponding uncertainties of the 3d approach (see the relevant subsections).
To determine the relative merits of the two approaches outlined in the previous section, in this section we critically examine a range of sources of theoretical uncertainty. Concrete numerical comparisons are made for the four thermodynamic parameters which are most important for determining the SGWB spectrum: T c , T p , α and β/H p . 15 Further, in order to compare the various sources of uncertainty, we combine these thermodynamic parameters into a single parameter, where Ω max,min are the maximum and minimum peak of the SGWB spectrum due to sound waves (sw) [1,94,147]. These are predicted by varying the thermodynamic parameters across the theoretical uncertainty band and utilising the following, where S sw (f ) encodes the frequency dependence; at the peak, S sw (f ) = 1. We estimate the outstanding factor in Eq. (3.2) , the timescale on which acoustic waves are active, from [147,148] Since we are interested in the variation of the theoretical predictions and not in the magnitude of the spectrum, the latter is of minor importance here. For those theoretical uncertainties which can be assessed numerically, Eq. (3.1) gives a relatively good measure of the corresponding uncertainty in predictions for upcoming GW experiments, such as LISA. We apply this measure in estimating the following theoretical uncertainties: renormalisation scale dependence in Sec. 3.1, gauge dependence in Sec. 3.2, high-T approximation with truncation of the 3d EFT in Sec. 3.3, higher loop orders in Sec. 3.4, and for some corrections to the bubble nucleation rate in Sec. 3.5.
The construction of the dimensionally-reduced theory naturally entangles three different sources of errors: (i) higher-dimensional operators, (ii) loop expansions, and (iii) hightemperature expansion. Concretely, the information from integrating out the hard scale fermionic and bosonic modes is distributed over all three of the above.
Not all the uncertainties that we consider can be reliably quantified. In particular, internal inconsistencies in the 4d approach cannot be estimated by simply varying a parameter. These stem from the relatively ad hoc implementation of thermal resummation in the 4d approach. The most notable of these internal inconsistencies arises in the bubble nucleation calculation. Section 3.5 discusses these, in particular: (i) double-counted degrees of freedom, and (ii) an uncontrolled derivative expansion. Therefore, at the level of principle, the 3d approach should always be preferred over the 4d approach, as it does not suffer from the same internal inconsistencies of the bubble nucleation calculation. A similar point could be made regarding gauge independence, which in principle should be maintained order-by-order in perturbation theory, as is only possible in the 3d approach.
Finally, we should note that any purely perturbative calculation suffers from an irreducible uncertainty due to the nonperturbativity of the IR modes of magnetic gauge bosons in the symmetric phase. We discuss this nonperturbativity in Sec. 3.6, and collect some estimates of its expected magnitude present in the literature. Table 3 summarises all the assessed sources of uncertainty and refers to the relevant sections and figures.

Renormalisation scale dependence
Perturbative approximations to physical results generally depend on the renormalisation scale, signalling a source of theoretical uncertainty in the approximation. Thus, a strong dependence on renormalisation scale, as in Fig. 1, reflects the inadequacy of the approximation 16 . Dependence on the renormalisation scale can be ameliorated by performing renormalisation group improvement [20,60,131], but it also can be exploited to probe the magnitude of higher order perturbative corrections. Without large hierarchies of scale in the problem, all logarithms can be made small by an appropriate choice ofμ, of the same order as the energy scale of the process under consideration. In that case, a variation ofμ by an O(1) factor will induce a correction which is formally of higher order in the perturbative expansion.
When studying a thermal first-order phase transition, it seems natural to choose a 4d renormalisation scaleμ which depends on temperature. A question arises though, as to determine the precise numerical coefficient,μ = #T . In the Matsubara formalism the temperature enters via the thermal frequencies ω n = nπT , with n an odd integer for fermions and an even integer for bosons. As a consequence, one might consider any of the following energy scales as possible choices forμ, T for the n = 0 modes, no factor of π is present, 05T the weighted sum of nonzero bosonic ω n .
The "weighted sums" here are those that arise within logarithms at one-loop order in a range of quantities, such as the free-energy. In a theory such as the SMEFT with a large and varied particle content, any of these choices for the renormalisation scale can be equally motivated. To estimate the magnitude of higher order corrections, we vary the renormalisation scale over an order of magnitude in the rangeμ = (0.5 . . . 2π)T . Figure 2 shows the result of this calculation for both the 4d and 3d approaches. An optimal choice of scale,μ opt can be found according to the "principle of minimal sensitivity" [100]. This principle demands that atμ =μ opt some approximation to a physical quantity, O, is independent of the renormalisation scale, In the 4d approach, choosing O = T c , this equation finds no solution forμ opt , T c (μ) being a monotonic function over the range ofμ considered. In contrast, in the 3d approach, solutions to Eq. (3.5) exist for the SMEFT for the whole range of M that we consider. In general, at finite order in a perturbative expansion, different physical quantities may satisfy Eq. (3.5) at different scales. In the 3d approach we find that for O = T c , Eq. (3.5) is satisfied around In both cases the dependence on M is mild. For O = α and O = β/H, the solution to Eq. (3.5) depends more strongly on M , but is still centred around which we choose as the default renormalisation scale in the 3d approach. Figure 3 shows the scale dependence of T c at the benchmark point M = 640 GeV. As seen in Fig. 2, the scale dependence in the 4d approach (shown as the grey bands) is very significant: about 20-30% for T c , 20-75% for T p , 200-800% for α and 40-200% for β/H. For the gravitational wave peak amplitude, the corresponding uncertainty is ∆Ω/Ω min = O(10 2 − 10 3 ). This shows that higher order perturbative corrections to the 4d approach are large, and the calculation is not well under control at this order. In notable contrast, the scale dependence of the 3d approach (shown as the blue bands) is much smaller: 2-5% for T c , 2-20% for T p , 10-70% for α and 10-60% for β/H. For the gravitational wave peak amplitude, the corresponding uncertainty is ∆Ω/Ω min = O(10 0 − 10 1 ). Thus, at the order that we work to, the 3d approach appears much better under control. The uncertainty in the 4d approach becomes more dramatic when ignoring the scale dependence of the couplings. This is a corner often cut in the literature. An example of the resulting scale uncertainty is shown in Fig. 1 in the introduction.
The 3d effective theory depends on its own additional renormalisation scale,μ 3 . Just as with the 4d renormalisation scale, dependence onμ 3 is unphysical. As the effective theory is applicable to energy scales of order g 2 T ≈ g 2 3 and below, we choose simplyμ 3,opt = g 2 3 as our default scale. Taking this default scale as the geometric midpoint, we vary the 3d renormalisation scale over the range g 2 3 / √ 10 to g 2 3 √ 10. The dependence onμ 3 is in all cases much weaker than that onμ shown in Fig. 2, so we do not plot theμ 3 -dependence explicitly.
The uncertainty related to theμ 3 dependence amounts to ∆Ω/Ω min = O(10 −2 − 10 −1 ), and increases monotonically with M . Note that, due to its numerical insignificance, we do not solve the renormalisation group equations forμ 3 , and neither do we use the more optimal choiceμ 3,opt = g 3 v 3(0) . For more discussion of this, see Ref. [60].

Gauge dependence
In perturbation theory a gauge-invariant treatment requires an -expansion. Therein the Nielsen identities ensure gauge invariance order-by-order [76,128,135,137,138]. By contrast, without performing an -expansion, perturbative results in general bear a residual gauge dependence. We adopt the class of general covariant (or Fermi) gauges, with gauge parameters ξ i . Since daisy resummation generally conflicts with a strict -expansion -as explained in Ref. [76] -in the 4d approach thermodynamic parameters depend on the choice of gauge. This can be traced to the Goldstone modes, which acquire an explicit gauge dependence which we derive for the SMEFT at 1-loop using a general covariant gauge in Appendix A.3. The 3d approach exhibits two sets of gauge-fixing parameters: one within the 4d theory itself and another within the dimensionally-reduced 3d EFT. Our matching relations for the SMEFT are gauge independent at order O(g 4 ), thus implying that one can consistently truncate the relations at this order. However, we choose to include a subset of the O(g 6 ) corrections, in particular the O(g 6 ) corrections to c 6,3 , as these are expected to be relatively numerically important given that c 6,3 is O(g 4 ) at LO (in fact we find this amounts to an O(10%) correction to c 6,3 ). These O(g 6 ) corrections show an explicit dependence on the 4d gauge parameters, though this gauge dependence is numerically very small -even for very large gauge parameters, as is shown in Fig. 4. To remove this gauge dependence would require matching to a complete operator basis in the 3d EFT, and perhaps even a complete matching at O(g 6 ). However, we choose not to, given the difficulty of such a calculation, as well as its numerical insignificance (cf. end of Sec. 2.2.1). Regarding the gauge parameters of the 3d EFT, since the thermal nature of the 4d theory manifests itself only in the matching parameters, computations within the 3d EFT are carried out in its vacuum, simplifying matters greatly. This allows a strict -expansion, as long as the tree-level potential of the 3d EFT qualitatively agrees with the full effective potential. For the specific 3d EFT we study this is indeed the case, since there is a (temperature-dependent) tree-level barrier between the phases. Therefore, Nielsen identities are recycled in the -expansion, as outlined in Sec. 2.2.2, ensuring independence of the gauge fixing parameters ξ i,3 of the EFT, order-by-order.
In a gauge-invariant analysis, the gauge parameters can take any value, as dependence on them cancels exactly. However, in a gauge-dependent analysis this is no longer true. Sufficiently large values of |ξ i | scale as inverse powers of the coupling constants and violate perturbativity [127,128,134]. Since generic gauge-dependent loop corrections take the form ∼ g 2 ξ i at zero temperature, this reduces to ∼ gξ i at high temperature (cf. Eq. (2.31)), where g is some dimensionless coupling. 17 Thus, perturbativity at high temperature constrains the gauge parameters by, dropping numerical factors. Hence, by varying the gauge parameter in a range around O(1), we can estimate the magnitude of the uncertainty in our results due to gauge dependence. A calculation of when perturbativity breaks down in general covariant gauges in the SU (2)-Higgs theory indeed demonstrated that it does so for ξ i ∼ O(1) [127]. The focus on small gauge parameters is also supported by various observations in specific models, in which gaugeinvariant analyses appear to agree well with the Landau gauge [77,131,135]. Ideally, we would like to compare the relative magnitudes of the uncertainty due to gauge dependence and renormalisation scale dependence, considered in Sec. 3.1. Both arise in perturbative calculations multiplied by coupling constants in essentially the same way, though the latter in logarithms. So for a relatively fair comparison, we vary the gauge parameters over the range ξ i = {0, 3}, which corresponds approximately with ln(μ max /μ min ) = ln(4π) ≈ 2.53.
Finally, for values of ξ i as large as 10, the imaginary part of the 4d effective potential exceeds the real part at temperatures close to the critical temperature. Consequently, the 4d approach breaks down completely there due to the argument of the square root in the Goldstone modes (Eqs. (A.39)-(A.40)) which grows with ξ i . Our calculation of the critical temperature may therefore underestimate the theoretical error, accounting for the behaviour in the top left panel of Fig. 4. Specifically, for large values of ξ i an extra minimum forms in the real part of the potential, and the interplay of the two minima accounts for the behaviour of T c . Fortunately, the imaginary part of the effective potential decreases with temperature and tends to be relatively small near the percolation temperature.
Thus we may conclude that while in the 4d approach gauge dependence represents a limiting theoretical uncertainty in the gravitational wave spectrum of the phase transition, it is conceptually absent in the 3d approach.

High temperature approximation
Both daisy resummation -à la Arnold-Espinosa -and dimensional reduction at least implicitly rely on the high-temperature approximation, as they are predicated upon a hierarchy of energy scales. The hard thermal scale ∼ πT is assumed much larger than the masses of bosonic zero modes. This assumption generically holds for thermally driven phase transitions near the critical temperature. This is ensured by the structure of the loop expansion near the critical temperature; see Eq. (2.32). However, for transitions dominated by vacuum (rather than thermal) physics, in which there is a tree-level barrier between phases at T = 0, a lot of supercooling can occur between the critical temperature T c and the percolation temperature T p . In this case the high-temperature approximation can break down at T p .  Figure 4. Gauge dependence of thermal parameters for the 4d approach atμ = T (black) and the 3d approach atμ = 2.2T (blue). In both cases the continuous lines denote ξ 1 = ξ 2 = 0 and the dot-dashed lines denote ξ 1 = ξ 2 = 3. At ξ 1 = ξ 2 = 10 the 4d approach breaks down, whereas in the 3d approach the artificial gauge dependence is largely indiscernible, even to very large values of ξ i . A dashed blue line demonstrates this at ξ 1 = ξ 2 = 100 for the 3d approach. Note that the residual gauge dependence in the 3d approach is an artifact of an incomplete operator basis in this EFT, and as such is not morally equivalent to the inherent gauge dependence in the 4d approach.
In the 4d approach, the high-temperature approximation enters explicitly in the thermal (Debye) masses, Sec. A.3, but also implicitly through the singling out of the bosonic n = 0 Matsubara modes for resummation. Hence our daisy resummation relies on the hightemperature approximation even though we numerically evaluate the full m/T dependence of the thermal functions in the effective potential, Eqs. (A.68) and (A.69). Alternative approaches which do not rely on the high-temperature approximation have been developed in e.g. Refs. [51,54,98].
While it is not straightforward to quantify the accuracy of the high-temperature approximation in the 4d approach, we can estimate its effect by the size of the error introduced by approximating a choiceμ = T for the RG-scale -are plotted in Fig. 5 as the dotted line. The discrepancy with the full line in that figure indicates our estimate for the size of the uncertainty in the 4d approach which stems from the high-temperature approximation. In this estimate, we do not also expand the fermionic thermal functions, as only the bosonic degrees of freedom are resummed, so the high-temperature approximation does not enter the fermionic sector. We find a difference in the overall gravitational wave spectrum of order ∆Ω/Ω = O(10 −1 − 10 0 ) with the larger uncertainties at smaller M . Within the 3d approach of dimensional reduction, the high-temperature approximation is closely related to the truncation of the 3d effective theory. The coefficients of higher dimensional operators in the 3d effective theory are related to the coefficients in the high-temperature expansion of the hard mode parts of thermal functions. Thus, large contributions from the addition of higher dimensional operators to the 3d effective theory signals the breakdown of the high-temperature approximation.  Figure 6. Estimate of the uncertainty in our 3d approach due to truncating the effective theory at (φ † φ) 3 . The cyan lines are the LO results in this truncation, the red lines include also the (φ † φ) 4 operator and the orange lines include both the (φ † φ) 4 and (φ † φ) 5 operators.
To estimate the leading corrections to our truncation of the high-temperature expansion, we match the scalar dimension-8 and -10 operators, (3.8) These operators enter the Higgs potential directly, and hence we can estimate the order of magnitude of their effects by analysing the tree-level potential in 3d with these operators. By a direct extension of the tree-level analysis in Sec. 2.2.2 and Appendix B.5, we calculate T c , T p , α, and β/H p including the operators in Eq. (3.8). The difference between this and the tree-level result without operators O 8 and O 10 estimates our uncertainty. Fig. 6 shows this uncertainty estimate for the four thermodynamic parameters. The effect of c 8,3 is relatively small, and the additional effect of c 10,3 even smaller, suggesting good convergence of the expansion. Note that the effect of these higher dimensional operators is seen to be larger for stronger transitions. Since the discontinuities of the scalar condensates are larger for stronger transitions, in general (φ † φ) n -operators give larger effects.
Higher dimensional operators in the 3d EFT do not just arise from the high-T expansion, but also arise necessarily at higher loop orders. Powers of m/T from the high-T expansion compete with powers of the coupling constants which arise from the loop expansion [149]. Thus, in a rigorous power-counting scheme, as one increases the accuracy of the calculation, and includes more loop orders, one will also need to match to correspondingly higher dimensional operators [59]. This is also necessary in order to render soft and ultrasoft observables finite [107].
While the high-T expansion of thermal integrals is utilised in our matching relations for the 3d approach and simplifies the matching significantly, it can be avoided at the cost of tougher master integrals; see Refs. [126,150]. Rather than attempting this, we utilise the 4d approach to estimate the error introduced in the 3d approach from the use of the expansion of thermal integrals in the matching relations. To do this we calculate the change in thermal parameters as calculated in the 4d approach when all bosonic and fermionic thermal functions are expanded (up to and including the logarithm); see Eqs.

Higher loop orders
At zero temperature, the convergence of the loop expansion is dictated by the smallness of the dimensionless coupling constants. Large zero-temperature couplings will correspond to large theoretical uncertainties. Being interested in studying finite temperatures, the convergence of the loop expansion is more delicate due to the Infrared Problem.
In the 4d approach, our calculations reach one-loop level. There is unfortunately no treelevel result to compare with our one-loop calculations, as the phase transition takes place due to thermal fluctuations, which appear first at one-loop. Already Refs. [51,52] extended the daisy-resummed approach to two-loops, though realistic Standard Model extensions were only recently tackled for the two-Higgs doublet model [54] and the SM with a lighter Higgs mass [76]. In the daisy-resummed approach, the extension to two-loops is complicated due to massive two-loop sum-integrals, as well as the large number of different particles contributing As such, we refrain from two-loop corrections in the 4d approach.
In the 3d approach, next-to-leading order (NLO) calculations become more amenable. This is because, for NLO matching, the only two-loop corrections required are the thermal mass corrections, and thermodynamic properties are analysed within the simpler 3d effective theory. This approach expands separately: (i) the matching relations in powers of g 2 and (ii) the 3d effective theory in powers of . In order to calculate some physical quantity to a given order in g 2 , both expansions should reach that same order. It is nevertheless instructive to test their convergence separately.
For the matching relations in the 3d approach, we compare three different approximations: LO matching, one-loop matching and NLO matching; see Fig. 7. The LO and NLO approximations are both consistent truncations of the perturbative series in powers of g 2 , and as such are gauge invariant 18 . The LO approximation consists of one-loop matching of masses, and tree-level matching of couplings. The NLO approximation consists of two-loop matching of masses, and one-loop matching of couplings, and this is what we utilise elsewhere in the paper. As argued in Ref. [63], the LO approximation is not expected to be quantitatively accurate as fractional corrections to observables from NLO corrections are O(1). One must carry out NLO matching for higher order corrections to observables to be perturbatively suppressed. In between the LO and NLO approximations lies the full one-loop approximation. This constitutes an incomplete calculation at any order in g 2 , and consequently is not gauge invariant (e.g. Eq. (B.37)), just like the 4d approach. For example, similarly to the 4d approach discussed in Sec. 2.2.1, a full one-loop matching includes some O(g 4 ) corrections to the effective potential, but misses others. However it is simpler than the full NLO matching and includes logarithms which cancel the running of couplings (but not masses); see Refs. [96,115,151,152].
The results for the SMEFT shown in Fig. 7 inspire some confidence in the NLO matching relations which we have utilised. For T c and T p , as one progresses from lower to higher order approximations, from LO to 1-loop and then NLO, the scale dependence bands shrink from ∼ 20 − 50% to ∼ 10% to ∼ 2%, while the centre of the bands, at the optimal scalē µ = 2.2T , change relatively little. This suggests good convergence as one increases the order of the approximation for the matching relations. It also shows the importance of higher order corrections for reducing unphysical scale dependence. For α and β/H, while the results at the optimal scale show good convergence, the LO approximation appears to break down at small M andμ = T /2. Nevertheless, that the NLO scale dependence bands lie entirely within the 1-loop bands still suggests relatively good convergence. We can naively estimate the size of the unknown NNLO corrections to T c as (1 − LO/NLO) 2 ∼ (10%) 2 ∼ 1%, with all quantities evaluated at the optimal scale. For T p this estimate for the NNLO corrections is ∼ 1 − 4%, for β/H it is ∼ 4 − 15% and for α it is ∼ 10 − 50%. These values naively estimate NNLO corrections to the matching relations to result in ∆Ω/Ω = O(10 0 − 10 1 ), with the larger uncertainties at smaller M . This estimate is supported by its agreement with the magnitude of the scale uncertainty of the NLO result.
In the 3d approach, we also assess the convergence of the -expansion within the 3d EFT. For purely equilibrium quantities, such as the transition strength evaluated at the critical temperature, α c , we can -expand up to O( 2 ), or NNLO. This utilises the two-loop effective potential computed in Appendix B.4. However, for the bubble nucleation rate, the spatial dependence of the bubble profile severely complicates the computation of higher loop orders, so we are only able to compute the LO in and estimate the NLO in . Due to this distinction, we defer the discussion of the nucleation rate to Sec. 3.5.
The presence of the single higher dimensional operator in the truncation of the 3d effective theory leads to an interesting structure of the perturbative series. Essentially, if c 6,3 is the smallest coupling in the EFT, it determines the convergence of the entire loop expansion Here the matching relations are fixed at NLO. The convergence of the loop expansion within the 3d EFT is fairly consistent across the considered range of M (cf. Fig. 8). From Eqs. (3.9) and (3.10) we would naively expect three loop (N 3 LO) and higher loop corrections in this expansion to give fractional corrections of order (1−LO/NLO) 3 ∼ (10%) 3 ∼ 0.1%. However, it should be remembered that, in the symmetric phase the ultrasoft spatial gauge bosons are nonperturbative. We discuss this issue in Sec. 3.6.

Nucleation corrections
The nucleation rate is by far the most technically challenging quantity for which to calculate higher loop corrections. In fact, even calculating the leading order self-consistently is nontrivial, something which is often unappreciated. It is the intersection of the following two points which lead to this technical difficulty: (i) The phase transition occurs due to thermal fluctuations. These appear first at one loop, hence the critical bubble (or bounce) cannot be solved for at tree-level. 19 (ii) The fields should be loop-expanded around the critical bubble, i.e. around an inhomogeneous classical background field, φ = φ(x).
These two points lead to a catch-22: the critical bubble only exists in the background of one-loop corrections and yet the one-loop corrections should be made around the background of the critical bubble. This, coupled with the inhomogeneity of the critical bubble, constitute the main technical challenges in consistently calculating the nucleation rate. An intuitive solution to (i) above is to use the perturbatively computed effective potential, rather than the tree-level potential, to compute the critical bubble [87,121]. Starting from the tree-level, vacuum bounce equation of Coleman [139], the suggested prescription to modify it for the thermal case is the following: This prescription gives exactly Eq. (2.28) once assuming O(3) symmetry and forms the basis of our nucleation calculations in the 4d approach, as discussed in Sec. 2.1. It is also the most common approach taken in the literature. While plausible, and clearly a step in the right direction, the naive replacement of Eq. (3.11) is not derived from first principles. It suffers from inconsistencies since, in general, the bounce action thus calculated is not the correct result at leading order in any expansion, as has have been discussed in Refs. [129,130,[153][154][155][156][157][158][159][160][161]. 20 The presence of the extraneous nonzero imaginary part of the effective potential, which must be discarded by hand, perhaps offers a clue that something has gone wrong.
In Eq. (3.11), by using the effective potential in the bounce equation as an attempt to solve (i), one evades the catch-22 by inconsistently integrating over the fields twice. Integrating first in generating the effective potential and then second in solving the bounce equation and evaluating the nucleation prefactor. In computing the effective potential all nonzero momentum modes of the fields are integrated out, including those of the scalar. The remaining degree of freedom, the constant mode of the scalar field, cannot describe a localised, inhomogeneous critical bubble (or bounce), hence the practical necessity for ad hoc promoting the constant mode back up to a full spatially-dependent field, to be integrated over a second time. Furthermore, in computing the effective potential, the spatial dependence of the scalar field, the basis of (ii), is still wrongly ignored. One might think that a derivative expansion could justify this, and that the effective action evaluated on a background bounce solution might be approximated by the naive bounce action, where we have indicated the size of the simplest derivative corrections, denoting by . . . all other possible terms. Here m φ is the mass of the φ particle, i runs over the modes with masses m i which have been integrated out in deriving V eff , n runs over 2, 3, 4, . . . and the C i,n are O(1) constants. The solution to Eq. (3.11), will exhibit a virial-type theorem implying that the bubble wall has a width of order ∼ 1/m φ , so that ∇/m φ ∼ 1 when evaluated on the bubble wall. By deriving the effective potential the φ field itself has been integrated out, so that m φ ∈ {m i }. Therefore, even ignoring the inconsistencies of double counting, the derivative expansion in this case is, at best, an expansion in powers of ∼ m φ /m φ = 1, and even worse if there are any bosons lighter than φ. 20 See also Ref. [162] for the resolution of an analogous issue at zero temperature.
A consistent resolution to these issues was given by Langer [89,153] in the context of classical statistical mechanics, which has been formulated in quantum field theoretic terms in Refs. [129,130,[156][157][158]160]. 21 The resolution depends upon the existence of a certain hierarchy of scales, with UV and IR modes well separated in energy scales, i.e. Λ IR ≪ Λ UV . In essence, one first integrates over UV modes, resulting in a temperature-dependent, coursegrained effective action, Γ IR , for the remaining IR modes. One can then derive a bounce equation for the IR modes, which stationarises this course-grained effective action. This resolves (i) without double counting. Further, as only the UV modes have been integrated out, the course-grained effective action is only nonlocal on length scales 1/Λ UV ≪ 1/Λ IR . Hence a derivative expansion of Γ IR is applicable, and amounts to an expansion in powers of Λ IR /Λ UV . This consistently resolves (ii). The action from which one determines the bounce equation is where φ IR are the IR modes and V IR and Z IR are the potential and field renormalisation for φ IR , after having integrated out only the UV modes. Again we have indicated the size of the simplest derivative corrections, with C n being O(1) Wilson coefficients. In this case the bubble wall width is of order 1/Λ IR and hence these higher order derivative corrections are suppressed by powers of Λ IR /Λ UV . These terms, as well as those omitted as . . . in Eq. (3.13), can thus be neglected when there is a sufficiently large separation of scales. Dimensional reduction is built around just such a separation of scales, summarised in Table 2. It therefore provides a natural framework for which to perform a self-consistent calculation of the bubble nucleation rate. After integrating out the modes on length scales Λ UV ∼ 1/πT and then ∼ 1/gT , we are left with an effective theory on length scales Λ IR ∼ π/g 2 T , the ultrasoft theory. Note that the explicit Z IR factors are removed by a redefinition of the infrared fields according to Eq. (B.56). Thus, by identifying φ IR with φ 3d of the 3d EFT and Γ IR with the tree-level action of the 3d EFT, we can self-consistently calculate the tunnelling action using Eq. (3.13), obtaining the correct result at LO in the 3d -expansion (see Eq. (2.50)), and to NLO in powers of Λ IR /Λ UV ∼ g/π. Terms with more powers of φ IR or more derivatives are suppressed by powers of the coupling, as long as the hierarchy of scales in Table 2 holds.
The resulting tunnelling action is independent of the gauge fixing within the 3d EFT, being the leading order in a consistent -expansion. This extends also to the nucleation prefactor, which is gauge invariant when evaluated on a solution to the tree-level equations of motion [137,138]. Therefore, if the matching relations are independent of the gauge fixing within the 4d theory (as, for example, we have shown them to be up to O(g 4 ) in the SMEFT), the calculation is then gauge invariant from end to end. Conversely the small gauge dependence introduced into the matching relations in the SMEFT due to the incomplete basis of operators will carry through to the tunnelling calculation.
The picture formed is that the hard, UV modes cause the 3d effective parameters to run with temperature, driving the light, IR modes through the transition. Note that it would be incorrect to equate V IR with the effective potential in the 3d effective theory, after having integrated out the ultrasoft fields, as that would amount to double-counting and an uncontrolled derivative expansion akin to Eq. (3.12). The correct identification for V IR is the tree-level potential of the ultrasoft EFT.
Returning now to the naive recipe, Eqs. (3.11) and (3.12): despite its inconsistencies, in certain circumstances it approximates Eq. (3.13). This is the case when the scalar undergoing nucleation is much lighter than all other particles in the theory and has much smaller selfcouplings than couplings to other fields. This occurs, for example, in the Standard Model in the region of parameter space where the Higgs is much lighter than the W bosons. In this case the derivative expansion of the effective action is justified for all diagrams except those containing scalar loops, however scalar loops are subdominant by assumption. Further, it is the scalar loops which lead to the erroneous imaginary part of the effective potential, and this again will be subdominant. Nevertheless, Refs. [129,130,133,155,157] have pointed out that rather than accepting a subdominant inconsistency in the calculation, one should instead omit the scalar loops in the potential used in the bounce calculation.
The NLO correction in to the bubble nucleation rate is the nucleation prefactor, the term A in Eq. (2.1). Its most important contribution is essentially the entropy of the critical bubble, given by integrations over all field fluctuations in the vicinity of the critical bubble. At this point the inconsistency of Eq. (3.11) becomes very clear, as all field fluctuations have already been integrated over in computing V eff , mixing the LO and NLO terms in . Hence it is not possible even to define the nucleation prefactor in our 4d approach, resulting in a fundamental theoretical uncertainty beyond just the breakdown of the derivative expansion implied by Eq. (3.12). As a consequence, most of our analysis of the nucleation prefactor is carried out in the 3d approach, though through this we are also able shed some light on theoretical uncertainties in the 4d approach.
The nucleation prefactor can be split into a product of two terms, called the dynamical and statistical prefactor respectively [164], (3.14) Both of these terms are independently technically challenging to compute from first principles, and we refrain therefrom in this paper. The dynamical prefactor requires a real-time computation in the presence of a thermal bath, while the statistical prefactor requires computing functional determinants. To estimate the magnitude of corrections related to the nucleation prefactor we compare various approximations present in the literature. This is discussed at length in Appendix C.
We consider two different approximations to the dynamical prefactor, A dyn . These are (A) an approximation in which the ultrasoft gauge bosons are assumed to dominate the dynamics and (B) a hydrodynamic approximation, which also assumes the bubbles are thin walled. Our default choice in the 3d approach is (A), the infrared gauge boson dominance approximation.
We also consider three different approximations to the statistical prefactor, A stat . For the statistical prefactor, the approximations are (a) a thick wall approximation, correct up to a multiplicative function of g ′2 3 /g 2 3 and λ 3 /g 2 3 , (b) a thin wall approximation, which correctly gives the dominant terms in the thin wall limit and (c) LO in an uncontrolled derivative expansion (comparable to Eq. (3.12)) with ad hoc treatment of the zero modes. This last approximation, (c), gives exponentially large errors in general but is expected to become more reliable in the limit that the Higgs is much lighter and more weakly coupled than the gauge bosons, and in the thin wall limit. Our default choice in the 3d approach is (a), the thick wall approximation.
The tunnelling action at nucleation is plotted in Fig. 9, together with prefactor corrections in our different approximations. These should be small corrections if the -expansion is under control. Further, agreement between different approximations to the dynamical and statistical parts of the prefactor respectively would suggest that theoretical uncertainties in these quantities are relatively under control. We remind the reader that the full prefactor is the product of the dynamical (orange) and statistical (cyan) parts.
The most obvious anomaly in Fig. 9, is the statistical prefactor in the uncontrolled derivative approximation, (c), which shows O(1) multiplicative deviations from the LO tunnelling action. This uncontrolled derivative approximation underlies the 4d approach, being essentially the same as that carried out in Eq. (3.12). The numerical results shown in Fig. 9 suggest the approximation breaks down badly, especially for larger M , where the transitions are weaker. Although the 4d approach relies on this derivative expansion, by incorporating these corrections directly in Eq. (3.11), rather than in the prefactor, it effectively resums the expansion, so mitigating (but not avoiding) the breakdown of the expansion. The consequences of this uncontrolled error cannot easily be quantified in ∆Ω/Ω min , so we denote them as unknown in Table 3. However, the systematic difference for β/H * with respect to the 3d approach may give some indication (see e.g. Fig. 2). Just this amounts to a systematic error ∆Ω/Ω min of around O(10 1 ).
For the 4d approach, leaving aside the unknown effect due to the application of the uncontrolled derivative expansion, the different approximations for the prefactor A result in a modest change in the effective action compared to the value recommended in the recent LISA review, log(A/T 4 ) ∼ −14 [6]. This in turn results in a modest change in the peak gravitational wave amplitude between O(10 −1 − 10 0 ). However, since the true error is dominated by the systematic error, we do not include this uncertainty in the final analysis.
On the other hand, for the 3d approach, we face no such difficulties and so can give a reasonable estimate for the theoretical uncertainty related to the nucleation prefactor. Fig. 9 shows that the thin and thick wall approximations to the statistical prefactor roughly agree, as do the two approximations of the dynamical prefactor. Better agreement than this is not The full prefactor factorises into dynamical (orange) and statistical (cyan) parts; see Eq. (3.14). More details of these approximations can be found in Appendix C. Note the large deviations from the LO result for the derivative approximation to the statistical prefactor. This aligns with the expected breakdown of this approximation, which underlies the 4d approach.
to be expected, as these are all crude estimates. However, their rough agreement encourages us to suggest that they could be used to estimate the magnitude of O( ) corrections in the 3d approach. Combining the various approximations to the dynamical and statistical prefactors (excluding the uncontrolled derivative expansion) in all possible ways leads to a range of different estimates. These different estimates amount to an uncertainty ∆Ω/Ω min for the 3d approach, which varies monotonically from O(10 −1 ) at M = 700 GeV to O(10 0 ) at M = 580 GeV.

Nonperturbativity
In the symmetric phase the perturbative mass of the magnetic gauge bosons is zero. As a consequence, the effective coupling of their zero Matsubara modes, ∼ g 2 T /m, appears to be infinite. This leads to IR divergences at finite loop order [49], and the consequent complete breakdown of perturbation theory. 22 Physically, these divergences in perturbation theory are softened by a nonpertubatively generated mass for the gauge bosons, of order g 2 T [9]. Based on power-counting arguments, one would expect the resulting softened IR divergences to contribute to the free energy density at O(g 6 T 4 /(2π) 4 ) [165,166], or the typical size of a 4loop order term. However, this estimate is perhaps misleading, as the effective gauge coupling becomes large for momentum-transfers Q 2 ∼ (g 2 T /π) 2 , so perturbative estimates cannot be relied upon. The nonperturbative nature of the symmetric phase leads to an irreducible uncertainty in our perturbative calculations which in principle only lattice simulations can resolve.
In the broken phase, by contrast, perturbation theory may work well, as the troublesome magnetic gauge bosons acquire a mass via the Higgs mechanism. For a sufficiently strong transition, this broken phase mass is large, and hence the effective coupling of the magnetic gauge bosons, ∼ g 2 T /m, is perturbative. Of course, both phases are relevant for calculating the thermodynamic properties of the transition. However, if for some thermodynamic quantity the contribution of the broken phase exceeds the contribution of the symmetric phase, then one would expect perturbation theory to be an at least qualitatively reliable guide. This scenario should take place for sufficiently strong transitions, in which the Higgs vev in the broken phase is large, and gives correspondingly large contributions to thermodynamic quantities.
Explicit comparisons between lattice and perturbation theory have indeed shown good qualitative agreement for strong first-order transitions in the xSM, 2HDM, ΣSM [24,83,118], Abelian Higgs model [131] and also in the SM [9,11,133] (with artificially light Higgs). In particular, the scan of the parameter space of the xSM in Ref. [83] showed good qualitative agreement between lattice and perturbative calculations for the phase diagram of the theory. For the two benchmark points in the 2HDM considered in Ref. [24] it was found that T c calculated in a dimensionally-reduced approach 23 differed by 4-7% from the lattice result, whereas α c differed by 5-25%. By contrast, calculations using a daisy-resummed approach showed discrepancies from the lattice result for T c of 20-45% and for α c of 45-75%. For the two benchmark points considered in the ΣSM in Ref. [118] it was found that in a dimensionallyreduced approach α c differed from the lattice result by 30-40%, whereas in a daisy-resummed approach the discrepancy was more than 50%. Bubble nucleation was studied on the lattice in the Standard Model with light Higgs in Ref. [133]. There it was found that perturbative calculations of the value of (m 2 3,c − m 2 perturbative calculations are significantly larger than one might guess based on perturbative power counting, highlighting the importance of carrying out nonperturbative computations.

Discussion
This work systematically examines theoretical sources of uncertainty in the finite-temperature calculation of thermal tunneling rates, and the resulting uncertainty in thermal parameters of a first-order phase transition. In particular, we have compared the sources of uncertainty in two different methods that address the breakdown of perturbation theory due to the longwavelength modes at high temperature: resummation of daisy diagrams in 4d, and dimensional reduction to 3d effective theory. The benchmark model we have used is the Standard Model aided with a dimension-six operator, O 6 = φ † φ 3 /M 2 , which we expect represents qualitatively a large set of models of EWSB, which is of particular interest for the planned LISA experiment. Section 2 and the Appendices include a thorough and comprehensive review of the 4d and the 3d approaches, and the calculation of the thermal parameters capturing the dynamics of the phase transition. We develop the 3d approach, building upon previous works utilising dimensional reduction, to include also a gauge invariant treatment up to O(g 4 ) of bubble nucleation and the calculation of thermodynamic quantities. The two approaches are outlined in general, and are fleshed out for the specific example of our benchmark model.
The main sources of uncertainty were categorised in Section 3, as follows: renormalisation scale dependence, gauge dependence, the high temperature approximation, the unknown contributions of higher loop orders, corrections to the bubble nucleation rate, and nonperturbativity or the Infrared Problem.
We include a qualitative discussion of the origin of these uncertainties, and where possible explicitly calculate their magnitudes over the relevant range in M in our benchmark model. We summarise the results of these calculations in terms of the parameter ∆Ω/Ω min defined in (3.1) in Fig. 10. Importantly, the extent of the uncertainty calculated in this paper can not be taken as a direct predictor of the correctness of the result, as is most obviously noticeable from the apparently vanishing uncertainty for certain M in some of our figures. However, we expect that the variation of the results may be incorporated as an integral part of the diagnosis of the appropriateness of the two methods.
From the analysis in this paper, we may draw several conclusions. Firstly, a direct comparison between the 4d and 3d approaches indicates that the latter generically implies a significantly smaller theoretical uncertainty. 24 Such a comparison is shown in Fig. 10, and yields a difference in ∆Ω/Ω min of several orders of magnitude. This may be taken as support for the choice of the 3d approach over the more commonly used 4d approach. Secondly, of those uncertainties which we were able to estimate, the most important driver of theoretical 24 As the calculations in both approaches here were done semi-analytically and did not utilise lattice Monte-Carlo simulations, there remains a nonperturbative uncertainty associated with the magnetic gauge bosons in the symmetric phase; see Sec. 3.6. gives the difference between the result including that parameter and the result with only c 6,3 . The purple line labeled "high-T " gives the variation due to the high-temperature expansion of the thermal functions J b/f . Note that theoretical uncertainties due to nonperturbativity and due to the breakdown of the derivative expansion in the 4d approach are not included in this plot, as we were unable to give reliable and detailed estimates of the magnitude of these uncertainties.
uncertainty in the 4d approach is the dependence on the choice of renormalisation scale. This effect dwarfs the uncertainty introduced by other effects, such as gauge dependence -as uncomfortable as the latter is. In light of the results presented here, we recommend forthcoming works to include an analysis of the variation of the renormalisation scale. Thirdly, the error introduced by applying the high-T approximation to thermal functions grows sharply for the strongest transitions, below about M = 590 GeV. This was traced back to the contribution of the top quark. As these strong transitions are those most relevant to gravitational wave experiments such as LISA, we suggest that at least the full, numerical m t /T dependence should be accounted for in future calculations of the SGWB of first-order phase transitions. In addition to those theoretical uncertainties for which we were able to give reliable estimates, we demonstrated that the 4d approach generically depends upon an uncontrolled derivative expansion. The breakdown of this expansion is remedied in the 3d approach, and hence may account for the systematic discrepancy in β/H * between the two approaches. We advise that care should be taken in future studies to avoid this particular stumbling block and to ensure that calculations of the tunnelling rate are self-consistent.
Prior to this work, there have been many studies of various different theoretical uncertainties in calculations of the thermodynamics of phase transitions; recent examples include Refs. [4,20,24,54,55,[76][77][78][79][80][81][82][83][84]. However, for the first time, in this work we have comprehensively analysed and compared a wide variety of relevant theoretical uncertainties across the full range from weakly to strongly first-order transitions. In doing so, we have focused on the implications for the resulting gravitational wave spectra. Previous works have almost exclusively focused on the theoretical uncertainties of the transition temperature and the corresponding vev, which can disguise the severity of the theoretical uncertainty in the gravitational wave spectrum. 25 For example, as α is inversely proportional to the fourth power of the percolation temperature, and the peak gravitational wave amplitude in turn depends quadratically on α, an apparently innocuous uncertainty ∆T p /T p = 0.1 in the percolation temperature will result in an uncertainty ∆Ω/Ω ≈ 1 in the gravitational wave spectrum. In fact this argument may significantly underestimate the uncertainty, as additionally the trace anomaly grows strongly as the percolation temperature decreases. To the extent that a comparison is possible, the results presented in this paper are qualitatively consistent with the existing literature.
The results in this paper will be of particular interest for gravitational wave studies in anticipation of the LISA experiment and other space-based interferometer experiments with sensitivity in the milli-Hertz range. In the benchmark model studied in this paper, the strongest phase transitions -leading to the gravitational wave spectra most likely to be observable -are found for smaller M . Importantly, this is also the range in which the diagnosis presented in this paper indicates the poorest theoretical control, cf. Fig. 10. This implies that the amplitude of the gravitational wave spectrum, but also its peak frequency, can only be predicted to a level of accuracy which depends on the method used. This has serious implications for the prospect of model differentiation and complementary studies with collider probes. The sheer magnitude of the theoretical uncertainty and the comparative success of the 3d approach at NLO, motivates its use in studying GW phenomenology in specific BSM models. Further, we advocate that the theoretical uncertainty should be taken seriously and analyzed in all future phenomenological gravitational wave studies. gravitational wave spectrum within the Abelian Higgs model. See also Ref. [84] in which the gravitational wave spectrum of the Z2-symmetric xSM was studied, and within which a ∆Ω/Ω of O(10 4 ) was observed in the comparison of different approximation schemes.

A SMEFT in four dimensions
This appendix collects multiple technical details of our computation in the four-dimensional SMEFT, as defined in Eq. (2.7). The conventions and notations of the SM parts follow Ref. [96]; see Sec. 2 of that reference.

A.1 Renormalisation: counterterms and running
Excluding the example in Fig. 1, all calculations in this paper include RGE running of the parameters. While our main focus is the SMEFT with only the inclusion of the sextic Higgs operator c 6 (φ † φ) 3 , below we do include similar operators of dimension-8 and -10 with coefficients c 8 and c 10 . For example, note that below the counterterm for the dimension-8 operator -which is non-zero even if the coefficient c 8 itself vanishes -is crucial in order to cancel divergences related to the dimension-6 coefficient c 6 . This is due to the non-renormalisable nature of the SMEFT, where new higher dimensional operators have to be included to cancel divergences related to operators of lower dimension. For the same reason, the running of higher dimensional operators is non-zero, so even if these operators vanish at some initial scale, with RGE running they are non-zero at other scales. For this reason, one has to ensure that this running does not ruin the numerical analysis in the EFT without these operators. However, in most of our numerical analysis we simply neglect dimension-8 and -10 operators altogether.
For renormalisation, we use dimensional regularisation in the MS-scheme, and for gauge fixing we adopt the class of general covariant (or Fermi) gauges, introducing gauge-fixing parameters ξ 1 for the U(1) and ξ 2 for the SU(2) sector. In Landau gauge ξ i = 0. The one-loop counterterms that are affected by c 6 , c 8 and c 10 read The corresponding β-functions read Counterterms and β-functions unaffected by c 6 , c 8 and c 10 are collected for example in Sec. 3.2 of Ref. [96]. Note that since the presence of a nonzero-temperature preserves the UV structure of the theory, these same counterterms and β-functions are used in the renormalisation of unbroken phase correlators for dimensional reduction and for the broken phase vacuum renormalisation calculation of pole mass corrections.

A.2 Relations between MS-parameters and physical observables
We relate the MS-parameters of the Lagrangian to physical observables, that serve as input parameters Note that the physically observed masses are the pole masses. These relations also depend on the new MS-scheme BSM parameters c 6 , c 8 and c 10 , which we also treat as input parameters.
For the values of the physical observables used in this work, we refer the reader to Table 1. We define the shorthand notation g 2 0 ≡ 4 √ 2 G f M 2 W for the tree-level coupling and v 2 0 ≡ 4M 2 W /g 2 0 ≈ (246.22 GeV) 2 for the tree-level minimum. At tree-and one-loop level only the Higgs mass parameter and self-coupling are affected by c 6 , and the tree-level relations can be solved from (V tree is defined in Eq. (2.9)) resulting in At tree-level, the relations for gauge and Yukawa couplings are unaffected by c 6 and read For an accurate numerical analysis of the thermodynamics, the above tree-level relations can be improved by their one-loop corrections (cf. Refs. [24,54,63]). These corrections are necessary for the complete O(g 4 ) accuracy of our 3d approach. Regarding the masses, this can be achieved with a standard one-loop pole mass renormalisation at zero temperature. With the Minkowski metric at zero temperature, propagators are schematically dressed as 1/ p 2 − m 2 + Π(p 2 ,μ) , where Π is a self-energy function with external momentum p,μ is the MS-scale, and m 2 is the MS-mass eigenvalue (see Eqs. (A.38)-(A.44)) at the tree-level minimum v 0 . Diagrammatically Π consists of The self-energy functions capture the momenta dependence of the two-point correlators, and in addition include one-particle reducible tadpole contributions. Those are generated at oneloop by the non-zero vev in the broken phase computation, since the minimisation condition is imposed only at tree-level. Identifying the physical pole mass at p 2 = M 2 leads to a condition m 2 = M 2 + Π(M 2 ). We denote physical masses with capital letters, and MS-mass eigenvalues -that are functions of running couplings -by lower case letters. Note that self-energies Π are functions of MS-masses m 2 , but one may linearise these pole equations by replacing MSmasses by the corresponding physical pole masses inside the one-loop computation, since the difference is formally of higher order. 26 We employ this linearisation, as it suffices to reach the desired O(g 4 ) accuracy. The one-loop correction δg 2 to the SU(2) gauge coupling can be obtained by computing the one-loop 4-fermion correlator related to muon decay µ − → e − +ν e + ν µ at zero external momenta. This correlator is proportional to Diagrammatically, there is a tree-level contribution where the W -boson self-energy appears via the one-loop improved propagator, and the following classes of one-loop diagrams: where blobs denote possible one-loop attachments. Note that all c 6 -dependent contributions cancel between the W self-energy piece in the above tree-level diagram of Eq. (A.21) and the first diagram class of Eq. (A. 22), such that the final result is equal to the SM. Furthermore, this correction is numerically small. We adopt the approximation of Ref. [63], where the tree-level value of the U(1) gauge coupling is not improved at one-loop, due to its small numerical significance. Note that one could include one-loop corrections to the U(1) gauge coupling for example by computing the correction to Thomson scattering (cf. Ref. [24]).
In total, we have the conditions: Here the top quark self energy consists of a vector part Σ v and a scalar part Σ s ; see Sec. 5 in Ref. [63]. Its axial and axial vector parts do not contribute to the pole mass condition.
From these equations we can now solve for the one-loop improved MS-parameters in terms of the physical parameters: The final task is to evaluate the self-energy functions. This calculation is straightforward since there are no new fields compared to the pure SM computation: there are no new diagrams involved, and c 6 only modifies the mass eigenvalues and vertices. By direct computation -and by adopting the shorthand The one-loop integral function F (k, m 1 , m 2 ) and its various limits are given in Eqs. (187)(188) in Ref. [63]. Also note that in the above self-energies, U(1) gauge contributions can be turned off by taking the limit z → 1 and M Z → M W . Apart from the new c 6 terms and U(1) gauge contributions, these results agree with Eqs. (184), (191)(192)(193) in Ref. [63], apart from a − 8 3 t 2 ln(h) term which is missing in the top quark self energy in Eq. (193) therein. 27 From the above self-energy functions and relations between MS-parameters and physical parameters we observe that g and g Y are independent of c 6 also at one-loop order. Even though self-energies for both the W -boson and the top quark have c 6 -dependent pieces, these cancel each other, and the 4-fermion correlator of muon decay is explicitly independent of c 6 .
Note that often [27,74] one-loop improved relations for the scalar mass parameter and self coupling are obtained from the conditions where the tree-level potential is accompanied with the one-loop Coleman-Weinberg potential at zero temperature (see Eq. (A.67)), and µ 2 h and λ are solved numerically from these equations. However, these conditions should be taken only as a heuristic approximation, since the physical pole mass lies at nonzero momentum and hence cannot be obtained from the effective potential (see also Ref. [26]). Figure 11 compares the scalar parameters as functions of c 6 = 1/M 2 , at tree-level versus one-loop -computed both from the Coleman-Weinberg potential and from pole conditions. This comparison shows that in fact the estimated 1-loop effect from the Coleman-Weinberg potential overestimates the difference with respect to the tree-level result, compared to the computation from pole masses that consistently include momentumdependent contributions and tadpole diagrams.

A.3 Mass eigenvalues and thermal screening
Calculating the effective potential to one-loop level by using field-dependent mass eigenvalues, requires fixing a gauge. In this context it is common to use an R ξ -style gauge. However this fixes a different gauge for each field value; regarding the gauge-dependent effective potential, it is not clear a priori that this is permissible [127,167]. We therefore follow Ref. [168] by using a general covariant (or Fermi) gauge,which does not include the vev in the gauge-fixing Lagrangian, to calculate the loop corrections to the effective potential. Gauge parameters are denoted as ξ 2 and ξ 1 for the SU(2) and U(1) fields respectively. In this case the field-dependent mass eigenvalues are given by where m 2 χ = µ 2 h + λφ 2 + 3 4 c 6 φ 4 is the Goldstone mode mass eigenvalue in the Landau gauge. Note the complicated gauge dependence in the Goldstone mass eigenvalues m 2± 1 , m 2± 2 , and m 2± 3 . A very illuminating and thorough derivation for the gauge dependent parts of these mass eigenvalues can be found in Ref. [169] (see also Ref. [136]).
Turning to high temperatures, the most straightforward way to obtain thermally resummed mass eigenvalues, is to calculate thermal corrections from hard modes to the zero modes of the fields in the gauge eigenstate basis of the theory; we demonstrate this explicitly in Appendix B.2. Gauge fields A a 0 and B 0 obtain Debye masses The number of kinematically active families is N f = 3 at electroweak-scale temperatures. For temporal gauge fields, the bilinear part of the Lagrangian reads where thermal corrections read However, due to thermal screening this is not an eigenstate basis, as in fact there remains a mixing term 54) i.e. the Z-boson and photon are not mass eigenstates in the heat bath. The actual resummed mass eigenvalues read where the rotation angle to define eigenstates Z ′ , A ′ depends now on m 2 D and m ′2 D and exhibits a more complicated form compared to the usual vacuum Weinberg angle. Note that if these mass eigenvalues are linearised with respect to the thermal contributions m 2 D and m ′2 D , one gets exactly the above resummed Z-boson and photon masses. In fact this difference is numerically small, which suggests that one could still treat the Z-boson and the photon as mass eigenstates. Also note that in the above Π Z,W,γ are for the longitudinal modes (perturbative Debye masses for the transverse modes vanish).
For the scalar doublet, thermal screening affects both its mass and self-coupling, with contributions , where cos(θ) = g/ g 2 + g ′2 , sin(θ) = g ′ / g 2 + g ′2 . and the resummed parameters become The above c 6 contributions are those given by the flower diagrams in Eq. (2.25). Resummed scalar mass eigenvalues are then obtained by substituting the mass parameter µ 2 h and selfcoupling λ for their resummed values in m 2 φ and m 2 χ .

A.4 One-loop thermal effective potential
This appendix composes the thermal effective potential for the SMEFT in four dimensions and implements the leading ring resummations, utilising the mass eigenvalues and thermal corrections of Appendix A.3. As a reminder, we fix the MS-parameters on the Z-pole as explained in Appendix A.2 and evolve them using the renormalisation group equations (RGEs) in Appendix A.1.
The one-loop correction to the effective potential can be computed as a sum over all oneparticle irreducible diagrams with a single loop and zero external momenta. This standard computation (see Ref. [170] for an illuminating derivation) results in the master sum-integral This can be divided into a zero-temperature Coleman-Weinberg piece and a thermal piece as in Eq. (2.19), and can be evaluated by using with D = 4 − 2ǫ. The divergent 1/ǫ terms are cancelled by the counterterm part The thermal functions from Eq. (2.19) are expanded at high-T in d = 3 − 2ǫ with finite parts where z = m/T , a b = (4π) 2 exp( 3 2 − 2γ E ) and a f = 16a b . Here z-independent terms do not contribute to the dynamics of the phase transition and can be dropped.
With these tools -together with the daisy prescription for resummation as explained in Sec. 2.1 -one can write down the familiar result as a sum of a temperature dependent and a temperature independent, Coleman-Weinberg piece where both terms in the one-loop correction rely on the field dependent masses. Explicitly, this is where scalars includes the Goldstone modes, gauge bosons, Higgs degrees of freedom, and The thermal part is accompanied with daisy terms for bosonic fields where m 2 res. denotes the mass eigenvalue with resummed parameters. Note that the daisy contribution of the photon is non-zero due to thermal screening, but this does not give any field-dependent contribution. The gauge dependence manifests itself through the field dependent mass eigenvalues m. These we give in Appendix A.3 together with the resummed mass eigenvalues.
Finally, an alternative form for the one-loop part of the effective potential is provided by separation of the master sum-integral into the soft and hard parts following Eq. (2.15), and upgrading the mass eigenvalues in the soft part to their resummed versions. The soft parts can be evaluated as T J soft (m res. ) (and noting that in the reduced dimension the overall factor for the gauge parts is d − 1 and not D − 1) and the hard parts in the high-T expansion are where we have dropped a mass-independent piece. The one-loop master sum-integrals are expressed for example in Refs. [120,171].

Correlators from the effective potential
We have carried out a crosscheck of our results, utilising the fact that the effective potential is the generator of renormalised 1PI correlators of the Higgs field. We utilised two different methods to computeΓ φ † φ ,Γ (φ † φ) 2 , andΓ (φ † φ) 3 . On the one hand, we computed them by taking derivatives of the effective potential with respect to the external scalar field. On the other hand, we computed the same correlators directly in terms of Feynman diagrams. For this we utilised computer algebra tools to automate the diagrammatic calculation, as discussed in Sec. B.2. We emphasise that this is a strong crosscheck, since for example the six-point correlator in the diagrammatic calculation is a sum of O(10 2 ) different diagrams, each with permutated external legs. 29 When the scalar field is shifted by the background field φ/ √ 2, one has at one-loop order where the one-loop piece expands in powers of the background field, also utilising the high-T expansion, as in renormalised perturbation theory, where the hat denotes the renormalised correlation function. The diagrammatic one-loop piece V 1 can be computed as (note that for simplicity we include only hard contributions -and drop O(v 0 ) terms) with ζ(n) the Riemann zeta function and employing the shorthand notation We have included terms with ζ(5) and ζ(7) that are used in computing some leading (1-loop) corrections to scalar 8-and 10-point correlators. The correlators,Γ, can then be extracted from the coefficients of the expansion in the background field. One obtains the same pure scalar correlators as in Appendix. B.2.1. As emphasised above, this is a strong crosscheck of the automated diagrammatic calculation outlined in the following section.

B Dimensionally reduced SMEFT in three dimensions
The construction of a dimensionally-reduced effective theory requires the computation of various correlation functions in the high-temperature, unbroken phase. For the Standard Model and its simplest extensions (e.g. 2HDM [113,116,120] or Z 2 -symmetric real-triplet extension [117]), it is possible to perform these computations mostly by hand with little computer assistance. However, for more complicated BSM theories -or in order to reach higher orders of perturbation theory in simpler models -automation is inevitable. In the past, dimensional reduction in electroweak theories has been performed in Landau gauge since the computational effort related with a general covariant gauge is usually immense. Landau gauge greatly reduces the effort, as propagators (and corresponding master sum-integrals) simplify and the number of contributing diagrams is reduced immensely. However, ideally we would carry out all computations in a general gauge, to explicitly test gauge invariance. In this work, using tools presented in Ref. [171], we have fully automated the dimensional reduction of the SM in a general covariant gauge, and furthermore we have included effects of the extra SMEFT dimension-six operator (1.1). Using the general covariant gauge throughout the computation, shows explicitly that the matching relations -for the super-renormalisable part of the EFT (c.f. discussion in Appendix. B.2) -are gauge invariant. This justifies previous computations conducted in Landau gauge [63]. Before diving into the details of our computation -a combination of manual work and automation -we give a pedagogic introduction to the different steps that call for automation in a dimensional reduction computation. For this, we use the example of the 2-point correlator or self-energy of the SU(2) gauge boson, with an explicit diagram-by-diagram calculation.

B.1 Dimensional reduction for beginners: electroweak Debye mass
This section gives an explicit tutorial for a typical dimensional reduction calculation, by computing the SU(2) gauge boson self-energy at one-loop order in the unbroken phase, or the A a µ A b ν -correlator. The results for this correlator can then be used to obtain the Debye mass of the temporal component, as well as the field normalisations for both temporal and spatial fields. Note that a classic reference for these computations is Ref. [63]. Here, we follow Ref. [96] in which Appendix C.1 gives a diagram-by-diagram result for this correlator in Landau gauge, and all required master sum-integrals are listed in Appendix B and Feynman rules in Appendix A. Note that the computation of the A a µ A b ν -correlator resembles the calculation the of corresponding correlator in high-T QCD which yields the thermal gluon mass. This computation is found in Sec. 5.4 of the textbook Ref. [108] using Feynman gauge.
The outline of a typical dimensional reduction computation is the following: (see Sec. 3 of Ref. [171]) Step 1: Choose a model Lagrangian.
Step 2: Derive corresponding Feynman rules in the unbroken phase. Note that the computation of correlators is done in the symmetric -unbroken phase of the theory -in the gauge eigenstate basis where all gauge bosons and fermions are massless.
Step 3: Generate all diagrams for each required correlator, and use Feynman rules to compose expressions for each individual diagram and compute related symmetry factors.
Step 4: Perform all algebra contractions of Lorentz, isospin, Dirac etc. indices and manipulate sum-integrals to express all diagrams in terms of a basis of master sum-integrals.
Step 5: Evaluate the basis of required master sum-integrals.
Step 6: Match finite parts of the correlators to solve for the desired quantities of the 3d EFT.
For the SU(2) gauge boson self-energy, there is no new contribution from the dimensionsix coupling of the SMEFT at one-loop order. Therefore the computation is the same as in the SM.
Step 1 The relevant Lagrangian of step 1 is of the form where definitions and the exact form of the gauge field strength tensor, fermion structures, covariant derivatives, definitions of group indices etc. are found in Sec. 2.1 of Ref. [96]. Here, we upgrade the computation to a general covariant gauge instead of using Landau gauge. This Lagrangian suffices for the computation of the SU(2) gauge boson self-energy outlined here, but of course for the full dimensional reduction, the photon field, corresponding ghost field and Higgs potential should be added.
Step 2 For step 2, the relevant Feynman rules are listed in Appendix A of Ref. [96], with the replacement of the SU(2) gauge field propagator in a general covariant gauge, In addition, we use the following short-hand notations: for gauge field cubic and quartic self-interaction vertices. The Lorentz indices and adjoint isospin structures are those of Eq.(A.3) and Eq.(A.5) of Ref. [96], respectively. For the generation of these rules one can resort to any preferred Feynman rule generating software (e.g. FeynRules [172]). However, since the computation is conducted in the unbroken phase, the isospin structure of the fields is comparatively simple, allowing the economical possibility of formulating Feynman rules for the full fields, rather than for their individual components. This significantly reduces the number of different diagrams, but results in a non-trivial structure of the isopsin indices for some vertices (cf. (B.3)). Our in-house software generates the Feynman rules starting from the model Lagrangian. After going over to momentum space and symmetrising over fields, the final rules are fully symmetric in group indices and momenta.
Step 3 Moving forward to step 3, the SU(2) gauge field 2-point correlator, A a µ (K)A b ν (−K), is evaluated at one-loop level with corresponding diagrams: For illustration, we evaluate the pure gauge bubble and ring diagrams that explicitly depend on the gauge parameter ξ 2 , as well as the fermionic diagrams. The ghost and scalar diagrams do not depend on the gauge parameters ξ i so results for them can be read from the Landau gauge calculation in Appendix C.1 of Ref. [96]. In renormalised perturbation theory, the UV-divergence is cancelled by a tree-level counterterm interaction diagram. The external momentum K can be set to be purely spatial, i.e. K = (0, k), and soft K ∼ gT . The softness of external momenta allows a series expansion in K, and its K 0piece gives a contribution to the thermal mass, while its quadratic piece contributes to field normalisations between 4d and 3d fields of EFT; see Ref. [63,173]. Note that since the heat bath breaks Lorentz invariance only in the temporal direction -leaving spatial Lorentz (or rotational) symmetry intact -only the temporal part of the fields can obtain a thermal mass. However, thermal screening does effect the spatial fields in ways that do not break rotational symmetry, affecting their couplings and field normalisations when the hard scale is integrated out. The following only includes contributions from the hard scale, as at one-loop level zero-mode contributions trivially cancel against 3d contributions when correlators are matched.
The pure gauge bubble diagram reads (with 4-momentum P in the loop) where s is the symmetry factor of the diagram. From a computational point of view, the required diagrams of each correlator are generated using e.g. qgraf [174]. All momenta are shifted onto a canonical momentum basis already at the diagram level. Only then Feynman rules are inserted to compose expressions for each individual diagram.
Step 4 For step 4, while contraction over adjoint isospin indices is trivial, contraction over Lorentz indices is laborious already at one-loop level. For adjoint isospin δ aa = N 2 2 − 1 = 3, where N 2 = 2 is the fundamental isospin dimension. For Lorentz indices δ µµ = D = d + 1 in D spacetime dimensions. The pure gauge bubble reads Next, we separate the result into temporal µ = ν = 0 and spatial µ = r, ν = s (r, s = 1, . . . , 3) parts -note that cross-terms vanish due to odd sum-integrations -and manipulate sum-integrals in order to express everything in terms of one-loop master sum-integrals In the case at hand, we need the relation I 4b α+1,β+2 = 1 − d 2α I 4b α,β . For the calculation at hand, we find all relevant reduction relations in Appendix B of Ref. [96]. In order to utilise these, we additionally scalarise integrals by e.g. Σ ′ P prps P 4 = δrs d Σ ′ P p 2 P 4 which follows from 3d Lorentz invariance, and use trivial manipulations such as p 2 = P 2 − P 2 0 . Note that P 0 is always non-zero, as we integrate over non-zero modes only. Eventually, we arrive at an intermediate result for the pure gauge bubble Then, turning to the fermionic diagram, we find evaluate Dirac trace and contract Lorentz indices , (B.9) which becomes (B.10) Due to the momentum being soft, the propagator can be K-expanded up to quadratic order as Therefore, the fermion ring diagram reads, up to quadratic order in soft momentum K Again, all these sum-integrals can be expressed in terms of I 4b 1 and I 4b 2 after scalarisation of integrals such as Σ {P } PrPs(P ·K) 2 P 8 P 8 and by using the recursion relations from Appendix B of Ref. [96]. Eventually, this simplifies to Here, the spatial part is individually transverse and without a K 0 contribution, as spatial fields do not generate a thermal mass. Next, the pure gauge ring diagram reads (we denote here Q = P + K) (B.14) Performing contractions and expanding in soft momentum yields multiple terms, but nothing qualitatively new compared to previous diagrams that we have already mastered. An intermediate result in terms of familiar master integrals reads In its spatial part, the ξ 2 -independent piece is not individually transverse, but together with the corresponding ghost diagram, the sum of diagrams is rendered transverse. Also the momentum independent part of the full spatial correlator vanishes. On the other hand, the ξ 2 -dependent part is transverse, because ghosts do not contribute to that. Both ghost and scalar diagrams can be evaluated in a manner similar to previous diagrams, and since they do not posses a ξ 2 -dependence, we can read their Landau gauge results from Appendix C.1 of Ref. [96]. Finally, the counterterm diagram contributes by so that our correlators become UV-finite in dimensional regularisation. We emphasise, that counterterms within a finite-T computation, are recycled from zero temperature, since the UV structure of the hard mode contributions remains unaltered by high temperature. This means that the dimensional reduction computation can obtain all required counterterms and eventually β-functions along with the construction of the high-T EFT.
All algebraic manipulations encompassed in this step are handled using our favorite kitchen knife FORM [175]. This includes contractions of Lorentz, group (colour, isospin), and Dirac indices. Additionally we manipulate sum-integrals to express all diagrams in terms of a basis of master sum-integrals, shift onto different integral sectors, and employ integration-byparts methods [176] using a standard Laporta reduction [177]. Their algorithmic implementation is documented in Sec. 3.4 of [171].
Step 5 Proceeding, step 5 evaluates the basis of master sum-integrals in dimensional regularisation. Due to the hierarchy of scales, these sum-integrals can be expanded in powers of m/T , and hence evaluated as massless sum-integrals, leading to significant simplifications.
At the order that we work, NLO in powers of g 2 , only one-loop sum-integrals are needed to match the fields and couplings. Therein, computations are straightforward and give results in terms of Gamma-and Zeta-functions, see for example Appendix B of Ref. [96].
Two-loop sum-integrals are required for NLO matching of the masses of scalars and temporal gauge fields. At two-loop order there is only one master topology, the sunset diagram. A direct evaluations of the sunset topology sum-integrals can be found in Refs. [51,54,178,179]. We also recommend an illuminating Ref. [180]. However, our streamlined use of IBP relations [171,176,181], reduces all massless two-loop sum-integrals to products of one-loop sum-integrals. No two-loop masters are needed.
Going beyond NLO dimensional reduction requires higher loop sum-integrals. In this case the evaluation of master integrals becomes the most non-trivial part in the dimensional reduction pipeline. Automation of such computations is still in its early stages. For relevant computations of master sum-integrals at three-loop and higher loop orders, see Refs. [62,99,178,[181][182][183][184][185][186][187][188][189]. Further, see Ref. [110], for a three-loop computation of the pressure in the electroweak theory utilising the master integrals of the above references.
Step 6 At last, step 6 retrieves the final form of the correlators. After the ǫ-expansion, 1/ǫ-poles have been cancelled by the counterterm contribution and Here it is crucial that Π ′ -parts encode both the explicit ξ 2 -dependence and RG-scale dependence. It is indeed these dependencies that cancel contributions from other correlators related to the matching of corresponding EFT parameters, namely the 3d gauge coupling g 2 3 and portal coupling h 1 between A a 0 and 3d Higgs. We will show this cancellation explicitly in Eqs. (B.60)-(B.75). Finally, let us inspect the actual matching formula for the electroweak Debye mass. The matching of correlators can be achieved by equating effective Lagrangians in both 4d and 3d theories In this schematic matching example, we have in fact included both correlators at 2-loop level. At this level, one has to resum the zero-mode of the temporal gauge field by its one-loop thermal mass, and introduce a corresponding resummation counterterm interaction to the Lagrangian. As a result, the soft/hard mixing contribution of the correlator in 4d vanishes identically, and the 2-loop soft contribution of zero-modes matches exactly the loop corrections to the 3d correlator, i.e. cancelling soft and 3d parts in the above equation. By accounting for the relation between 3d and 4d fields one can finally solve Here we have highlighted how gauge dependence drops out between 2-loop hard contributions and the 1-loop field normalisation contribution, and how the final 3d parameter is (in a power counting sense) independent of the RG-scaleμ, as the running of the LO piece is cancelled by logarithms of the NLO piece. However, an important comment is necessary here: for the electroweak phase transition the NLO piece of m 2 D is not needed, as it contributes to the Higgs effective potential (or free energy) at O(g 5 ) and hence is of higher-order than we work. 30 The transition is driven by the Higgs field at the ultrasoft scale, and in the fact temporal A a 0 field can be integrated out in the second step of dimensional reduction going from the soft to the ultrasoft scale. Therefore, only the LO piece of m 2 D is needed and reads The explicit parameters of fundamental isospin dimension N 2 = 2 and fundamental colour dimension N c = 3 are set to their integer values henceforth. Note that the result (B.23) is not RG-invariant at O(g 4 T 2 ) due to the absence of the NLO contribution, but again this only contributes to the Higgs effective potential at order O(g 5 ). As a final remark, we note that for studies of the EWPT although 2-loop level matching for m 2 D is unnecessary, 2-loop level matching is important for the Higgs mass parameter, as it contributes to the Higgs effective potential at O(g 4 ). Indeed later this appendix encounters a similar matching relation for the Higgs thermal mass, with similar qualitative features of the cancellation of the gauge parameter and RG-scale.
As a final illustration, we show all the diagrams contributing to the Higgs thermal mass at one-loop level originating from the (φ † φ)-self-energy We highly recommend interested readers to embark on this computation themselves starting from the above diagrams. At leading order in the high-T expansion, this correlator readŝ and will result in a familiar one-loop thermal correction to the Higgs mass parameter, c.f. Eq. (A.57). The two-loop result for this correlator -with NLO mass correction -will be presented in Eq. (B.55). Next, we proceed from this pedagogic computation of thermal masses to the full dimensional reduction and construction of the 3d EFT that includes thermal screening effects to all EFT parameters.

B.2 Results for dimensional reduction of SMEFT
Dimensional reduction is performed by matching the infrared parts of correlators of the full 4d theory with those of the effective 3d theory. For leading order dimensional reduction, the 3d mass is computed at 1-loop and couplings at tree-level, i.e. couplings are only scaled by temperature. In terms of our power counting (g ′ ∼ g, g Y ∼ g, λ ∼ g 2 and c 6 ∼ g 4 ) viz. Eq. (2.13), this corresponds to O(g 2 ) accuracy for parameters at the soft scale after the first step of dimensional reduction, and O(g 3 ) at the ultrasoft scale after the second step of dimensional reduction. At next-to-leading order, i.e. O(g 4 ) accuracy at the soft scale, the 3d mass is computed to 2-loop accuracy and couplings to 1-loop.
Effective theory parameters are expected to be independent of the choice of gauge fixing parameters [60]. This is true at least up to O(g 4 ). In order to prove this, we perform dimensional reduction in a general covariant gauge and show below explicitly how the gauge parameters cancel. Additionally, 3d parameters are independent of the 4d RG scale in terms of our power counting. Any leftover scale dependence is of higher order than O(g 4 ), and if numerically important, signals a bad convergence of perturbation theory. The cancellation of gauge parameters and RG-scale provide a very powerful cross-check for the validity of the calculation.
The 3d effective theory which we match to is (note the Euclidean metric as this is a purely spatial theory), The τ a are the Pauli matrices. For simplicity we do not use a different notation for 3d and 4d fields, following the convention of Ref. [63]. The scalar potential in the soft scale 3d theory reads 31 and together with 3d gauge couplings g 3 and g ′ 3 , our task in dimensional reduction is to find these parameters in terms of 4d parameters and the temperature. Here g 2 3 , g ′  and, in Secs. B.4 and B.5 we adopt the Landau gauge, ξ 3,1 , ξ 3,2 → 0 + (taking this limit explicitly avoids certain IR divergences [127,135]). At the ultrasoft scale the temporal scalars A 0 , B 0 and C 0 are heavy and are integrated out. The remaining ultrasoft EFT parameters are differentiated from the soft EFT parameters with a bar, i.e.ḡ 3 ,ḡ ′ 3 ,μ 2 h,3 ,λ 3 ,c 6,3 ,c 8,3 , andc 10,3 . Although the ultrasoft theory differs only by these aforementioned changes, for clarity we show the ultrasoft action in full: The implicit gauge couplings areḡ 3 for the SU(2) andḡ ′ 3 for the U(1) sectors. The ultrasoft scalar potential reads and the gauge fixing corresponds to the soft scale in Eq. (B.28).
The dimensionally-reduced operator basis which we use is complete at dimension [GeV 2 ], but is incomplete at dimension [GeV 3 ], at which it contains only a single interaction operator, 3 . This is because this operator appears also in the 4d SMEFT, and hence c 6,3 receives a tree-level contribution of order O(g 4 ), whereas all other coefficients of dimension [GeV 3 ] operators start at one-loop level, at order O(g 6 ). Some examples of missing operators at dimension [GeV 3 ] in the dimensionally-reduced SMEFT are the following gauge invariant Higgs kinetic terms (we use a Warsaw-like basis, cf. Ref. [65,66,72]) There are also multiple dimension [GeV 3 ] operators of the Higgs coupling to temporal scalars at the soft scale, but these are also of order O(g 6 ).

B.2.1 Results for correlators
At NLO in dimensional reduction, higher dimensional operators with c 6 , c 8 and c 10 coefficients only affect pure Higgs correlators, while all other correlators are unaffected and therefore are already found in the literature. However, in order to demonstrate the cancellation of gauge parameters in parameter matching, we list all required correlators below. As explained in the previous section, we use an automated, computer algorithm computation of all correlators in the unbroken phase: diagram generation and computer algebra can be used to automate Lorentz contractions, Dirac traces and IBP-methods to reduce sum-integral structures to a basis or set of master integrals. This calculation is implemented via qgraf for diagram generation and FORM for algebraic manipulations and IBP reduction. Note that at NLO in dimensional reduction the scalar propagator can be treated as massless within twoloop order terms. This demonstrates the power of IBP reduction, since all two-loop integrals reduce to products of one-loop master sum-integrals.
We denote four-point correlators at zero external momentum by Γ, and two-point correlators by Π for the part evaluated at zero external momentum and Π ′ for the coefficient of the term which is quadratic in external momenta. 32 For the spatial gauge fields this is the transverse part of the correlator. For all correlators we only list contributions from the hard modes, as the soft contribution cancels with corresponding 3d contribution in the matching, as Sec. B.2.2 demonstrates. Note that all soft contributions at one-loop are finite, and do not require renormalisation. At two-loop, the scalar correlator requires resummation and soft contributions require renormalisation.
All correlators are computed in renormalised perturbation theory, and are therefore finite, which we denote by hats inΠ andΓ (note that we omit writing corresponding isospin, colour and Lorentz index structures explicitly). An illustrative example of renormalisation is the mixed scalar-gauge (φ † φAA)-correlator where the unhatted Γ sums all contributing Feynman diagrams, excluding counterterms. 33 Similar relations hold for other correlators, and emerge naturally in renormalised perturbation theory where counterterms are treated as interactions.
We start from one-loop contributions quadratic in the soft external momentum, and these read (L b/f are defined in Eqs. (A.77)-(A.78) and N f = 3 is the number of fermion families) (B.37) Again, due to the heat bath breaking the 4d Lorentz invariance, contributions to spatial and temporal gauge fields differ. Crucially, these quadratic momentum contributions to the SU(2) gauge field and to the Higgs explicitly depend on the gauge, and are required to render the final EFT parameters gauge independent in Sec. B.2.2. These contributions are evaluated at one-loop at O(g 2 ) order -which reaches the desired O(g 4 ) accuracy for the EFT parameters -and at this order there is no c 6,8,10 -dependence from higher dimensional operators.
Turning to contributions at zero external momentum, the pure Higgs correlators read Other correlators are independent of c 6,8,10 -couplings, and read The computation of the Higgs two-point correlator at two-loop level is the most complicated part of dimensional reduction. To illustrate this -and the power of our automated computation -we depict all the occurring diagrams in Fig. 12. Note that as this correlator is computed at zero external momentum, in Landau gauge many of these diagrams vanish trivially, significantly reducing the computational effort -and allowing even a manual, non-automated computation such as the one performed in Ref. [120] in case of the Two-Higgs-Doublet Model. By denoting where γ E is the Euler-Mascheroni constant, we obtain the result Again, the leading c 8,10 pieces are denoted with flowers and they appear at three-and four-loop orders. In order to obtain this expression, resummation is required to cancel the IR sensitive mixed hard/soft modes, that are non-analytic in the mass parameter µ 2 h . This resummation can be performed by adding and subtracting one-loop corrections to the mass parameter, in which case terms with plus signs contribute to the mass in the scalar propagator, while terms with minus signs are treated as (counterterm-like) interactions. For details of this procedure, see Refs. [63,120]. However, a shortcut to circumvent this procedure is provided by using IBP reduction to evaluate sum-integrals: since in dimensional reduction at NLO two-loop mass effects are of higher order, all two-loop sum-integrals can be treated as massless. Therefore, all mixed soft/hard contributions vanish trivially in dimensional regularisation, and non-analytic IR structures in need of resummation never appear.
Additionally, renormalisation of Π φ † φ at two-loop level is more complicated than for other correlators at one-loop. In fact, if one includes only hard contributions to this correlator, a divergence proportional to T 2 remains. This leftover divergence is cancelled by the divergence in the soft part of this correlator, and on the 3d theory side it corresponds exactly to the 3d mass counterterm, c.f. Eq. (B.123). What remains after cancellation of this divergence is the logarithm ofμ 3 visible in Eq. (B.55).

B.2.2 Parameter matching and gauge invariance
The matching of the parameters can be performed by equating effective vertices in both 4d and 3d theories. For this, the general relation between generic 4d and 3d fields ψ reads This relation can be derived from the condition that the physical, screened masses of the fields agree between the 4d and 3d theories [63]. For an illustrative example of matching, let us consider the scalar quartic coupling. The renormalised (φ † i φ j φ † k φ l )-correlators read (here we write isospin structure ∆ ijkl ≡ δ ij δ kl + δ il δ jk explicitly) The contribution from soft modes equals the 1-loop 3d contribution, i.e. Γ soft (φ † φ) 2 = Γ 3d (φ † φ) 2 so these contributions cancel each other. Effective vertices are then which leads to the following NLO matching relation once the relation between the 4d and 3d fields is taken into account. In a similar manner, we obtain a complete list of the other matching relations in terms of other n-point correlators: The key feature of these formulas is that all 3d parameters are gauge independent up to O(g 4 ): at one-loop level Debye masses and quartic temporal scalar self-interactions are immediately gauge independent since the correponding correlators are (field normalisation contributes to these parameters only at NNLO). In the other parameters we observe an explicit cancellation of gauge parameters between correlatorsΠ,Γ and field normalisationsΠ, when we insert the corresponding expressions. However, there is a notable exception, the higher order corrections we include for c 6,3 , c 8,3 and c 10,3 are gauge dependent, at orders O(g 6 ), O(g 8 ) and O(g 10 ) respectively. For e.g. for c 6,3 the corresponding correlatorΓ (φ † φ) 3 in Eq. (B.39) is gauge-dependent at subleading order, as the gauge-dependent parts proportional to ζ(3) are not cancelled by the field normalisation piece. A possible explanation of this leftover gauge-dependence is that the operator basis is incomplete. We have not included higher dimensional kinetic scalar operators of the schematic form (φ † Dφ) 2 where D represents covariant derivative. Our guess is, that when a complete operator basis is used, one can perform a field redefinition into a basis that is manifestly gauge invariant. However, we do not tackle this problem further here, but leave this topic for future research.
The final results for the 3d parameters at the soft scale read These formulas emphasise that the LO pieces run in terms of 4d RG-scaleμ. By applying corresponding β-functions, one can observe that this running cancels the explicit logarithmic scale dependence of the L b/f -terms. However, there remains a scale dependence which is formally of O(g 6 ) for the parameters of the Higgs and spatial gauge bosons, as discussed in Sec. 3.1. There is also a scale dependence at O(g 4 T 2 ) for the Debye masses m 2 D , m ′2 D and m ′′2 D . To cancel this scale dependence, the Debye masses should be evaluated at two-loop order. However, this scale dependence only contributes to the Higgs effective potential at O(g 5 ), and a full O(g 5 ) calculation goes beyond the scope of this paper; see for example Refs. [109,125]. In practice this leftover scale dependence is numerically insignificant since the running of the gauge couplings is small.
The above Eqs. (B.78) and (B.79), employ the following shorthand notation for terms proportional to ζ-functions (B.94) These contributions -with leftover gauge dependence -are formally O(g 8 ) and O(g 10 ) respectively, and we include them as they contribute at leading (one-loop) order to the respective 3d parameters. Finally, when the soft temporal scalar scalars are integrated out in the second step of dimensional reduction (c.f. Ref. [63]), the action of the final ultrasoft scale EFT is given in Eq. (B.29). The parameters of this ultrasoft EFT read (B.97) In this step of dimensional reduction, the two-loop matching of mass parameters requires resummation (by adding and subtracting the one-loop contribution to ultrasoft mass). Consequently the soft Higgs mass parameter does not appear inside the two-loop piece -despite mixed diagrams involving both Higgs and soft temporal scalars -since the effect of the Higgs mass is formally of higher order. We point out that this technical detail was overlooked in Ref. [120] in Eq. (3.45). For higher dimensional operator couplings, for simplicity we neglect contributions from B 0 and C 0 fields, as they are numerically insignificant. In fact, already corrections from A a 0 are heavily suppressed. At the ultrasoft scale we havē (B.99) (B.101) These formulas complete our construction of the 3d EFTs that we use to incorporate higher order resummations in our perturbative analysis of the phase transition in the SMEFT.

B.3 The 3d perturbative expansion parameter
For the case of the 3d EFT considered in this paper (cf. Eq. (B.29)), we show that theexpansion for the phase transition equals an expansion in powers of √c 6,3 . Near the critical temperature two minima are separated by a barrier, and are of similar heights. For the treelevel Higgs potential to show this structure, all three terms must be of the same order. From this one can derive that the broken minimum scales as φ ∝ −λ 3 /c 6,3 . To expand around the broken phase or around the critical bubble, we shift the Higgs by a background field, φ → 0, Φ/ √ 2 + φ, a saddle point of the tree-level scalar action, where {ϕ α } = {φ, A r , B r , η,η} runs over all of the fields, and all factors ofc 6,3 have been made explicit. This shows that the effective loop expansion parameter is proportional to c 1/2 6,3 . This perhaps surprising observation, thatc 6,3 controls the loop expansion in the coupled gauge-Higgs theory (even e.g. the (A r ) 4 interaction!), is a special property of the truncated 3d theory we are considering, and would not hold in a more general truncation including, for example, higher order derivative terms. As regards theḡ 2 3 ,ḡ ′ 2 3 andλ 3 dependence of physical quantities computed in the 3d loop expansion, by scaling out overall powers ofḡ 2 3 to fix dimensions, one can see that only the ratiosλ 3 /ḡ 2 3 andḡ ′ 2 3 /ḡ 2 3 may arise to modify the expansion parameter. In our analysis of the SMEFT,c 6,3 is naturally the smallest parameter, and hence the 3d loop expansion will generally converge well, at least around the broken phase.
The structure of Eq. (B.105) is not modified in the symmetric phase, meaning thatc 1/2 6,3 acts as the loop counting parameter there too. However, in the symmetric phase the treelevel mass of the 3d gauge bosons is zero, which leads to IR divergences, and consequently nonperturbativity. A fuller discussion of this is given in Sec. 3.6.

B.4 Two-loop thermal effective potential
For the SU(2) gauge theory, with a Higgs field in the fundamental representation, the effective potential has been calculated to two-loop order in Refs. [60,127,128]. Ref. [60] also generalises this to include the U(1)-hypercharge. For the SMEFT, the (φ † φ) 3 term introduces corrections to the Feynman rules. However, there are no new connected vacuum diagrams to be added, as we now show. Using the usual topological identities for connected graphs (see for example the chapter on divergences and regularisation in Ref. [190]), one can derive the following equation where δV is the vacuum counterterm, and all counterterms arise at two-loop order. Due to the presence of the dimension-6 operator, the 3d EFT is not super-renormalisable, and the 2-loop counterterms are not exact, as in the pure SM 3d EFT. Note that as we utilise the gauge-invariant -expansion for our 3d computations, we are free to fix a gauge in calculating the 3d effective potential. In Landau gauge, the one-loop The two-loop contribution to the effective potential in the 3d EFT is straightforward to include, by following Refs. [60,118,127]. Since new contributions fromc 6,3 appear only through the modified mass eigenvalues m χ,3 and m φ,3 , and through pure scalar vertices, where the latter affects only the (SS) and (SSS) topology classes below. where different topology classes are illustrated in Fig. 13 and their results read (again, in Landau gauge)

C Estimates for the nucleation prefactor
In a semiclassical evaluation, the bubble nucleation rate takes the following form, where S c is dimensionless and A has mass dimension 4. The nucleation prefactor, A, is discussed in Sec. 3.5. We do not calculate it explicitly in this paper, but instead give some rough estimates for it. The prefactor naturally splits into the product of two parts, a dynamical part A dyn and a statistical part A stat 34 , as shown in Eq. (3.14). The statistical part of the prefactor has mass dimension 3. Although it is difficult to calculate, the definition of this part is agreed upon in the literature. It is given by a ratio of functional determinants of the second-order fluctuations around the critical bubble and the symmetric phase, where S is the action, α and β run over the fields which fluctuate about the critical bubble, and we are using DeWitt notation [191] for the functional derivatives. In the 3d theory α and β run over {φ, A a r , B r , η a ,η a , η,η}. At this stage the ratio of determinants is formal, as there are zero and negative eigenvalues which must be dealt with separately. The imaginary part arises due to the presence of a single negative eigenvalue, for which the corresponding integral must be carried out by analytic continuation [88].
The dynamical part of the prefactor arises for thermal, and not for vacuum, transitions and is essentially an inverse timescale for the critical bubble to grow, This should depend both on the exponential growth rate of undamped linear perturbations to the bubble radius, and on the damping due to the thermal bath. However, its precise definition is not widely agreed upon (see e.g. Refs. [86,121,133,164,192,193]).

C.1 Dynamical prefactor
We consider two different estimates of the dynamical prefactor. The first is essentially parametric, and follows from the dominance of infrared gauge bosons in the time evolution of the critical bubble. The second is a detailed formula which relies both upon on the thin wall approximation and a hydrodynamic approximation.

C.1.1 Infrared gauge boson dominance
Here we follow Refs. [133,192,194,195]. The exponential growth of the critical bubble is checked by the parametrically slower evolution of the infrared modes of the gauge bosons, as shown concretely in lattice simulations (see Fig. 11 in [133]). Thus, the dynamical prefactor should be of order g 4 T , the inverse timescale for the evolution of the infrared modes of the gauge bosons [194]. This is the same reason why the sphaleron rate in the symmetric phase is O(g 10 T 4 ) 35 . More precisely, up to an O(1) multiplicative factor, the result is where σ SU(2) ≈ 0.477 T is the "colour" conductivity of the weak sector. The addition of powers of 4π here follows Ref. [133]. We use this for our first estimate of the dynamical prefactor.

C.1.2 Hydrodynamic approximation
Refs. [104,105,164,196] adopted a coarse-grained, hydrodynamic approach to calculating the nucleation prefactor. For the purposes of describing bubble nucleation, it was assumed that the dynamical degrees of freedom could be modelled by the macroscopic energy density, rather than working directly with the quantum fields. Further the thin-wall approximation was adopted.
With the proceeding caveats and approximations, the expression for the dynamical prefactor obtained is [104,105], where σ is the surface tension, η S is the shear viscosity, r b is the bubble radius and ∆w is the change in the enthalpy, or latent heat, of the transition. Here the bubble radius can be consistently replaced with its thin-wall expression, r b = 2σ/∆V . We have used η ≈ 82.5T 3 following Ref. [105]. Note that in Refs. [104,105] the statistical prefactor was also estimated, however in this case the contribution of the IR geometric deformations of the bubble was wrongly dropped (see Sec. C.2). In the SMEFT in the 3d approach, the remaining pieces of the expression are, σ T =λ In the 4d approach, these expressions are replaced by expressions which are evaluated numerically. Again, one must ad hoc throw away the imaginary part of the effective potential.

C.2 Statistical prefactor
The statistical prefactor has been computed numerically for the thermal phase transition of a real scalar field by several groups [145,146,159,197,198], and various approximation schemes were proposed in Refs. [145,155,198]. For more complicated theories, such as the electroweak sector, related computations have been performed for the sphaleron rate [192,[199][200][201] and the Higgs vacuum decay rate [162,202]. A complete calculation of the statistical prefactor for the bubble nucleation rate in the electroweak sector, or related BSM extensions, is made difficult by the radiatively induced nature of the transition (see Sec. 3.5), and associated infrared divergences [203]. To our knowledge, there is no complete calculation in the literature of the statistical prefactor in the electroweak sector, or related BSM extensions. However approximations to it have been proposed in Refs. [104,105,203]. The statistical prefactor can be further broken up into the contributions from the zero modes, the negative mode, and the positive modes, stat . (C.8) The contributions from the zero modes, and the negative mode can be calculated using standard methods. It is A (+) stat which is most difficult to calculate. The factor of one over the volume of space is included in the zero mode contribution.

C.2.1 Zero modes
Some knowledge about the statistical prefactor can be gained simply by knowledge of the zero modes of the fluctuation spectrum about the critical bubble (see e.g. Ref. [129]). In particular, we will be able to determine the dependence of the prefactor onc 6,3 , which is, in all cases, the smallest parameter in the 3d effective theory.
Zero modes arise due to global symmetries of the theory which are broken by the bubble configuration. Integration over the zero modes can be carried out with the method of collective coordinates [204] (see also Refs. [192,199] for similar computations). There are 3 zero modes corresponding to the breaking of spatial translations. These are well known and integration over them results in the following volume and Jacobian factors, where on the left hand side we have shown the form of the measure on these modes before the collective coordinate transformation and V is the volume of space. There are also 3 zero modes due to the breaking of the global symmetry part of SU(2) × U(1) → U(1). Note that these global symmetries are not broken by our choice of general covariant (or Fermi) gauge, Eq. (B.28), unlike in the case of R ξ gauges. Integration over these zero modes has been carried out in Appendix C of Ref. [130], and results in, For bubbles which are not in the thin wall regime [139], the only length scale entering the bounce solution is 1/μ h,3 , which is in turn of order 1/μ h,3,c ∼ 2 √c 6,3 /λ 3 when the supercooling is not parametrically large. Noting this, and scaling the operators in the ratio of determinants byμ h,3 , we arrive at where the function κ is given by (C.14) Here the eigenvalues of the operators have all been scaled byμ 2 h,3 to be dimensionless and the dash on det ′ denotes that the six zero modes are removed. In the second line S ,W G and S ,ZG refer to the second derivative matrices in the subspaces spanned by the W bosons and their Goldstone modes and the Z boson and their Goldstone modes respectively. In going from the first to the second line we have used that, in the Landau gauge, neither the ghost propagators nor the photon propagator depend on the background Higgs field. Thus, only the physical Higgs particle, the W and the Z bosons and their respective Goldstone modes contribute to the statistical prefactor [162,205].
Our first approximation to the statistical prefactor is simply κ = 1, which we call our thick wall approximation. It is accurate up to the multiplicative function, κ, which should be of O(1) if the ratios of couplings are themselves of O(1). In reality these coupling ratios lie in the region ∼ 0.1 − 0.4, and hence one might expect corrections to the nucleation rate of a few orders of magnitude. Note that we have not included ac 6,3 dependence for κ, as thec 6,3 dependence of the prefactor is fixed by our assumption thatμ h,3 is parametrically of the same order asμ h,3,c , i.e. that the supercooling is not parametrically large. At larger supercooling thec 6,3 term becomes irrelevant to tunnelling, as the scalar field only tunnels to φ ∼μ h,3 / λ 3 , much short of the broken minimum, hence κ has noc 6,3 dependence in this case too.

C.2.3 Nonzero modes: thin wall approximation
For very small supercooling from the critical temperature, the bubble radius, r b , will be much larger than the thickness of the bubble wall, r w ∼ 1/μ h, 3 . In this case one can make a thin wall approximation. The large hierarchy of scales, r w /r b ≪ 1, leads to a multiplicative correction to the statistical prefactor unique to the thin wall limit. This is made up by lowenergy geometric deformations of the bubble, with eigenvalues which scale as ∝ 1/r 2 b [122,[206][207][208]. For thermal transitions in the 3d effective theory it is which when combined with the contribution of the zero modes in the thin wall approximation, where we dropped the O(1) constant, given the uncertainty in the contribution of the positive modes. Note that this result differs from Ref. [104] (and consequently Refs. [6,105]) for two reasons: first, we have included the additional zero modes from the breaking of custodial symmetry and second, we have included the parametric contribution of the positive low-lying modes, i.e. those considered in this Ref. [206]. These both affect the statistical prefactor at leading order. Our result agrees with Refs. [154,209].
Eq. (C.17) is applicable only after performing dimensional reduction. The reason is that only spatial deformations of the critical bubble are included in A (+) stat . A similar analysis is carried out included the deformations of the bubble in the compact thermal direction. This was carried out in Ref. [122] where it was shown that they contribute exponentially, This result suggests that, at least in the thin wall approximation, the nucleation prefactor gives a much larger correction to the nucleation rate if one does not perform dimensional reduction.

C.2.4 Derivative expansion
Modes with wavelengths much shorter than the critical bubble allow for a derivative expansion of the fluctuation determinant. In this case, the leading order approximation takes the background Higgs field to be locally constant. The wall of the critical bubble has a width of order the inverse Higgs mass. As such the prefactor of particles which are much heavier than the Higgs can be well approximated in the derivative expansion.
The magnetic gauge bosons can never be heavier than the Higgs in the symmetric phase, as they are massless, at least in perturbation theory. However, for |λ 3 | ≪ḡ 2 3 , the gauge bosons are much heavier than the Higgs in the broken phase. For sufficiently strong phase transitions, the broken phase contributions are expected to dominate over the symmetric phase ones as discussed in Sec. 3.6, and hence one might expect the derivative expansion to give a reasonable approximation for the contribution of the W and Z bosons to the nucleation prefactor, as in the following, (C.19) where the sum over α runs over the modes of the W , Z and corresponding Goldstone bosons. This approach has been taken for the sphaleron rate in Refs. [199][200][201]210], and was shown to give a good approximation, at least for the case |λ 3 | ḡ 2 3 (see in particular Fig. 1 in Ref. [201]).
On the other hand, the derivative expansion is never formally justified for the prefactor of fluctuations of the Higgs particle itself, or for lighter particles. If one were to apply the derivative expansion to the Higgs itself, the negative effective mass near the top of the potential barrier would result in a spurious nonzero imaginary part (not systematically related to the expected imaginary part from analytic continuation [88]). Nevertheless, ploughing on and applying the derivative expansion to the Higgs particle fluctuation determinant, setting imaginary parts to zero by hand, one can derive a rough order-of-magnitude estimate of the statistical prefactor. The leading order contribution is taken to be, where we have followed Ref. [145] in extracting, ad hoc, the zero-mode contribution, Eq. (C.11), and have extracted powers of the massμ h,3 to make up the dimensions. Note the necessity of the introduction of the real part, justified on practical grounds. An analogous expression was shown in Ref. [145] to result in a good approximation of the statistical prefactor for a real scalar field in the case of thin-walled bubbles.