Robust approach to thermal resummation: Standard Model meets a singlet

Perturbation theory alone fails to describe thermodynamics of the electroweak phase transition. We review a technique combining perturbative and non-perturbative methods to overcome this challenge. Accordingly, the principal theme is a tutorial of high-temperature dimensional reduction. We present an explicit derivation with a real singlet scalar and compute the thermal effective potential at two-loop order. In particular, we detail the dimensional reduction for a real-singlet extended Standard Model. The resulting effective theory will impact future non-perturbative studies based on lattice simulations as well as purely perturbative investigations.


Introduction
A strong first-order cosmic phase transition (SFOPT) is a violent process that can trigger the generation of a primordial gravitational wave (GW) background (cf. [1,2] and reviews [3][4][5]). Gravitational waves from astrophysical sources have been detected by Earth-based detectors LIGO and VIRGO for binary black hole [6][7][8] and neutron star mergers [9][10][11]. Their success and the mission to probe evidence of relic gravitational waves from the early Universe have sparked interest for space-based gravitational wave observatories such as LISA [12], BBO [13], Taiji [14], and DECIGO [15]. A detection of such a relic GW background could scope the underlying theories of particle physics complimentary to collider physics [16][17][18][19].
The electroweak phase transition (EWPT) is a smooth crossover in the minimal Higgs sector of the Standard Model (SM) [20][21][22][23][24]. Therein, the observed Higgs mass of 125 GeV [25,26] exceeds the requirements for a SFOPT which precludes both the production of a cosmic GW background and electroweak baryogenesis [27]. The latter is a mechanism to produce the baryon asymmetry during the electroweak phase transition [28,29]. New beyond the Standard Model (BSM) physics can alter the character of the electroweak symmetry breaking towards a SFOPT. To this end, new particles need to be sufficiently light in the vicinity of the electroweak (EW) scale and strongly enough coupled to the Higgs. This indicates that such BSM theories offer theoretical targets to guide future high-energy collider experiments [30].
Thermal field theories are plagued by the infrared problem [100]. Their perturbative description of long distance modes is invalidated at high temperatures due to high occupancies of bosonic modes. Nevertheless, perturbation theory is still widely used when reorganising the perturbative expansion by resummation, such as in hard thermal loop perturbation theory [101] and daisy resummation [102].
Dimensional reduction implements the required resummations automatically upon perturbatively constructing the 3d EFT. Nonetheless, it is customary to study the EWPT in terms of the thermal effective potential [159,160] computed directly with other resummation schemes [102,161]. While improved two-loop computations exist [102,[162][163][164] (cf. also refs. [165,166]), it is typical for recent EWPT literature to implement a daisy-resummed thermal effective potential only at one-loop level. However, the infrared problem persists and these fully perturbative studies of the EWPT are severely limited with their setbacks often underestimated. Even their qualitative description can -and often will -fail. In contrast to lattice studies, transitions are often realised as (weak) first order since a crossover character is incompatible with perturbation theory.
Describing the EWPT thermodynamics all the way by a non-perturbative simulation poses a formidable task. This roots in analytical challenges related to the construction of required 3d EFTs and foremost excessive computational cost of simulations. Still, many perturbative studies of the EWPT could be significantly improved by employing perturbation theory within the dimensionally reduced 3d EFT. While this approach is entirely perturbative and hence incapable of solving the IR problem, it allows for systematic resummation and straightforward computations at two-loop level [167,168]. Thereby it supersedes the one-loop daisy-resummed thermal effective potential; see refs. [152,156,157] for recent direct comparisons. Indeed, this approach to perturbation theory was advocated already in ref. [135].
The dimensional reduction can be largely automated and the careful matching to multiple individual BSM theories streamlined. The task has been tackled recently [157,169] and in this work at hand. As a consequence one can exploit the universality of the resulting dimensionally reduced EFTs to efficiently examine the parameter space of different BSM theories. This article combines these recent developments to extend previous work [149] for the xSM -a flagship model that is attractive for particle cosmology due to its minimal nature. Based on the construction of the 3d EFT of the xSM, its applications [170,171] chart a course of a state-of-the-art analysis of the EWPT thermodynamics. Thereby, perturbative scans that utilise a 3d EFT approach can guide non-perturbative simulations that finally solve the IR problem.
This article is organised as follows. Section 2 reviews the computation of thermodynamics in generic scalar driven phase transitions, in particular focusing on the use of dimensionally reduced effective theories. Section 3 is a pedagogic tutorial to the construction of such a 3d EFT and computes the thermal effective potential at two-loop order, in the simplest case of a real singlet scalar. In sec. 4, for the first time, we generalise the dimensional reduction to the xSM. Finally, sec. 5 discusses our results and outlook, while several technical details relevant to our computation are collected in the appendices.

Thermal phase transitions
The focus of this paper is to take steps towards a state-of-art determination of the cosmic phase transition thermodynamics for individual BSM theories with non-minimal Higgs sector. Direct ab-initio lattice calculations are, however, not feasible for a completely realistic 4d description of the thermodynamics of electroweak phase transitions [172] due to problems related to chiral fermions. 1 One alternative approach are non-perturbative simulations of dimensionally reduced effective theories. Following this idea, we survey the required technology on a generic level in the following section.

Down the pipeline
Several steps have to be considered to accurately predict gravitational waves from cosmological phase transitions. To this end, we illustrate a "pipeline" ranging from the collider phenomenology of BSM particle physics models to a primordial, stochastic gravitational wave background. Following a comprehensive ref. [174], we display different steps of this pipeline in fig. 1 (ibid. ref. [174]).
From a theoretical standpoint, it is natural to start by defining the Lagrangian of the corresponding BSM theory. In our case of interest, the BSM field content enters as a nonminimal Higgs sector which contains one or more scalar fields. In general, the scalar fields can occur in any representation of the SU(2) symmetry and possess other symmetries and couplings to new gauge field or fermion content in a dark sector. However, several alternatives are conceivable (see refs. [174,175]), such as models for holographic phase transition [176][177][178] or Composite Higgs scenarios [179][180][181]. The pipeline constitutes the following steps: Step (A): Relating collider signatures and BSM theory Lagrangian parameters. The Lagrangian (running) parameters are related to physical observables such as pole masses and mixing angles in zero-temperature perturbation theory. Then, actual collider signatures include production cross-sections and a relative shift in the Higgs couplings from their SM predicted values, 2 and can constrain the available parameter space for phase transitions. The relation of the EWPT and collider physics is further discussed in ref. [30].
Step (B): Equilibrium thermodynamic properties as a function of BSM theory parameters. The former include the character of transition (crossover, first-order etc.), the critical temperature (T c ) and latent heat (L/T 4 c ). They are encoded in the free energy of the system which is associated with the thermal effective potential in perturbation theory. Due to IR sensitivities at high-T , this step requires non-trivial resummations compared to perturbation theory at zero temperature, and eventually non-perturbative techniques.
Step (B) is the main focus of the remaining sections of this article.
Step (C): If the phase transition is of first order, it proceeds by nucleation and expansion of bubbles of the broken phase in the presence of a surrounding plasma [183][184][185][186]. The bubble nucleation rate can be computed in a semi-classical approximation from the effective action which includes quantum and thermal corrections. The relevant quantities [187][188][189][190][191] are the Hubble parameter (H * ) or temperature (T * ) when the phase transitions completes, its inverse duration (β), strength (α) at T * and the bubble wall velocity (v w ). The exact definitions and derivation of these quantities are detailed in e.g. refs. [157,174,192], and in particular [193][194][195][196][197] for the bubble equations of motion and v w . Also non-perturbative methods for nucleation have been developed [198,199] as an alternative to perturbation theory.
Step (D): Numerical, large scale lattice simulations of relativistic hydrodynamics; cf. refs. [200][201][202][203][204][205]. The parameters that describe the phase transition dynamics, (T * , α, β/H * , v w ), are input to simulations of colliding bubbles, cosmic fluid and sound waves after the phase transition completes. These determine the GW power spectrum. In practice, the approximate, analytical power spectrum has been solved from such simulations in terms of a generic set of input parameters. For an application of this, ref. [174] has devised the online tool PTPlot.
Step (E): A detectable GW background signature depends on the architecture of the detector in addition to the predicted stochastic GW power spectrum. The determination of the signal-to-noise ratio for a predicted signal at LISA is specified in ref. [174].
Step (F): A necessary condition for the EW baryogenesis [27] are first order phase transitions occurring via bubble nucleation. For reviews cf. [28,29,206]. The generation of a baryon asymmetry could be realised when new BSM sources of C and CP violation are invoked and baryon number violating sphaleron transitions are sufficiently suppressed in the broken phase. The latter can be associated with sufficiently strong transitions.
Next, we detail step (B) starting by a brief summary of the technique of high-temperature dimensional reduction.

Dimensional reduction for a high-temperature 3d effective theory
High-temperature dimensional reduction encodes the IR physics of the high-temperature plasma in an effective three-dimensional theory to describe long wavelength phenomena. In the context of electroweak theories, classic references are [109][110][111] but we also refer [207][208][209].
The equilibrium thermodynamics of a thermal field theory is described by an evolution in imaginary time (τ ). Therein, bosonic (fermionic) fields satisfy (anti-)periodic boundary conditions with period τ = 1/T and can be decomposed into bosonic and fermionic Matsubara [210] modes where k is a three-dimensional (3d) momentum. In other words the resulting theory is a 3d one with an infinite tower of modes each carrying a mass ω 2 n corresponding to the Matsubara frequency of mode n.
This system can be studied in an effective theory formulation. In that EFT the central degree of freedom is the static bosonic 3d zero mode (ω B n=0 ) of the original four-dimensional (4d) field. The remaining non-zero modes of scale ∼ πT can been integrated out. This is the dimensional reduction step which is based on the high-temperature scale hierarchy In the scale hierarchy we introduced a power counting parameter g defining the hard (πT ), soft (gT ), and ultrasoft (g 2 T /π) scales. While the scaling of the hard scale is a direct consequence of the Matsubara decomposition, the soft and ultrasoft scale are pertinent to collective plasma effects. Based on this hierarchy one can invoke a high-temperature expansion m ψ /T ≪ 1 for generic scalar fields ψ, whereas gauge bosons and fermions are massless in the unbroken phase.
In hindsight of the ensuing studies of a real scalar field, we establish the formal scaling λ ∼ g 2 for the scalar quartic coupling λ which is based on their appearance at one-loop. In gauge field theories this power counting parameter is often set to be the gauge coupling. For a scalar field we assume the original mass squared parameter (µ 2 ) to behave as µ 2 ∼ λT 2 ∼ (gT ) 2 which implies that the mass of the 3d soft mode (µ 2 3 ) is thermally corrected by µ 2 3 ≃ µ 2 + (gT ) 2 at leading order. Phase transition physics can often be studied at the ultrasoft scale by a simplified 3d EFT, where the soft scale has been integrated out. In fact, the transition point resides near a vanishing µ 2 3 where thermal loop corrections cancel the tree-level part. At this point the 3d mass scale is formally of the next natural order which is the ultrasoft one µ 2 3 ∼ (g 2 T ) 2 where soft modes are screened. The corresponding soft degrees of freedom are the temporal (adjoint) scalar fields which are remnants of the zero components of gauge fields and induced by the broken Lorentz symmetry from the heat bath. They remain soft in the vicinity of the transition point. For this second step of dimensional reduction, see ref. [109].
The EFT is constructed by determining the operator coefficients of the effective Lagrangian. In practice, these parameters follow from matching correlation functions of both the fundamental 4d theory and effective 3d theory. The generic rules of this procedure were established in refs. [109][110][111] and applied recently [157,158]. This construction of the 3d EFT by dimensional reduction is completely infrared-safe. In the matching of correlators, the IR and 3d contributions cancel each other and only the hard scale (non-zero Matsubara modes) contributes. The corresponding sum-integrals over non-zero modes are IR-regulated by non-vanishing Matsubara frequencies at high temperature. Hence, the dimensional reduction defers the IR problem of high-temperature bosonic perturbation theory to the 3d EFT.
At next-to-leading order (NLO) dimensional reduction, couplings are matched at oneloop and masses at two-loop level. This ensures a O(g 4 ) accuracy in the established power counting. To fully match this accuracy, the running parameters have to be related to physical observables at one-loop order in zero temperature perturbation theory. In 3d perturbation theory the effective potential is computed at two-loop order. Notably, the frequently used 4d daisy-resummed thermal effective potential at one-loop includes some -but crucially not all -O(g 4 ) contributions [102].
Instead of detailing the generic procedure in later sections, we choose an alternative approach. In an explicit hands-on demonstration, we dimensionally reduce a scalar field theory in sec. 3, and generalise it to the xSM in sec. 4 and appendix A.

Approaches to thermodynamics of thermal phase transitions
Let us assess the main approaches to access the thermodynamics of the thermal phase transitions in electroweak theories. See also similar summaries in sec. 2 of ref. [135], sec. 2.2 of ref. [140] and sec. 1 of ref. [158]. The following approaches are illustrated in fig. 2: Perturbative effective potential with daisy resummation.
Non-perturbative lattice simulation of 3d EFT. Robust approach combining perturbative dimensional reduction and non-perturbative (Monte Carlo) methods.
The individual steps encompass: (a) Relating physical parameters (such as pole masses) and Lagrangian (running) parameters at zero temperature. Often the "4d approach" uses only tree-level relations (e.g. refs. [65,211]), but in order to match the accuracy of dimensional reduction at NLO O(g 4 ), one-loop vacuum renormalisation is required [109,156,165].
(b) Perturbative computation of the thermal effective potential [159]. Frequently performed at one-loop, with leading order daisy resummation [102,161]. Two-loop computations Monte Carlo simulation Figure 2: Three different approaches towards the thermodynamics of the electroweak phase transition. We focus on the purely perturbative 3d approach with steps (d) and (e) for a real scalar theory and the Standard model supplemented by a real scalar singlet.
are discussed for e.g. in refs. [162,163,165,166]. This computation suffers from the IR problem the most, and additionally can contain a dramatic artificial RG scale dependence if two-loop thermal masses are unaccounted [171].
(c) Computation of thermodynamics. At the critical temperature the minima of the effective potential are degenerate and thermodynamic quantities are obtained by differentiation with respect to temperature, viz. latent heat. Model-independent tools to locate degenerate minima have been implemented numerically in software including CosmoTransitions [212], BSMPT [211], and PhaseTracer [213]. It is worth noting that location of minima of the effective potential are not gauge invariant. Thus, this computation frequently introduces unphysical estimates for thermodynamics as discussed in refs. [157,166,214] and also ref. [215] (in 3d EFT context).
(d) Dimensional reduction to a 3d EFT. See refs. [109][110][111] and also recent refs. [157,158]. It is perturbative and IR-safe, since only the hard scale is integrated out, and systematically implements all required resummations. Furthermore, dimensional reduction at NLO is analytically independent of the 4d renormalisation scale up to that order, see ref. [109,171]. This decreases the theoretical uncertainty in perturbation theory. A concrete computation is displayed in secs. 3, 4 and appendix A.
(e) Computation of 3d effective potential, see refs. [167,168,215]. The computation in the 3d EFT simplifies significantly compared to 4d because sum-integrals are replaced by vacuum integrals in d = 3 − 2ǫ spatial dimensions. Hence, it straightforwardly extends to two-loop order, cf. sec. 3.4. For recent applications, see refs. [152,157]. Even the three-loop effective potential has been computed for a pure scalar theory [216] and applied recently [158].
(f ) Computation of thermodynamics from 3d effective potential. Again a pathological gauge-dependent analysis can be based on degenerate minima at the transition point. However, also a gauge invariant treatment is possible, in terms of gauge invariant condensates [103,157] or the pressure in -expansion. However, IR divergences arise at two-loop order for a radiatively generated transition [104], compromising the analysis [215]. On the other hand, these IR singularities are avoided in presence of a barrier at tree-level, and a manifestly gauge invariant treatment for the thermodynamics can be obtained in perturbation theory; see ref. [157].
(g) Lattice-continuum relations; see refs. [103,104,158,217,218]. The Lagrangian parameters of the lattice discretisation need to be related to those of the continuum theory. Thus, the results of Monte Carlo simulations can be associated with the 3d continuum theory and via dimensional reduction to temperature and physical parameters. This can be done by computing and equating effective potentials in both discretisations, to two-loop order. In super-renormalisable theories without higher dimensional operators, all divergences arise at finite loop order and hence relations between continuum and lattice are exact. However, this aggravates in the presence of higher dimensional operators as the 3d theory retains renormalisability but loses super-renormalisability. It remains a future challenge to overcome this technical issue. Note that in lattice gauge theories there is no need to fix the gauge, and the treatment is automatically gauge invariant by construction [219].
(h) Monte Carlo lattice simulations of spatial 3d EFT on finite volume and lattice spacing. 3 Arbitrary field configurations are evolved -usually by a colourful cocktail of update algorithms -to form a Markov chain converging to a Boltzmann probability distribution. Thereof, physical quantities can be measured such as scalar and gauge condensates, and correlation lengths. Many autocorrelation times are measured to ensure that statistical errors remain small. At the transition point, the system is equally likely to occur in any of the phases. Multi-canonical methods in first order transitions ensure that the system can efficiently sample all phases while not getting stuck in one.
(i) Extrapolate simulations of finite volume and fixed lattice spacing to the continuum. This corresponds to infinite volume and vanishing lattice spacing, thermodynamic and continuum limits, respectively. In practice, several lattice spacings are needed, each with several different volumes. This rapidly becomes computationally expensive and even a single parameter space point requires a large number of individual simulations. Furthermore, oftentimes manual effort is required to fit a proper continuum extrapolation to the data instead of an elephant.
This article focuses specifically on steps (d) and (e) which are detailed for a real scalar theory in sec. 3 and a real singlet scalar coupled to the SM in sec. 4 and appendix A. The full nonperturbative path of the real scalar theory is presented in ref. [158], where the corresponding results are compared with three-loop 3d perturbation theory. Finally, let us summarise the different approaches and describe some of their merits.

4d approach
The 4d approach is the accustomed "bread and butter" approach with the advantage of its conceptual simplicity. At one-loop order, a closed form expression for the effective potential is available in terms of mass squared eigenvalues and one-loop thermal mass corrections, which straightforwardly automates to different models. In addition, numerical tools for thermodynamics (e.g. minimisation of potential) have been developed [211][212][213]220] and can scan large regions parameter space of BSM models. However, the 4d approach suffers from the IR problem of perturbation theory [100] and is often plagued with large inaccuracies and theoretical uncertainties [21,103,104,156,157]. In particular, weak transitions are poorly described by perturbation theory and are sometimes even qualitatively mistaken. Especially, crossover transitions are not predicted at all and it is not expected to determine the critical temperature accurately since it is highly IR-sensitive. Conversely, large couplings are often required for strong transition and can compromise the perturbative expansion, even at zero temperature [156]. Consistent ( -)expansions leading to gauge invariant results are oftentimes unavailable, since they require the knowledge of higher order contributions. Furthermore, a truncation of the computation already at one-loop order omits important thermal mass contributions at two-loop order. In turn, this causes a large leftover renormalisation group (RG) scale dependence, see ref. [171].

Perturbative 3d approach
Also this method still suffers from the IR problem of perturbation theory. We emphasise that for the perturbative effective potential itself, there is no real quantitative difference between 4d and 3d approaches, provided that in both cases the computation is performed to the same order in both coupling expansion and high-T expansion. However, dimensional reduction systematically accesses higher order resummations and it is customary to include a consistent O(g 4 ) accuracy by a two-loop level computation, which yields a reduced RG scale dependence. A gauge dependence of the analysis can still be a theoretical blemish, but a gauge invariant treatment is possible by employing a -expansion and computing gauge invariant condensates [103,157]. Although radiatively induced transitions suffer from IR divergences at O( 2 ) [215].
As a downside, the perturbative 3d approach is harder to automate and streamline compared to 4d approach due to additional steps. Although the computation of the 3d effective potential of a 3d EFT can be related to known topologies arising at two-loop order, the automation of dimensional reduction and the construction of the 3d EFT are still not common standard. For developments, see [157,169,221,222]. In the future, automated dimensional reduction for multiple BSM theories could permit a perturbative 3d approach to be implemented to software that currently relies on the 4d approach.

Non-perturbative 3d approach
This method solves the IR problem, by treating perturbative (hard and soft) modes perturbatively while non-perturbative ultrasoft modes are analysed by lattice simulations. Furthermore, lattice simulations provide manifestly gauge invariant results. While this approach is very technical and computationally slow and demanding, it is still straightforward compared to direct 4d simulations (see ref. [173] and references therein). Recent attempts [152,156,158] simulate phase transitions in a limited number of BSM setups and benchmark points. The hope is to expose general trends regarding accuracy and reliability of simpler tools in perturbation theory. However, model-independent or conclusive results are unavailable so far and similar investigations are actively continued in the future. Finally, we highlight that the 3d EFT approach is also an applicable and attractive framework for non-equilibrium physics of phase transition, such as bubble nucleation and sphaleron rate; see refs. [198,199,223].

Dimensional reduction with a real scalar: A tutorial
The following tutorial constructs the dimensionally reduced 3d EFT of a single real scalar field. 4 The machinery presented builds upon classic literature [109][110][111]167] and generalises straightforwardly to more complicated BSM theories with non-minimal Higgs sector. To guide upcoming generalisations of complicated models, we detail step-by-step derivations that can be used for future crosschecks. Furthermore, based on advances of automation in thermal field theories [157,169,221,222], we implemented in-house software in FORM [224] and applied qgraph [225] for diagram generation. Integration-by-parts reductions (IBP) follow a standard Laporta algorithm [226] adapted for thermal integrals [221]. This algorithmic 4 An independent computation [158] treats masses and tadpole as interactions in strict perturbation theory and in the unbroken phase. The dictionary {µ1 ↔ σ, µ 2 σ ↔ m 2 , µ3 ↔ g 2 , λσ ↔ λ 6 } gives a direct comparison, wherein g denotes the cubic coupling of the real scalar not to be confused with a formal power counting parameter or the SU(2) gauge coupling.
perturbative treatment fully automates the computation of correlation functions within the unbroken phase and their matching. This software can tame the ever increasing complexity of computations in future models with multiple interacting new BSM fields.
The real singlet scalar model demonstrates all details of dimensional reduction. Without coupling it to the SM with the Higgs doublet, gauge fields, and fermions, this model poses an ideal starting point. The following computations employ explicit resummation to cancel delicate soft/hard mixing contributions at two-loop order. Practically, these IR contributions are trivially dropped [158] in strict perturbation theory. As an instructive crosscheck, we perform the computation both in the broken phase using the effective potential, the generator of correlators, as well as the unbroken phase computing correlators directly.
Section 4 focusses on full xSM -where a real singlet scalar is coupled to the SM Higgsand presents the definition of the EFT including results. Details of this full computation are relegated to appendix A and results of appearing (sum-)integrals are collected in appendix B. For the derivation of such integrals, we refer refs. [221,227,228]. Our notation follows ref. [149] in which dimensional reduction for the xSM was initially discussed. This reference deferred the case of a light ("soft") singlet which remains dynamical in the 3d EFT to our computation.

Model and parameter matching
Consider the theory of a single real scalar field σ, given by bare 4d Lagrangian (with Euclidean metric) in the imaginary time formalism Definitions of bare quantities in terms of their renormalised versions and counterterms in renormalised perturbation theory are found in sec. 2.1 of ref. [149]. Note that we choose a general renormalisable theory with linear and cubic terms without Z 2 -symmetry σ → −σ.
Conveniently, in the simple case of a real scalar (without gauge fields) the 3d EFT and 4d parent theory bear the same form. With the exception that after dimensional reduction the couplings and field live in a spatial 3d theory. We organise our perturbative expansion by establishing the following formal power counting where g is a formal power counting parameter that corresponds to the weak coupling at zero temperature. Once this theory couples to the SM in sec. 4, the formal power counting parameter g is identified as the SU(2) gauge coupling. Within the power counting (3.2), we aim for a dimensional reduction at NLO with O(g 4 ) accuracy. At loop-level, this requires one-loop accuracy for cubic and quartic couplings, and two-loop order for tadpole and mass parameter. Figure 3: Illustration of a O(g 4 ) NLO matching of correlators between 3d and 4d theories with Z 2 -symmetry. A full non-Z 2 -symmetric case is analogous. Blobs present the sum of hard contributions to one-and two-loop diagrams in perturbation theory, and the differentiation (prime) acts upon the external soft momentum K = (0, k).
We emphasise that the above formal choice for the scaling of the cubic coupling leads to a peculiarity illustrated along the tree-level quartic interactions induced by the cubic coupling. Its contribution parametrically dominates over the corresponding quartic coupling, and could even compromise perturbativity at zero temperature. In practice for dimensional reduction, this causes no complications as the above interaction is 1-particle reducible. Hence, it is absent in Green's functions that are matched during dimensional reduction for the ultrasoft (light) field in 3d EFT. By enforcing a different scaling, namely µ 3 ∼ g 2 T , the contribution (3.3) formally scales as the corresponding quartic coupling. However, this suppresses almost all contributions of µ 3 in the matching relations at O(g 4 ).
Hence, our strategy is the following: for generality we indeed install a scaling of µ 3 ∼ gT and include all contributions of cubic couplings in our matching relations. These contributions can always be trivially dropped if an extra suppression is assumed. This choice allows us to more widely illustrate different aspects of the dimensional reduction procedure, such as effects from field normalisation. Indeed, a non-Z 2 -symmetric theory can demonstrate the high-temperature screening of the fields by the hard scale via ring topology diagrams such as in eq. (3.7). These are absent in a Z 2 -symmetric case. The motivation to include these contributions is the presence of similar diagrams in theories with gauge fields and fermions, even if the scalar sector is Z 2 -symmetric. Figure 3 illustrates the NLO matching of the parameters. For references with explicit matching examples, see refs. [109-111, 149, 157, 158]. Loop corrections from the 3d side and soft contributions match exactly. Hence both drop out trivially in the matching, leaving only hard contributions. At two-loop order mixed soft/hard terms are cancelled by countertermlike resummation interactions at one-loop [109,155]. We demonstrate this in sec. 3.2.3.
For example, the matching of correlation functions for the quartic self-interaction in both theories yields Note that here the correlators equal minus the sum of the tree-level vertex and Feynman diagrams at higher orders. For a consistent matching, the resummation of parameters (see details in eq. (3.16)) for zero modes allows to identify 3d loop corrections with soft terms in the 4d computation, and these two (IR) terms cancel. Furthermore, resummation ensures a cancellation of all mixed soft/hard mode contributions in one-and two-point correlation functions at two-loop order. For details, see eqs.
and relating 3d and 4d fields by the first line of fig. 3 leads to a O(g 4 ) result Matching relations of other effective parameters are obtained analogously. The 4d and 3d fields at one-loop level are related by a computation of the ring topology diagram at non-zero external static momentum K = (0, k) with soft |k| = k ∼ gT . Denoting sum-integrals according to appendix B, a series expansion to quadratic order yields where we interchangeably denote Π n = σ n as correlation functions and utilise the intact d-dimensional rotational symmetry (k · p) 2 → 1 d k 2 p 2 , trivial manipulation p 2 = P 2 − P 2 0 and integration-by-parts relations Σ P The effective potential generates all correlation functions at zero external momenta. In cases where an explicit momentum dependence exceeds the accuracy of the computation, accessing correlators simplifies greatly by starting from the effective potential. Shifting the scalar field σ → σ + s with a real background field s, the effective potential reads where the ellipsis truncates potential higher dimensional correlators that lead to O(g 5 ) marginal operators in the EFT.
The complete matching relations that define the dimensionally reduced 3d EFT, encode the thermodynamics of the original 4d theory of eq. (3.1) at O(g 4 ). They read where we indicated the required loop order for a O(g 4 ) accuracy. The coefficients V merely contain hard contributions and correspond to correlators via eq. (3.8).

Computation of correlators
The computation of n-point correlation functions is expounded in two different ways. In the following, we discuss their merits. In the broken phase, the computation uses the mass eigenstate basis and we employ the effective potential which is the generator of the correlation functions. Therein, the scalar field is shifted by a classical background field. In the unbroken phase, in the gauge eigenstate basis, we compute all correlators directly diagram-by-diagram. These two approaches give rise to an equivalent result. Technically, for a real scalar theory the broken and unbroken phase computations vary marginally. However, subtleties of these two approaches become more prominent for gauge field theories with (multiple) scalars in different representations of the underlying gauge symmetry group.

Broken phase: Correlators from the two-loop effective potential
Within the broken phase computation, the Feynman rules for vertices read (3.14) Denoting the four-momentum by P = (P 0 , p), where P 0 = 2πnT for each bosonic Matsubara mode, the free scalar propagator is and employs the notationδ ,0 denotes the Kronecker delta for vanishing zero mode. The broken phase mass parameter m 2 = µ 2 σ + 2sµ 3 + 3λ σ s 2 appearing in the propagator corresponds to the squared mass eigenvalue of field s. The resummation of the zero mode σ 0 writes 5 where Π n contain hard mode corrections. Terms with plus signs resum the zero mode mass and terms with minus sign act as interactions. In particular, we have a quadratic resummation interaction for the zero modes. We also have a UV counterterm interaction for all modes where δm 2 = δµ 2 σ + 2sδµ 3 + 3δλ σ s 2 and δZ σ = 0 at one-loop level. After resummation the propagator reads Perturbation theory is organised order-by-order. Thus, thermal corrections Π n at one-loop are needed explicitly for resummation at two-loop level. Therefore, we already quote the result 5 This (order-by-order) resummation identifies IR contributions by relating soft 4d loop contributions with 3d ones. The soft/hard mixing terms cancel explicitly at two-loop order. For gauge fields this procedure becomes technically complicated (cf. ref. [155]) and ref. [109] indicates that such an explicit resummation is somewhat cosmetical. While IR contributions in the matching must be identified, their specific expressions are obsolete. Therefore, soft/hard mixing terms never appear in the matching of strict perturbation theory [158], following refs. [110,111]. The effective potential including two-loop level reads and even though counterterm and resummation diagrams are one-loop topologies they contribute at equal order as two-loop topologies. Figure 4 illustrates the corresponding two-loop level diagrams. The separate terms in the potential yield where m 2 3d corresponds to the mass eigenvalue in the 3d theory. Therein, all master integrals are defined in appendix B in the high-T expansion and in dimensional regularisation utilising the MS-scheme. On the UV side, all T 2 -independent 1/ǫ and 1/ǫ 2 poles cancel in dimensional regularisation. On the IR-sensitive side, non-analytic, mixed soft/hard terms ∝ m 2 3d cancel due to resummation.
Expanding the effective potential in ǫ and the background field s (cf. eq. (3.8)) results in employing abbreviations for thermal logarithms in which Λ is the 4d renormalisation scale and γ E the Euler-Mascheroni constant. The uncancelled T 2 -dependent divergences correspond to the two-loop 3d counterterms in eqs. (3.59) and (3.60).
One hallmark of the broken phase computation is its diagrammatic simplicity: the combinatorics of permuting external legs is intrinsic in the effective potential. As a drawback one has to evaluate massive sum-integrals at two-loop level to generate the dependence on the background field. Even though reaching O(g 4 ) the mass parameter µ 2 σ itself will not appear within two-loop pieces of the matching relations. This detail facilitates the unbroken phase computation in the next section. As another drawback, in models with multiple scalars, multiple background fields appear and it can be tedious to obtain an analytic series expansion for the effective potential in these background fields. This poses a complication, since an expansion in background fields, analogous to eq. (3.8), is needed to extract the correlators.

Unbroken phase: Correlators from the diagrammatic approach
An alternative approach computes the correlation functions directly diagram-by-diagram (cf. footnote 4). The downside of this approach is its large number of diagrams with several permutations of external legs. Conversely, its extension to more complicated models is conceptually straightforward and even multiple gauge fields and scalars coupling to them can be handled algorithmically. This poses an advantage compared to the aforementioned complications in the broken phase where series expansions in (multiple) background fields were needed. In turn, at two-loop level one can set all propagators massless for the NLO dimensional reduction at O(g 4 ) in analogy to strict perturbation theory (cf. ref. [158]).
Within the unbroken phase computation, the Feynman rules for vertices read Note that the tadpole µ 1 never contributes to 1PI diagrams required for the matching. The scalar propagator reads As mentioned, aiming for O(g 4 ) accuracy allows to treat propagators inside two-loop diagrams as massless. This provides the correct hard mode parts and non-analytic IR sensitive contributions vanish trivially in dimensional regularisation due to an absent mass scale. The 3-point and 4-point correlator consist of the following diagrams including their results in terms of master integrals (cf. appendix B) where we indicated symmetry factors and combinatorial factors related to permutations of external legs. The tadpole (1-point) correlator up to two-loop level yields (c.7) Finally, the diagrammatic expressions of the self-energy (2-point correlator) up to two-loop level read After summing individual diagrams for each correlator, we apply integrals of appendix B and recover the correlators of eqs. (3.28)-(3.31). Recall that the correlator itself is minus the sum of diagrams and within our convention V n = σ n /(n − 1)! given in eq. (3.8).

Cancellation of mixed hard/soft terms
As mentioned earlier, explicit resummation is obsolete since all propagators at two-loop diagrams are treated massless. To this end, we demonstrate how resummation unfolds as a cancellation of IR sensitive mixed hard/soft terms while keeping sum-integrals massive. For simplicity, we discuss the Z 2 -symmetric case and employ the resummation From the quadratic part we can read off the resummed propagator 39) and the resummed vertex for the pure zero modes becomes λ σ → λ σ + Π 4 . Also resummation interaction terms are introduced for the zero mode one observes that masses are expanded in the pure hard terms but are kept in the mixed soft/hard terms. In fact, the resummed 3d mass equals the one-loop dimensionally reduced mass parameter. The last line above reintroduced numerical factors and the scalar selfcoupling. Similarly, the bubble integral (d.4) yields where the pure 3d vertex is resummed and corresponds to the 3d effective one. The mixed mode contributions (A), (B) and (C) are illustrated in fig. 5, together with counterterm interaction diagrams that read .
In practice, the unbroken computation can be conducted by dealing a zero-mass to propagators in two-loop diagrams even though resummation is conceptually indispensable. As a result, IR-divergent contributions vanish in dimensional regularisation which, in turn, obscures the need for explicit resummation.

Matching relations for 3d parameters
The ensuing matching relations define the dimensionally reduced 3d EFT based on the parameters of the fundamental 4d theory and temperature. These follow the explicit computations of the previous sections: indicating contributions originating from field normalisation (f.n.), one-loop, and two-loop level. Note that the high-T expansion gives rise to NLO terms at one-loop which are µ 2 σproportional. Importantly, we explicitly denoted the 4d scale dependence (Λ) in terms of which the running of LO terms produce contributions at NLO (O(g 4 )). In addition, we indicated that the 3d tadpole µ 1,3 and mass parameter µ 2 σ,3 run with the 3d renormalisation scale Λ 3d . An exact dependence on Λ 3d is presented even though it includes higher contributions than O(g 4 ). This exact dependence can be solved due to the super-renormalisability of the 3d EFT; see sec. 3.4.
By applying β-functions of appendix A.1 and ref. [149], we immediately observe that all 3d parameters are independent of the 4d renormalisation scale Λ at O(g 4 ): For example, the temperature-dependent scale dependence of the 3d mass parameter µ 2 σ,3 in eq. (3.48) arises via its one-loop running contribution 1 4 T 2 λ σ (Λ) and cancels upon its twoloop logarithmic term ∝ T 2 L b . As a general feature for other scale-dependent terms, this renormalisation scale dependence is discussed in ref. [171]. It is worth to point out thatas depicted in eqs. (3.47)-(3.50) -the tadpole and mass parameter are running in terms of the 3d renormalisation scale Λ 3d whereas couplings do not. This dependence of 3d RG scale cancels in computations within the EFT, as we illustrate in the next section.

Two-loop effective potential in 3d EFT
To conclude this section, we illustrate the computation of the two-loop thermal effective potential, within dimensionally reduced 3d perturbation theory. This corresponds to step (e) of sec. 2.3. At two-loop level, the 3d effective potential composes of The 3d mass is given by the mass eigenvalue: m 2 3d = µ 2 σ,3 + 2µ 3,3 s 3 + 3λ σ,3 s 2 3 , where s 3 is a background field of the 3d theory. The corresponding vertices read The individual pieces of the potential read with the respective master (loop) integrals collected in appendix B. In the 3d EFT divergences stemming from the field dependence appear only at two-loop level wherefore one-loop diagrams with counterterms contribute at three-loop level. The only parameters in need of renormalisation are the mass, tadpole, and field-independent vacuum counterterm. Their counterterms are The upper two of these counterterms are exact due to super-renormalisability of the 3d EFT; new divergences are absent for tadpole and scalar mass at higher loop orders. The fieldindependent vacuum counterterm gets contributions up to four-loop order [158]. Hence, we can solve an exact renormalisation scale (Λ 3d ) dependence of the 3d parameters by requiring that the bare parameters are scale-invariant: Their solution yields where temperature-dependent initial conditions are parametrised by Λ 0 ≡ 3T e c and coefficients C 1,2 are fixed to reproduce the O(g 4 ) hard mode contributions of the matching relations at Λ 3d = Λ 4d . The above renormalised parameters µ 3,3 in eq. (3.49) and λ σ,3 in eq. (3.50) are RG-invariant. In total, the effective potential at two-loop has a compact result We observe that the tree-level running of the 3d parameters compensates the scale dependence in the two-loop logarithmic terms. This RG-improved effective potential (and even the threeloop effective potential) is compared to non-perturbative lattice simulations in ref. [158]. The latter agrees surprisingly well even with large expansion parameters in perturbation theory. This concludes our instructions to the dimensional reduction of the real scalar field theory.

Dimensional reduction of the real-singlet extended Standard Model
This section details the dimensionally reduced 3d EFT for the SM coupled to a real scalar singlet at NLO. This is the novel result of this article. Applications of this 3d EFT to study of electroweak phase transition in the xSM are presented in refs. [170,171]. When extending the Standard Model by a real scalar singlet, the position of that scalar in the high-temperature hierarchy is a priori undetermined. If the scalar singlet assumes a hard (or "superheavy") scale it is integrated out entirely during the dimensional reduction [149]. 6 The resulting version of the SM 3d EFT encodes effects of the singlet merely in its matching relations.
The following analysis relaxes this assumption and performs the dimensional reduction with a soft (or "heavy") singlet. The singlet remains a dynamical field in the 3d EFT and can eventually become ultrasoft (or "light"). Such a configuration allows for dynamical transitions with two consecutive steps which are a viable candidate for EWPT with SFOPT in this model (cf. singlet refs. in sec. 1). During such a dynamical two-step transition, first the singlet acquires a non-zero vacuum expectation value (vev) at high temperatures which is followed by a SFOPT in Higgs-direction once temperature is lowered further.
The scalar sector of the 4d Lagrangian reads which notationally aligns with ref. [149] (see sec. 2 ibid.) except the opposite sign convention for the 4d Higgs mass parameter µ 2 h . We assume the following formal power counting, or scaling in powers of the SU(2) gauge coupling g: Since both mass parameters are soft, this leads to an EFT with two dynamical light scalars φ 3 and σ 3 . The scaling of dimensionful couplings is more delicate (cf. sec. 2.1.2 in ref. [149]), where for the tadpole and cubic couplings we assume In analogy to eq. (3.3) this formal choice leads to similar peculiarities. Contributions parametrically dominate over the corresponding quartic couplings. We reiterate the strategy below eq. (3.3): for generality we install µ m , µ 3 ∼ gT and include all contributions of cubic couplings in our matching relations. Note, that contributions proportional to µ 2 3 and µ 4 3 are further numerically suppressed by extra powers of 1/(4π). In the chosen formal scaling, the one-loop contributions of the tadpole correlator ∼ µ m,3 × T 2 are formally of the same order as the tree-level tadpole. Therefore, the β-functions of tadpole and mass parameters (cf. appendix A.1), are partly needed at two-loop level for a O(g 4 ) accuracy.

Effective 3d theories
The corresponding effective 3d Lagrangian is with r ∈ {1, ..., d} and the 3-subscript denoting parameters in the dimensionally reduced 3d EFT. 7 In addition, we include the following (non-kinetic) pure scalar marginal operators The nomenclature of the effective field theory [229][230][231] classifies these operators as S 6 . We omit classes with higher dimensional kinetic operators such as D 2 S 4 and D 4 S 2 where D formally presents a derivative operator and classes with gauge fields such as F 3 where F presents a field strength tensor. This choice is purely practical: The matching of class S 6 is straightforward since corresponding correlators can be computed at zero external momenta both from the effective potential in the broken phase and even diagrammatically directly in unbroken phase. On the other hand, the derivative structure of kinetic operators requires a computation with explicit external momenta dependence. We defer this challenge to a future comprehensive analysis of the numerical relevance of different higher dimensional operators. However, in the presence of large portal couplings it is natural to expect the class S 6 to numerically dominate over other classes that are always suppressed by g 2 .
The breaking of Lorentz symmetry by the heat bath induces temporal scalars. These are remnants of the temporal gauge field components [209] and obtain Debye screening masses at the soft scale. In analogy with "electrostatic" QCD [111], the Lagrangian composes of wherein τ a denote the Pauli matrices and the covariant derivatives act on the adjoint scalars as D r A a 0 = ∂ r A a 0 + g 3 ǫ a bc A b r A c 0 and D r C α 0 = ∂ r C α 0 + g s,3 f α βρ C β r C ρ 0 . Several interaction terms among adjoint scalars were omitted in the temporal Lagrangian since they are of secondary interest in our computation [149]. Among these omissions are operators with an odd number of temporal fields such as σφ † A a 0 τ a φ. These only appear in the presence of a finite chemical potential due to the breaking of parity [125]. Exceeding the accuracy of our analysis, we exclude higher dimensional operators involving temporal scalars since they are numerically suppressed compared to large scalar portal couplings [109]. Besides, the effect of the temporal sector is numerically subdominant in strong phase transitions driven by ultrasoft (light) scalar fields. This suppression is often empirically observed in BSM theories since physically the temporal scalars are screened at length scales much shorter than those relevant to the phase transition dynamics.
At the ultrasoft scale, the dynamics of the soft (heavy) temporal scalars A a 0 , B 0 and C α 0 has been integrated out by the second step of dimensional reduction [109]. Remember that their masses are at the soft scale and are not dynamical in the vicinity of the transition. The resulting Lagrangian resembles eqs. (4.6) and (4.7) but all parameters are denoted with a bar.

Integrating out the hard scale
The first step of the dimensional reduction occurs from hard to soft scale by integrating out all hard, non-zero Matsubara modes. Here, the dimensional reduction is performed at NLO (O(g 4 )), which means at one-loop in the couplings, two-loop in the tadpole and masses, and one-loop in the field renormalisations. This chapter merely quotes final results and matching relations, casting details of the computation of correlators to appendix A. We employ a general covariant gauge with gauge parameters ξ 1 for U(1) Y , and ξ 2 for SU(2) opposed to Landau gauge in ref. [149].

Normalisation of fields
The relations between 4d and 3d scalar fields are (for generic field ψ) where primed correlators are differentiated with respect to the external momentum squared. The hat denotes the correlator in renormalised perturbation theory with implicit counterterms δZ ψ in appendix A.1. The corresponding renormalised 2-point correlation functions yield where n f = 3 is the number of quark and lepton families and N c = 3 the number of colours. The U(1) Y hypercharges are for which we abbreviate recurring sums as The matching between the 4d and effective 3d theory relates their couplings at one-loop order whereΓ represents a n-point correlation function with subscript corresponding to external fields. The mass parameters and tadpole at two-loop match according to The matching of marginal operators (4.7) is analogous and we do not explicate their formulas. 8 We showcase the computation of the required correlators in appendix A.4. As an example, for the quartic Higgs coupling, we obtain where we indicated the origin of individual contributions, displaying cancellations between correlators and field normalisation (f.n.) for the RG-scale and gauge dependence. To witness a cancellation of the RG-scale Λ at O(g 4 ) order, we include the one-loop β-function (A.8) for the tree-level piece, which cancels the logarithmic scale dependence in L b,f . In these cancellations, the contribution from field normalisation is essential. Note that the coefficient of ln(Λ/T ) (in L b,f ) matches the β-function for the tree-level parameter. The matching of other parameters yields analogously We observe an important -and bluntly surprising -exception for the Z 2 -symmetric Higgssinglet portal coupling The remaining gauge dependence in λ m,3 originates from the correlatorΓ φ † φσ 2 (cf. appendix A.3) and is uncancelled by the field renormalisation contribution which is proportional to the portal coupling λ m instead. Since the matched parameters are merely 3d effective Lagrangian parameters, they are not directly associated with physical observables of the 3d theory and may well depend on the gauge fixing of the 4d theory. A similar discussion in ref. [222] addresses the role of higher dimensional operators in hot QCD. We defer the topic of gauge dependence in matching of hot electroweak theories to future research. Meanwhile, an immediate solution (used in ref. [170]) is provided by a stricter power counting µ m ∼ g 2 T for the cubic portal coupling. Hence, by sticking to O(g 4 ), all gauge-dependent contributions are cast to higher orders.
The singlet interacts with temporal scalars through the couplings (4.37) The results for the two-loop mass parameters and tadpole coupling read and We explicated the dependence of the 4d RG-scale Λ in terms where the running can be verified to cancel the logarithmic scale dependence, rendering the 3d parameters Λ-independent. The running is determined by the β-functions listed in appendix A.1. In the xSM the cancellation of the RG-scale is slightly more subtle than in the SM, since the running of the masses starts at order O(µ 2 m , µ 2 3 ). Therefore, the running of the one-loop mass corrections is crucial. In the above formulas, the 3d running is in fact not exact, as we did not include two-loop terms with singlet-temporal scalar couplings x 3 , x ′ 3 , y 3 and y ′ 3 . However, this approximation is justified since temporal couplings only arise at O(g 4 ) because the singlet couples only indirectly to gauge fields. Hence, the omitted terms proportional to x 2 3 , . . . are numerically negligible. The SM result for the three-dimensional mass parameter reads ). Therefore, also the one-loop mass correction in the second line of µ 2 h,3 is affected by it which is required for the cancellation of the RG-scale. This result can be found in refs. [109,155] where the latter neglects the two-loop contributions involving g ′ . The matching relations of the marginal operator coefficients are listed in appendix A.3. Additionally also the couplings between temporal scalars and the Higgs doublet receive singlet-induced corrections: Due to the absence of singlet contributions, other 3d parameters in the dimensionally reduced Lagrangian (4.8) agree with the Standard model. For completeness, we collect these results below, with n f = 3 the number of fermion generations We point out, that eq. (3.87) in ref. [149] misprints the gluon Debye mass m ′′ D 2 , which we corrected above in eq. (4.51).
Two-loop electroweak Debye masses m D and m ′ D Above, we quoted the one-loop electroweak Debye masses m D and m ′ D , which are standard, as their two-loop corrections are of higher order at final ultrasoft scale EFT. To reach a consistent O(g 4 ) order at the soft scale, these Debye masses should be computed at two-loop order, and at this order they receive contributions from the singlet. We point out, that the QCD Debye mass is independent of singlet contributions at two-loop order. We omit its two-loop result here, since these are further suppressed at the ultrasoft scale than two-loop contributions to the EW Debye masses because the Higgs couples only indirectly to the gluon sector.
The singlet-induced two-loop corrections to the electroweak Debye masses read

Integrating out the soft scale
The second step of dimensional reduction integrates out heavy temporal scalars at the soft scale. The resulting simplified 3d EFT at the ultrasoft scale assumes light, dynamical doublet and singlet scalars. Since the singlet couples to gauge fields only indirectly, its couplings to temporal scalars A a 0 and B 0 are suppressed already at leading order. Hence, we include only one-loop effects of temporal scalars in ultrasoft singlet parameters instead of including twoloop corrections for tadpole and singlet mass parameter. All correlators are then encoded in the one-loop contribution of the effective potential (see appendix B for the one-loop master integral) Denoting the background fields v 3 for the doublet and s 3 for the singlet, the (3d) background field-dependent mass eigenvalues read and give rise to the (one-loop) matching relations (4.65) Therein soft corrections stem only from correlators since field normalisations contribute at a higher order due to non-existing tree-level contributions. In particular, all soft contributions at leading order are gauge-independent, as there are no gauge field propagators involved at one-loop order. The remaining parameters are referred from [109,157,170] (4.68) (4.69) In general the ultrasoft Higgs self-energyμ 2 h,3 receives contributions from interactions with singlet and temporal scalars. Even though these are two-loop topologies, we discard them due to the suppression of x 3 , x ′ 3 , y 3 , y ′ 3 in analogy with discarding contributions with quartic self-interactions of temporal scalars. This is apparent for κ 3 , κ ′ 3 , κ ′′ 3 that lack a tree-level contribution O(g 2 ) and consequently their leading contribution is O(g 4 ). For simplicity, we drop corrections from temporal scalars to marginal operators due to their numerical insignificance.
These relations complete our construction of the high-T 3d EFT of the SM augmented with a real scalar singlet. As an effective theory, it can be used to examine the thermodynamics of the electroweak phase transition of the fundamental model (cf. refs. [170,171]). In particular, ref. [170] showcases the computation of the two-loop thermal effective potential in the 3d EFT of the xSM constructed in this section. This is analogous to our section 3.4.

Discussion
The pipeline between collider phenomenology of BSM theories and their implications to early universe cosmology and the potential birth of stochastic GW background convolves multiple complicated stages. One goal of this article is to take steps towards that, on theoretical grounds, uncertainties related to the prediction of the thermodynamics are not the largest in this pipeline.
Concretely, we gave a fresh qualitative review of the thermodynamics of the electroweak phase transition and focused on scalar extensions of the SM. In particular, we concentrated on the framework of high-temperature dimensional reduction. As an automatic all-order resummation scheme it perturbatively defers the infrared problem of thermal field theory. Thereafter, the IR sensitive physics is encoded in a dimensionally reduced EFT that can be studied non-perturbatively on the lattice. However, the constructed 3d EFT is powerful already in perturbative studies. The tutorial-styled computations explicated in sec. 3 aim to make this technique more accessible with emphasis on a scientific community studying the thermodynamics of the EWPT for a wide variety of BSM theories.
The majority of studies of the electroweak phase transition in BSM setups are restricted to perturbative computations of the thermal effective potential. They are often limited to oneloop order and naive leading order resummation usually of Arnold-Espinosa type. Therefrom, a straightforward extension to a non-perturbative treatment is less apparent and the IR problem remains in its core. On perturbative grounds, important contributions in the weak coupling expansion are missed which roots in a misalignment of loop and coupling (power counting) expansions. This omission causes a residual, artificial RG scale dependence which cannot be compensated by RG-improvement at one-loop. A recent study [157] concludes that this kind of leftover artificial renormalisation scale dependence can pose a dramatic two to three orders-of-magnitude theoretical uncertainty for subsequent analyses of the cosmic gravitational wave background originating from cosmic phase transitions. Such an uncertainty for thermodynamic parameters can compromise predictions for e.g. the signal-to-noise ratio for LISA and other future GW experiments.
Tools to automate dimensional reduction are much needed to handle large numbers of Feynman diagrams that arise at multi-loop orders. By adopting sophisticated tool from zero temperature, developments towards such automation have been taken recently [157,169]. As a concrete application of dimensional reduction, we derived for the first time the hightemperature 3d EFT of the real-singlet extended Standard Model (with a dynamical singlet) -one of the most widely studied BSM models in particle cosmology. This poses the main technical part of our investigation displayed in sec. 4. Perturbative studies of this 3d EFT to scrutinise the EWPT in this model are presented in refs. [170,171].
We conclude by envisioning specific but also model-independent future avenues: (i) The derived 3d EFT of the xSM is indispensable for subsequent studies. Lattice simulations can probe its equilibrium thermodynamics and in particular expose the character of the phase transition and determine parameter regions that admit SFOPT. Additionally, out-of-equilibrium properties of the phase transition such as bubble nucleation rate can be investigated by non-perturbative studies of the 3d EFT.
A leftover gauge dependence indicates an incomplete basis of higher dimensional operators in the dimensionally reduced theory. In general, it is interesting how their effects influence the IR dynamics of the system [222].
(ii) The real singlet scalar model (not coupled to SM) offers a testing platform for different approaches. Implications could be drawn for dark sector phase transitions, by determining the mapping of 4d parameters and temperature to its 3d phase structure. The latter was comprehensively analysed in ref. [158].
Since a real scalar theory is purely bosonic, it evades problems of discretising chiral fermions on lattice. Hence, a comparison between 3d EFT and full 4d lattice simulations is feasible (see related ref. [148]) also when including higher dimensional operators. In this context, its lattice-continuum relations with higher dimensional operators are needed. It remains to be determined which lattice measurements and extrapolations are required to extract the continuum physics from simulations in presence of higher dimensional operators.
(iii) It would be worth investigating how large parameter space scans of past EWPT studies (using a one-loop thermal potential) are affected when complete O(g 4 ) effects are included. We advocate pre-existing software to implement the perturbative dimensionally reduced 3d EFT approach. An example are parameter space scans using CosmoTransitions [212], BSMPT [211], PhaseTracer [213] to examine the phase structure of individual BSM models.

A. Detailed computation in the xSM
This appendix collects details of the dimensional reduction computation in the xSM from sec. 4 and extends our results of the 3d parameters to marginal operators defined in eq. (4.7).

A.1. Counterterms and β-functions of the 4d theory
One-loop counterterms and β-functions are listed e.g. in sec. 3.2 in ref. [149]. Using field renormalisations Z φ for the Higgs, Z q the left handed quark doublet and Z t the top quark, we define the bare top Yukawa parameter This convention for g Y and its counterterm δg Y align with eqs. (C. 22) and (C.29) in ref. [155] in contrast to eqs. (2.20) and (3.38) of ref. [149]. Since these references use Landau gauge, we merely display the ξ-dependent counterterms in general covariant gauge: where hypercharges are defined in eq. (4.16). Essentially also the unphysical gauge fixing parameter receives renormalisation: The two-loop computation of tadpole and mass parameters receives contributions that require new counterterms: and their corresponding β-functions are 16) which are necessarily gauge-independent. The remaining β-functions are listed in sec. 3.2 of ref. [149].

A.2. Correlators from the one-loop effective potential in general covariant gauge
In the background field method, the scalar fields are shifted by φ i → φ i + δ i,2 v/ √ 2 for i = 1, 2 and σ → σ + s, around real background fields v and s. We can read the scalar correlators from the effective potential expanded in these background fields up to dimension-6 terms. In our convention, the coefficients V relate to correlators appearing in the matching relations of eqs. (4.19)-(4.30) as and similarly for 1-and 2-point correlators and marginal operators. In Landau gauge, the background-dependent mass eigenvalues can be solved from the mass matrix constructed from the coefficients of the bilinear parts of the φ i -and σ-fields that mix in the broken phase. By employing the shorthand notation for the parameters of the shifted theoryμ the scalar mass eigenvalues read where the Goldstone mass eigenvalue m 2 G is triple degenerate. Since the singlet does not couple to the gauge fields or top quark, their mass eigenvalues align with the SM In general covariant gauge, the three Goldstone mass eigenvalues are replaced by [168,232] where m ± 2 2 is double degenerate. Based on these background field dependent mass eigenvalues, the one-loop effective potential becomes And by comparing to the expansion (A.17), we can solve for the desired correlators. For a crosscheck, we also compute all one-loop correlators directly in the unbroken phase.

A.3. Matching of marginal operators
Marginal operators of eq. (4.7) arise at O(g 5 ) and O(g 6 ), at which contributions from the field normalisations are absent in their matching. We get Notably, all coefficients related to operators with Higgs field are gauge dependent, in analogy to the Higgs-singlet portal coupling in eq. (4.35). The class of topologies responsible for this uncancelled contribution can be exemplified by the Higgs-singlet portal interaction and sextic Higgs marginal operator At zero external momenta the first diagram vanishes only in Landau gauge due its transversality, since it is proportional to Σ P P µ P ν D µν (P )/P 8 = ξΣ P 1/P 6 . Where D µν (P ) is the gauge field propagator in covariant gauge from eq. (A.38). Identically, the second diagram contributes to the Higgs 6-point correlator (A.29) where also a leftover gauge dependence remains. Similarly, other marginal operators with external Higgs legs are ξ-dependent. This leftover O(g 6 ) gauge dependence of the sextic correlator was pointed out in ref. [157]. Strikingly, the gauge dependence of the Higgs-singlet portal coupling λ m arises already at O(g 4 ).
This underlines the subtlety of the power counting of the cubic portal coupling µ m as discussed in sec. 4. Higher dimensional operators can be used to estimate the accuracy of the dimensional reduction by adding their effect at tree-level to 3d effective potential. It remains to be understood how a gauge-invariant analysis is to be performed as some of these operators are explicitly gauge-dependent. We leave this endeavour as a future challenge, and note that at this stage these operators can be used as mere numerical estimates of the convergence of perturbation theory.

A.4. Two-loop computation of correlators
This appendix documents diagram-by-diagram the results for the two-loop correlators used in sec. 4.2 in terms of master sum-integrals of appendix B. Despite being an algorithmic loop-diagrammatic exercise, we believe that this explicit documentation can facilitate future endeavours of dimensionally reduced high-T effective theories. Especially as one might find this to be the non-trivial part of a dimensional reduction computation. The computation was performed in general covariant gauge, where gauge parameters enter via gauge field propagators for SU(2) and similarly for other gauge propagators of U(1) and SU (3). The transverse projector is defined as P T µν (P ) ≡ δ µν − P µ P ν /P 2 . This section displays results compactly employing Landau gauge (ξ = 0) where gauge propagators are transversal. This transversality decimates the number of integrals. For the Higgs self-energy, we only list new contributions of the singlet scalar.
The corresponding Feynman rules and conventions (in Landau gauge) are outlined in ref. [149]. The pure singlet diagrams are listed above in sec. 3.2.2, wherefore here we only include diagrams wherein the singlet couples to the SM. Note, that in this case also the pure singlet counterterm diagrams contains SM contributions, and similarly the Higgs counterterms have singlet contributions. Again the correlator is minus sum of Feynman diagrams, and at two-loop level (two-loop diagrams, one-loop counterterm diagrams) we employ massless propagators sufficient for a NLO dimensional reduction; see discussion at the end of sec. 3.2.2.

Electroweak Debye masses and gauge couplings at two-loop
The singlet contributions to the gauge field self-energies are displayed in eqs. (4.55) and (4.56). Diagrammatically the SU(2) gauge field self-energy composes of which is identical for U(1) when replacing the external legs: Their corresponding SM contributions to the Debye mass align with refs. [127,128]: (A.113) Additionally, we show the two-loop singlet contributions to the gauge couplings: These contributions are formally of higher order, i.e. O(g 6 ) in our power counting.

B. Collection of integrals
This appendix collects definitions and results of sum-integrals encountered in our computation. Ref. [228] and references therein further showcase many explicit derivations. We use dimensional regularisation in D = d + 1 = 4 − 2ǫ dimensions in the MS-scheme with renormalisation scales Λ in 4d and Λ 3d in 3d. Euclidean four-momenta are denoted as P ≡ (ω n , p) with the bosonic Matsubara frequency ω n = 2πnT . We define the d-dimensional integral measure as where a primed integral denotes the absence of a zero mode. For fermionic counterparts, we employ the definition of ref. [149].
In pure 3d, we encounter the following integrals where we define the shorthand notation In the 4d computation, we encounter sum-integrals parameterised by where the Matsubara four-momenta have implicit fermion signature P 2 = (2n+σ i )πT 2 +p 2 with σ i = 0(1) for bosons(fermions). Below we list some of the recurring integrals in d = 3−2ǫ (B.20) Here the power of integration-by-parts reduction (cf. ref. [221] and in particular sec. 3.4 of ref. [169]) is obvious. All massless two-loop sum-integrals reduce to one-loop masters. Similar sum-integral structures are listed in appendix C of ref. [155] and can be used to compute pure SM contributions to the Higgs self energy at two-loop. In broken phase computations, we need massive sum-integrals (cf. ref. [109]), expanded at high-T setting m/T ≪ 1 Note that if one uses resummed propagators, all 3d integrals in these expressions obtain a resummed 3d mass instead of a 4d mass. In D SS and D SSS all mixed soft/hard terms (that are non-analytic in m 2 ) will be cancelled in resummation at O(g 4 ). For the latter integral, we do not write down terms of O(m 3 ) explicitly since these are absent in resummation at O(g 4 ).
Whereas the O(m 4 ) hard contribution is required for the matching of µ 4 3 terms. Indeed, in addition to O(m 2 ) terms we need O(m 4 ) terms due to cubic interaction µ 3 which can be compared eqs. (81) and (88) of the classic ref. [109].

Sunset sum-integrals
Finally, let us inspect the high-T expansion of the massive two-loop sunset sum-integral. See also ref. [233]. The result given in (B.24) shows almost surprising cleanliness in the terms of the mass expansion. In particular, it is very pleasant that the coefficients of even powers of mass are given in terms of full sum-integral structures, as opposed to linear combinations of (massless) mixed soft/hard and hard sum-integral structures. Essentially, this can be shown to take place by carefully expanding separately the mixed modes and hard modes of the original massive sum-integral. The arithmetic challenge arises from a proper treatment of the mixed modes, as a naive mass expansion can only occur in propagators with non-vanishing Matsubara index (hard scale contribution). Hence, the proper order of the mass expansion is chosen to be isolated with an iterative approach. This operation is explicitly described below for the special case with three degenerate masses. This simplifies the book keeping of the computation without any loss of information, due to the obvious symmetry between the three masses.
It is well-motivated to symmetrise the computation as far as possible. Thus, we choose to consider the hard modes by setting each propagator structure of the sunset integral identical, with non-vanishing thermal component: . (B.26) We can expand this in a trivial manner up to O(m 2 ) to find the following non-vanishing contributions: where we took note that all scaleless spatial integrals vanish, and more notably the massless sunset integral S 3 (B.11) vanishes [227,234]. The first integral expression mixes soft and hard modes. This naturally follows from the symmetrisation of the propagators, which essentially adds contributions from the mixed part of the expansion. Hence, we reinsert it in the second structure of interest.
Let us start the expansion towards O(m 2 ) by solely considering the propagators with nonvanishing thermal scale, which yields (B.29) In both of these sum-integral terms it suffices to remove the mass terms inside the integral, as we only wish to find contributions exactly at O(m 2 ). Thus, we find the quadratic coefficient of the expansion as where S 4 is defined by eq. (B.12). We follow a similar procedure to find the coefficient of O(m 4 ), with the slight difference of using the results found above as the iterative subtraction element for the mixed element expansion. This time, the hard expansion yields (B.33) In order to fully extract the mixed contributions, we again expand the suitable propagators and follow-up with a removal of the results of previous orders. This results in three separate computations: where S 5 is defined in eq. (B.13) and S 6 in (B.14).