Combining thermal resummation and gauge invariance for electroweak phase transition

For computing thermodynamics of the electroweak phase transition, we discuss a minimal approach that reconciles both gauge invariance and thermal resummation. Such a minimal setup consists of a two-loop dimensional reduction to three-dimensional effective theory, a one-loop computation of the effective potential and its expansion around the leading-order minima within the effective theory. This approach is tractable and provides formulae for resummation that are arguably no more complicated than those that appear in standard techniques ubiquitous in the literature. In particular, we implement renormalisation group improvement related to the hard thermal scale. Despite its generic nature, we present this approach for the complex singlet extension of the Standard Model which has interesting prospects for high energy collider phenomenology and dark matter predictions. The presented expressions can be used in future studies of phase transition thermodynamics and gravitational wave production.


Introduction
Understanding the thermal history of electroweak symmetry breaking is one of the central endeavours of next generation experiments, with both next generation colliders [1][2][3] and gravitational wave detectors [4] potentially giving definitive answers. Furthermore, if electroweak symmetry breaking occurs via a strong first-order phase transition, it could answer one of the central issues of cosmology -why is there more matter than anti-matter in the present day universe [5,6]. Since the Standard Model (SM) of particle physics predicts a smooth crossover transition [7][8][9][10][11], modifying this expectation necessitates new beyond the Standard Model (BSM) states at the electroweak scale . In the presence of a global symmetry, among such new states could in fact be a dark matter candidate [28,39,[50][51][52], giving a minimal unified explanation for the origin of matter.
Describing the electroweak phase transition (EWPT) perturbatively is a theoretical challenge. In non-Abelian gauge theories at finite temperature, perturbation theory fails completely at high enough orders [53]. Even at lower orders, the effective expansion parameter can be large despite a weakly coupled zero-temperature theory. Physically, this is a consequence of enhancement of light bosonic modes at the infrared (IR), and for a perturbative description to converge requires thermal resummations. The thermal plasma exhibits a class of mass hierarchies at high temperature, and in perturbation theory this can be described by integrating out contributions of the heavy thermal scale to parameters of an effective theory at lower scales [54]. This idea of effective descriptions is systematised by high-temperature dimensional reduction [55,56] to three-dimensional effective field theory (3d EFT) [57,58]. Dimensional reduction allows to by-pass problems of perturbation theory at the infrared using non-perturbative lattice simulations of the 3d EFT [7,59]. Due to the excessive cost and technical effort of such simulations, the virtue of perturbation theory within the 3d EFT [60] is to guide thermodynamic investigations.
Past decade EWPT analyses most often focused on minimising the thermal effective potential that adopts the resummation of "the most dangerous daisy diagrams". At leading order, this is achieved by resumming one-loop thermal corrections to masses of zero Matsubara modes and subsequently computing the one-loop effective potential for the resummed zero modes. The 3d EFT picture translates this to a dimensional reduction at leading order (one-loop in effective masses, tree-level in couplings) and a computation of the one-loop effective potential within the EFT. On the other hand, this description fails to include several important two-loop contributions [54,61]. Without these two-loop contributions, leading renormalisation group improvement is incomplete [54,60,62]. Such an improvement cancels the renormalisation scale between the running of leading-order contributions and explicit logarithmic terms at next-to-leading order (NLO). Indeed, due to the slower convergence of perturbation theory at high temperature, this cancellation requires a two-loop level computa-tion as a subset of leading logarithms only appears at this loop level. Recent work [44,[62][63][64] has stressed the numerical importance of a consistent perturbative computation to account for these effects in BSM scenarios of the electroweak phase transition.
transition. This renders the model an appealing minimal BSM candidate for both dark matter and a potentially strong electroweak phase transition. Recent literature has debated of how to organise the perturbative computation of cxSM thermodynamics. Ref. [28] includes leading daisy resummations, 2 but lacks gauge invariance, while ref. [39] maintains gauge invariance at the cost of resummation (also cf. [79] and a more recent investigation [80]). Similar confusion can in principle arise in any BSM model. Our computation combines the benefits of both previous studies and also improves them by removing their aforementioned shortcomings. Consistent with [62,63], we find significant error arising from omitting two-loop corrections to thermal masses. The takeaway message of this article is that such corrections should not be omitted. We furthermore demonstrate how to implement them in a gauge-invariant manner.
The structure of this article is as follows. Section 2 summarises our approach, introduces the model and its 3d EFT at high temperature, and describes the computation of thermodynamics. In addition, we review the perturbative expansion for the effective potential. In section 3 we numerically demonstrate our approach at a representative benchmark point. Section 4 summarises our findings and discusses their application for potential future parameter space scans. Appendix A collects technical details and explicates formulae of our analysis. Finally, appendix B discusses renormalisation group improvement and an epilogue therein exemplifies the importance of including thermal masses at two-loop order, in a toy model of the SM with an artificially light Higgs.

Thermodynamics: resummation and gauge invariance
This technical section details our computation. In fig. 1, we schematically outline the computation of thermodynamics, adapted from [81]. The novel prescription suggested in this article combines (1) Dimensional reduction to 3d EFT at NLO [58] which determines thermal masses up-to two-loo level.
(2) Computing the effective potential within the 3d EFT [60] at one-loop level and inexpansion [66]. For a radiative barrier, heavy field contributions at one-loop level are resummed to the potential at LO [72,73,82,83].
(3) Determining critical temperatures numerically from a condition that the value of the effective potential at different phases is degenerate [65] as well as determining gaugeinvariant condensates and latent heat [59,69].
For (1), it is crucial that thermal resummation by dimensional reduction is performed at NLO, instead of LO. The NLO corrections are sizable, which is already indicated by a large RG scale dependence of the LO computation [62,63]. Typical computations of the thermal effective potential at one-loop level resum also masses at one-loop. Such resummation matches the accuracy of the LO dimensional reduction, and is hence insufficient to eliminate large RG scale dependence.
In (2), we have chosen practicality over full RG improvement. The computation of the 3d effective potential is straightforward at one-loop, as one only needs background field dependent mass eigenvalues. However, RG improvement related to the 3d EFT RG scale can only be eliminated at two-loop order [60]. Naturally, such a computation is more challenging (cf. [63,66,69,81]). However, we demonstrate that in practice the dependence on an uneliminated 3d RG scale is less severe as the aforementioned dependence on 4d RG scale -a trend demonstrated already in [63].
In (3), our numerical implementation follows [65,69] and differs from the strategy in [63,66]. There, the critical temperature itself is expanded in and solved order-by-order.
The remainder of this section details our conventions for the complex-singlet extension of the Standard Model (cxSM), the structure of the corresponding 3d EFT, a review of the perturbative expansion, the thermal effective potential within the EFT and a gaugeinvariant computation of the thermodynamic quantities of interest. Several formulae are collected in appendix A, for a reader to replicate our analysis. While our discussion and concrete expressions focus on the cxSM, we present the underlying ideas generically, for them to be applied in other BSM theories. In particular, we advocate that our prescription could be implemented in software for analysing the electroweak phase transition, such as CosmoTransitions [84], BSMPT [85,86], and PhaseTracer [87].

A complex-singlet extended Standard Model
We use the model of ref. [39] with the Standard Model augmented by a complex singlet scalar S, abbreviated as the cxSM. The scalar sector of the full 4d Lagrangian in Euclidean space The scalar fields are parametrised as where v 0 and v S 0 are zero-temperature real vacuum-expectations-values (VEV) and G − = (G + ) † . After electroweak symmetry breaking, this model has three scalar eigenstates: two CP-even eigenstates resulting from mixing h and s, and one CP-odd state, A. All complex phases that can mix s and A are assumed to be zero. The SM-like Higgs boson is one of the CP-even fields, and A is the dark matter candidate. The fields G ± and z are the usual Nambu-Goldstone bosons as in the SM. A tree-level analysis of relations between the Lagrangian parameters above (in MS -scheme) and the parameters of the mass eigenstate basis can be found in [39] and we merely collect the relevant relations in appendix A.1. We follow conventions of [81,88] for the SM field content. The one-loop renormalisation group equations for the MS -parameters are collected in appendix A.2. Below, we compute selected thermodynamic properties in terms of MS -parameters of the Lagrangian. Only when turning to numerics, we express these parameters in terms of physical quantities.

High temperature 3d EFT
At sufficiently high temperatures, the equilibrium thermodynamics of a weakly interacting quantum field theory (with weak coupling constant, g) in the Matsubara formalism can be described by a dimensionally reduced effective theory (cf. [57,58]). Therein the threedimensional zero Matsubara modes at the soft/light scale (with masses ∼ gT ) are screened by the hard/heavy non-zero Matsubara modes (with masses ∼ πT ). Technically, we shall work at the ultrasoft scale (g 2 T ) EFT. This means that in addition to the hard non-zero Matsubara modes, also the soft temporal scalar fields have been integrated out 3 and only the ultrasoft fields remain. These include fields driving the transition (doublet and singlet) and spatial gauge fields.
The dimensionally reduced high-temperature effective theory for the cxSM is parametrically of the same form as its 4d parent theory in eq. (2.1). We do not restate the scalar Lagrangian for the doublet and singlet fields that defines the 3d EFT, but merely label its parameters with the subscript 3. Quartic couplings in the 3d EFT have dimension of mass (T ), and fields have dimension T 1 2 . We present an effective theory valid at next-to-leading order in dimensional reduction, namely one-loop in couplings and fields and two-loop for masses (except one-loop for parameter b 1,3 ). This corresponds to the formal power counting that is accurate to O(g 4 ). In the absence of cubic couplings, the tadpole parameter (a 1,3 ) does not receive any loop corrections. Accuracy higher than O(g 4 ), namely O(g 6 ), would require additional higher dimensional operators [58]. The derivation of the dimensional reduction matching relations for generic models -i.e. presenting 3d parameters as a function of temperature and parameters of parent theoryhas been demonstrated for example in [58] in particular for the SM and for a real singlet in a recent tutorial [81] (cf. also [64,89]). By extending this computation to the case of a complex singlet, we derive the matching relations presented in appendix A.4. We have explicitly verified their gauge invariance: matching relations at O(g 4 ) are independent of the gauge fixing choice of the parent 4d theory. We emphasise, that in general dimensional reduction at NLO (and at LO) is a gauge-invariant by construction; see [58,60,63,73,81,90]. For the past two decades, the popularity of the 3d EFT approach for the EWPT phase transition thermodynamics has been superseded by that of lower order computations. One presumable bottleneck denying more applications, is the inability to perform dimensional reduction at NLO. In this regard, we advocate upcoming software [91] that helps to avoid biting the bullet and automates such computations for user-defined models.
To capture non-perturbative IR effects related to the ultrasoft scale, lattice Monte Carlo simulations of the 3d EFT [7,60,89] are required. For the cxSM, these simulations are out of scope of this work. The remainder of the article, concentrates on the perturbative computation of thermodynamics. While inferior, and crucially lacking even a qualitatively correct description for certain cases, 4 these perturbative studies are computationally much less expensive and can be used as a first approximation and valuable guidance for future simulations. Nonetheless, the 3d EFT mapping presented in appendix A.4 is one of the key ingredients for such simulations. In addition, one needs lattice-continuum relations (cf. [92,93]), simulation code including proper Monte Carlo update algorithms (for a recent application, see [89] and references therein), and a computation cluster.

Thermal effective potential in perturbation theory
The central quantity for equilibrium thermodynamics is the free energy of the system, or the effective potential for homogeneous field configurations. We begin by recalling the perturbative expansion of this quantity, along the lines of [62] (also cf. [60], and further e.g. [70,[94][95][96][97] in the context of the QCD and electroweak pressure).

Perturbative expansion
At zero temperature, the effective potential admits a formal expansion in g 2 , resulting in 5 We indicate the loop order at which each contribution arises. In this case, the loop expansion aligns with the expansion in g 2 . Furthermore, we emphasised the leading dependence on the renormalisation scale µ (in dimensional regularisation): at tree-level an explicit dependence on this scale does not occur, but couplings are implicit functions of µ. This scale dependence, or running, is governed by renormalisation group equations, or beta functions µ d dµ g 2 ≡ β(g 2 ). At one-loop order, these beta functions match the coefficients of the logarithmic terms in This is the renormalisation group (RG) improvement at one-loop level. It manifests when couplings are solved from their one-loop beta functions. At higher loop levels, similar cancellations occur when including higher-order beta functions and logarithmic terms. The key feature of the perturbative expansion at zero temperature is that the full O(g 4 ) accuracy, with corresponding RG improvement, can be achieved by a mere one-loop computation. Furthermore, an error made by truncating at one-loop is parametrically of O(g 6 ).
The situation is more challenging in high-temperature perturbation theory. Due to the enhancement of IR bosonic modes and subsequent resummations in perturbation theory, the formal expansion of the effective potential reads The potential is parametrically slower convergent than at zero temperature as an expansion is in g instead of g 2 . Hence, to achieve the same accuracy as at zero-T requires the computation of higher loop levels. In more detail LO: g 2 -terms arise both at tree-level and from one-loop thermal corrections to masses.
NLO: g 3 -terms arise solely at one-loop level in the 3d EFT for soft contributions.
NNLO: g 4 -terms comprise of: 6 (1) one-loop hard contributions in quartic terms, (2) one-loop field renormalisation contributions, (3) two-loop hard contributions to masses and one-loop mass corrections in the high temperature expansion, (4) two-loop soft contributions within the 3d EFT. N 3 LO: g 5 -terms appear at three-loop level in 3d EFT [98].
The NLO one-loop term matches usual "daisy resummation" when only the mass parameters are resummed by hard thermal corrections. Practically, in the 3d EFT approach also higher order resummations are automatically included since also couplings are resummed and masses include two-loop corrections. Truncated terms at O(g 6 ) include two-and three-loop hard contributions, contributions to higher dimensional operators in the 3d EFT, and four-loop soft contributions within 3d EFT. 7 Considering the above breakdown, typical one-loop approximations of the thermal effective potential are lacking in accuracy (cf. e.g. [18,24,34,[99][100][101][102][103][104][105][106]). They are only fully correct at O(g 3 ), and while they do include a subset of g 4 -contributions related to one-loop potential (both zero-temperature and hard mode contributions), the computation is incomplete at that order. It has an O(g 4 ) parametric error. Such an inaccuracy is much worse than the O(g 6 ) achieved at zero temperature at one-loop level. In particular, important logarithmic terms are missing at O(g 4 ) (such as those related to two-loop thermal mass) and this causes a large residual RG scale dependence. Using this feature one can probe intrinsic, theoretical uncertainties in one-loop analyses of the phase transition thermodynamics. For some scenarios it was reported to be alarmingly large [44,62,63,107].
A full O(g 4 ) accuracy is achievable in the 3d approach with NLO dimensional reduction (two-loop for thermal masses) and within the 3d EFT perturbation theory at two-loop level; cf. [44, 58, 60, 62-64, 69, 81, 108, 109]. In this case, the error is of O(g 5 ) which is still parametrically worse than one-loop level at zero temperature. Only by including the three-loop effective potential within the 3d EFT, one can reach O(g 6 ) accuracy -the same as the oneloop accuracy at zero temperature. For a real scalar field this has been computed in [98], while for the SM or its BSM extensions this computation remains outstanding.

Radiatively generated barrier
Before focusing on a concrete computation of the effective potential, we comment on phase transitions with a radiatively generated barrier. Since there a barrier is absent at the treelevel potential (in 3d perturbation theory), the barrier is provided by one-loop corrections. Schematically, the one-loop potential for a generic 3d field Φ 3 with background field φ 3 is of the form where the background field dependent mass eigenvalue m 2 (φ 3 ) ≃ µ 2 3 + 3λ 3 φ 2 3 is related to the field Φ 3 itself, and the mass eigenvalue M 2 is related to another, soft/heavy field with M ∼ gT . The other field could be a gauge field with M 2 ≃ (g 3 φ 3 ) 2 or a second scalar with mass parameter ν 2 3 ∼ (gT ) 2 and portal coupling a 2, For there to be a first-order phase transition with different minima separated by a barrier at leading order, the cubic term has to be of same order as quadratic and quartic terms. This can be achieved in the regime that admits a power counting µ 2 3 ∼ g 3 T 2 and λ 3 ∼ g 3 T , resulting in [71,73] Note, that the term m 3 ∼ g 4.5 T 3 is of higher order and does not appear in the leading-order potential. The formal perturbative expansion of the effective potential for a radiative barrier reads [82,83] V 3d eff ≃ α 3 g 3 + α 4 g 4 + α 4.5 g 4.5 + α 5 g 5 + O(g 5.5 ) . (2.9) In this case, convergence is even slower than in eq. (2.6) and the expansion is formally in √ g (instead of g). In fact, there are two expansions: one related to the heavy field, where each higher contribution is suppressed by g compared to previous order, and one related to the light field where each higher contribution is suppressed by g 1.5 . The latter expansion is responsible for non-integer powers in eq. (2.9). Due to this structure of the expansion, the g 3.5 -term is absent. The NLO g 4 -term arises from two-loop contributions of heavy fields. The NNLO g 4.5 -term arises at one-loop order for this light field, as we have seen above. The N 3 LO g 5 -terms are sourced from three-loop diagrams of the heavy field.
The work at hand performs computations in the 3d EFT merely to one-loop order. 8 Consequently, we can only access the leading-order potential for radiative transitions.

Computation of the effective potential
Next, we compute the effective potential to one-loop order within the 3d EFT [60]. Doublet and singlet fields can be shifted by the real background fields v 3 / √ 2 and s 3 / √ 2, 9 respectively. As a result the potential takes the form which introduced as a formal loop counting parameter within the 3d EFT. Tree-level and one-loop contributions read where the one-loop master integral (with MS scheme dimensional regularisation) in three dimensions has a simple result (2π) d and the last equality holds for d = 3 − 2ǫ dimensions. Here µ 3 is the renormalisation scale of the 3d EFT and γ is the Euler-Mascheroni constant. The background field dependent mass eigenvalues are functions of 3d parameters and are collected in appendix A.3. In the above one-loop expression, we used general covariant, or Fermi, gauge with gauge parameters ξ 2 for SU(2) and ξ 1 for U(1) fields. Gauge dependence appears only in the Goldstone mass eigenvalues m 1,3,± in eq. (A.22) and m 2,3,± in eq. (A.23).
The tree-level term captures the hard mode contributions at O(g 2 ) and O(g 4 ). The oneloop term matches the conventional daisy resummed cubic terms at O(g 3 ), while furthermore including a subset of higher order resummations: all 3d parameters in V 1-loop 3d are resummed at O(g 4 ), while in typical LO daisy resummation only mass parameters are resummed, at O(g 2 ). Importantly, in the 3d effective potential, via two-loop matching of the mass parameters, we include the thermal masses and hence quadratic terms at two-loop order. In contrast, typical direct computations of the thermal effective potential only include these terms at oneloop. Including two-loop thermal masses improves the 3d EFT based approach by reducing renormalisation scale dependence, as was pointed out in [62] and further demonstrated in sec. 3. For readers unfamiliar with the 3d EFT approach to resummation of the thermal effective potential, see appendix A in [62] for a comparison to a typical thermal effective potential computed directly in 4d parent theories.
Finally, we inspect the LO effective potential for a singlet field that undergoes the transition in the presence of parametrically heavier Higgs and gauge fields. In practice, this happens for the first step of a two-step phase transition. By construction, the heavy Higgs field does not acquire a non-zero background field. For simplicity, we assume here a Z 2 -symmetric case, such that the singlet tadpole does not enter. The result reads where "rb" stands for radiative barrier and m 2 φ,3 ≡ µ 2 h,3 + 1 4 δ 2,3 s 2 3 is the 4-degenerate doublet mass eigenvalue in the presence of a non-zero singlet background field. The last term is the one-loop contribution from the heavy Higgs, that is formally O(g 3 ). Note that gauge fields do not contribute at LO since they do not couple to the singlet. This result for the effective potential can be interpreted to correspond to an alternative EFT where the heavy Higgs field has been integrated out [72,73].

Thermodynamics
The pressure encodes the information of equilibrium thermodynamics. It is related to the effective potential as p = −V 4d eff = −T V 3d eff and in particular the pressure differences between different phases that are described by the minima of the effective potential. At the critical temperature, the pressures of two phases are equal ∆p = 0 which translates to the condition of degenerate minima in the effective potential ∆V 3d eff = 0. Here and below, we denote ∆(. . .) ≡ (. . .) low − (. . .) high for the difference between the low and high temperature phases. Hence, we do not consider those contributions to the pressure that are present when both background fields vanish [70,96,97], since these are equal in the symmetric and broken phase.
A sketch of the effective potential for a two-step phase transition [49,69,110,111] is presented in fig. 2. It focuses on the second transition from singlet to Higgs phase. At higher temperatures, the extremum at the origin becomes the global minimum and the system is in Figure 2: Schematic illustration for two-step phase transition, showcasing the effective potential near the critical temperature (T c,φ ) of the second phase transition from singlet to Higgs phase. At T > T c,φ (left), the minimum -depicted with a black dot -in singlet direction is global. At T < T c,φ (right), the Higgs phase becomes energetically favorable and the system undergoes a second phase transition. Since the singlet and Higgs minima are separated by a barrier, the transition is of first order. This barrier exists already at tree-level and therefore the second transition can be strong and relevant for gravitational wave production. the symmetric phase. In a first-order transition, the scalar condensates defined as [59] act in analogy to order parameters since they can be discontinuous at the critical temperature. The factor two in eq. (2.16) is a consequence of the chosen normalisation in eq. (2.1) for the singlet mass term. To measure the strength of the phase transition, we compute the latent heat L = T ∆p ′ . Here, the prime denotes a temperature derivative. In terms of the effective potential it reads In a naive perturbative treatment, the effective potential is simply minimised at different temperatures to find the phases, to determine the critical temperature and strength of the transition, as described above. However, this treatment leads to a subtlety related to gauge invariance: the value of the effective potential at its extrema, and therefore also at the minima, are gauge invariant. However, the minima themselves, i.e. values of the background fields, are gauge dependent. Furthermore, outside the extrema, the value of the effective potential is also gauge dependent and a numerical minimisation to find the minima of the potential inherits an artificial, residual gauge dependence of the gauge fixing parameters. Therefore, care has to be invested into this issue, as done in [65,66,71].

Gauge invariant computation
The Nielsen identities [67,68] guarantee that the value of the effective potential in its minima are gauge invariant order-by-order in perturbation theory. In the -expansion [65,66] .7)) and the potential evaluated at the minima expands as The generic form of the O( 2 ) correction can be found in [69], but we do not include it in our analysis as it requires a two-loop computation of V eff 3d . At O( 2 ) there would be additional contributions involving derivatives of the tree-level and one-loop pieces with respect to the background fields, as well as the two-loop potential itself. As we truncate our computation at O( ), the only difference to the effective potential in terms of generic background fields is that we evaluate both tree-level and one-loop parts at the tree-level minima. In Fermi gauge and at one-loop level, the gauge fixing parameters appear solely in Goldstone mass eigenvalues m 1,±,3 and m 2,±,3 (eqs. (A.22) and (A.23)). These vanish at the tree-level minima v 3,0 , s 3,0 which in turn provides a gauge-invariant treatment at O( ).
The procedure described above is improved compared to the PRM-scheme proposed in [65]. It consistently resums hard thermal loops (in 3d EFT parameters) to next-to-leading order (i.e. O(g 4 ) in a formal power counting in g) while maintaining the gauge invariance. As described earlier, this ensures partial RG improvement related to the hard thermal scale with consistent resummation, and reduces the intrinsic uncertainty of the computation [62,63]. In practice, we can find gauge-invariant critical temperatures in analogy to [65]: by determining values of the effective potential of eq. (2.18) in each minimum as function of temperature, and determining when the curves intersect (cf. fig. 3 (left)). On algorithm level [65], this is more efficient than the numerical minimisation of a complicated two-(or multi-) variate function. We also emphasise, that the condensates in eqs. (2.15)-(2.16) are gauge invariant when evaluated using eq. (2.18).
While gauge invariance is manifest in this treatment [66], one subtlety remains: radiatively generated or loop-induced transitions require additional care. In fact, for an expansion in to be meaningful, the leading order minima v 3,0 and/or s 3,0 have to exist. For the Z 2 -symmetric case (without mixing between scalar mass eigenstates) tree-level broken minima exist only if the (3d EFT) mass parameters become negative. In practice [69] this leads to a condition that the broken minima are global minima immediately when these mass parameters change sign. Therefore in this approach, solutions for the critical temperature are determined only from the condition that the corresponding 3d EFT mass parameter vanishes. This is not the physically correct picture. Another reason to be cautious, is that at two-loop order, the effective potential, or its derivatives related to gauge-invariant condensates, have a spurious IR divergence at a vanishing mass parameter which renders the computation further nonpredictive [69]. Similarly, in the original ref. [66], an expression for the critical temperature itself was given in -expansion and found to diverge at O( 2 ). Here, this problem haunts one-step transitions from the symmetric to the Higgs phase, as well as the first step of a twostep scenario from the symmetric to the singlet phase. For the latter, we further demonstrate this issue below for a concrete numerical example.
A cure of this technical problem is the resummation of a subset of one-loop contributions to the leading-order potential. Contributions from the heavy field give rise to a barrier at leading order (cf. eq. (2.14)). Consequently, minima of the leading-order potential are away from the point where the 3d EFT mass parameters vanish and the problematic IR behaviour described above can be avoided. However, as we only work with the leading-order potential for a radiatively generated first step of the transition, we are not able to directly demonstrate if higher order corrections are indeed free from spurious IR divergences. A related and detailed discussion on this topic can be found in [83].

Numerical analysis
To demonstrate numerically the setup for thermodynamics described in the previous section, we investigate a single parameter space point: with m H1 = 125.1 GeV identified with the observed Higgs boson mass. The pole masses of the other scalar bosons, m H 2 and m A , are related to MS parameters in appendix A.1, and the singlet portal coupling and self-interaction coupling are treated as input parameters. This parameter space point matches the S2 scenario studied in [39], and admits a two-step phase transition scenario. In fact, this point belongs to a subset of parameter space for which the tree-level potential eq. (2.1) is Z 2 -symmetric S → −S, viz. a 1 vanishes. In the same parameter space point also both the singlet VEV, v S 0 , and the mass parameter b 1 vanish. This simplifies the expressions in appendix A.1. For the replicability of the analysis, we explicitly write the MS parameters at the initial scale M Z : which are obtained by solving the MS parameters as a function of the input parameters of eq. (3.1) using the relations of appendix A.1. Here, we displayed rounded up numbers while our analysis uses higher decimal accuracy. For comparison, we describe below both the gauge in-and dependent determinations of our numerical analysis. It proceeds in the following steps: 1. For a fixed input parameter space point, solving the MS parameters by tree-level relations (cf. appendix A.1).

Solving RG running of the MS parameters from the one-loop beta-functions (cf. appendix A.2).
3. For fixed temperature T , the parameters of the 3d EFT are obtained from the matching relations of appendix A.4 in terms of the MS parameters that are run to a chosen Tdependent scale.
4. For fixed temperature T , the thermal effective potential within the 3d EFT is constructed as described in the previous section either in the gauge-invariant -expansion at leading order minima, or in Fermi gauge as a function of generic background fields.

5.
Steps 3 and 4 are repeated for different values of T and critical temperatures T c , condensates and latent heat are determined, as described in the previous section. We further demonstrate this in figures below in this section.
Partial RG improvement (related to the hard thermal scale) manifests through steps 2 and 3, when the 4d RG scale cancels at O(g 4 ). Unlike in [39], where running is only implemented to the tree-level part of the effective potential, we use running couplings consistently everywhere. This will induce contributions that are formally of higher order than the accuracy of our computation. In spirit of [2,62,63,[112][113][114], this indicates the intrinsic uncertainty of our analysis. We emphasise that, in a consistent cancellation of 4d RG scale at O(g 4 ), two-loop pieces of 3d mass parameters -or two-loop thermal masses -are essential [62]. 10 This technical detail has been overlooked by all existing literature on the cxSM thermodynamics. Our computation is the first one to account and demonstrate the importance of this partial RG improvement. By working at O( ) or one-loop within the 3d perturbation theory, we do not have full, consistent leading RG improvement at O(g 4 ). However, we can estimate the magnitude of the missing two-loop contributions by including the 3d running of mass parameters and varying the corresponding 3d RG scale.
The gauge invariance of our analysis is guaranteed by the gauge-invariant matching relations of step 3 and the -expansion in step 4. For comparison, we also perform the gaugedependent analysis in the Landau gauge which has been one of the most common choices in the literature for the EWPT in different BSM extensions.

Critical temperature
We first determine the different phases as a function of temperature from the value of the effective potential at different local minima of the leading-order potential. For our benchmark point, the formulae of the tree-level minima admit a simple closed form: In other parameter space points, that are not considered here, solutions for minima can be functionally more complicated. In terms of the above expressions, we denote The left panel of fig. 3 plots these expressions as a function of the temperature (cf. similar fig. 3 in [65]). In this figure, we present the value of the effective potential in units of g 6 3 ; the gauge coupling squared has dimension of mass in the 3d EFT. Further, we used fixed RG scales µ = 1.25πT and µ 3 = g 2 3 . Later in this section, we investigate how results change by varying these arbitrary RG scales.
The intersection points of V sym eff , V S eff , and V φ eff determine critical temperatures. On algorithmic level, finding intersections of these curves is significantly faster than minimising the  (gauge-dependent) effective potential which is a function of two variables. 11 Since V sym eff and V sym eff are numerically indiscernible within the chosen plot ranges, we do not visualiseV sym eff . The global minimum is easily identified with the lowest value and we denote it by a solid line, while dashed lines indicate that each local minimum is no longer the global one. Each of the three minima are global in turn, signalling a two-step phase transition. The intersection points of the global minima reveal the critical temperatures -illustrated by vertical lines -for both phase transitions. As already discussed above, the determination of the critical temperature of the first transition, T c,S , is cumbersome as only after this point (i.e. at lower temperature) the singlet minimum exists (at higher temperature V S eff has an imaginary value), and T c,S is therefore determined by the condition b 2,3 (T c,S ) = 0. This would lead to an IR divergence for the singlet condensate already at two-loop level [66,69] (also cf. [115]), and further signals that interpreting this point as the physical critical temperature is incorrect. The physically correct picture is provided by solving the broken singlet minimums 3,0 from eq. (2.14), i.e. d ds 3 V eff,rb (s 3 ) = 0 , (3.11) at s 3 =s 3,0 . This equation has four different solutions for the broken phase extremas 3,0 = 0 and in the case of interest, the broken minimum is described bȳ For a radiatively generated barrier, the functional forms for the different extrema become seemingly more complicated. However, the minimum at each temperature can still be identified by a derivative test of a single-variable function. We define the notationV S eff ≡ V eff,rb (s 3,0 ) andV sym eff ≡ V eff,rb (0). The critical temperatureT c,S is determined from the temperature value whenV S eff becomes real andV S eff <V sym eff . SinceT c,S is different from T c,S in which µ 2 S,3 = 0, the former solution for the critical temperature does not necessarily lead to spurious IR divergences for condensates or latent heat at higher loop levels.

Gauge invariant condensates
The gauge-invariant condensates are analogous to order parameters such that for a first-order phase transition they are discontinuous at the critical temperature. As an analogy to a typical gauge-dependent analysis in terms of gauge-dependent background fields, we investigate the expressions [64,113] v phys (3.14) The condensates have discontinuities at the critical temperatures, and the Higgs (singlet) condensate is negative outside of the Higgs (singlet) phase. Therefore these expressions are imaginary outside of their respective phase. In fig. 3  and witness that sufficiently belowT c,S these results overlap until a second phase transition. We fixed RG scales at µ = 1.25πT and µ 3 = g 2 3 . The same plot presents the gauge-dependent background fields at the global minima of the effective potential (eq. (2.10)) in Landau gauge ξ 1 = ξ 2 = 0 (dashed lines). From this comparison, we observe that outside of the critical temperature of the first transition, both approaches agree rather well. 12 We interpret this as an echo of common lore in the literature (cf. e.g. [112]) that the Landau gauge results -albeit inherently gauge-dependent, and therefore unphysical -are numerically close to gauge-invariant results, perhaps encouraging their use as a practical guide.
On purely theoretical grounds, gauge-dependent predictions for physical quantities are unsatisfactory and inherently incorrect. Hence, it is invaluable to develop and use a theoretically sound technique to analyse phase transition thermodynamics. The Landau gauge results have been hoped to be useful in practice, with the expectation that the error related to unphysical gauge dependence is smaller than other uncertainties of the problem (cf. discussion on varying renormalisation scales below). However, the gauge fixing parameter is still an arbitrarily valued input. In practice, larger values could lead to larger numerical errors additional to the theoretical blemish.

Latent heat and renormalisation scale dependence
The accuracy of our analysis is investigated by varying the input renormalisation scales, similar to recent studies [62,63]. The left panel of fig. 4 plots T c,φ -the critical temperature of the second transition -as a function of the 4d RG scale µ. Similarly, fig. 4 (right) evaluates the latent heat at T c,φ . We show five different lines: lines A, B and C (black) are based on the gauge-invariant -expansion while lines D and E (red) use a direct numerical minimisation of the effective potential in Landau gauge. Similar to fig. 3 (right), the two methods agree closely. Concretely the lines contain: A,D: one-loop level dimensional reduction, B,E: two-loop level dimensional reduction, C: as B, with varying µ 3 = ( µ πT )g 2 3 .
The difference between the two levels of accuracy in dimensional reduction arises from twoloop contributions to the thermal masses. By comparison, we observe that without two-loop thermal masses, the results are strongly RG scale-dependent. As reviewed in sec. 2, the twoloop logarithmic terms compensate the leading running of the one-loop thermal corrections to the mass parameters in quadratic terms [62]. If these two-loop contributions are omitted, the resulting uncertainty completely overshadows the ambiguity related to the numerical difference of gauge-invariant and gauge-dependent (Landau gauge) approaches. Finally, along the black dashed line C, we also varied the 3d RG scale µ 3 (whereas for other lines it is fixed . This variation signals the importance of two-loop contributions at O( 2 ) within 3d perturbation theory that we omitted in eq. (2.18). Note that mass parameters of the 3d EFT are running in terms of this 3d scale and logarithms of the 3d RG scale only appear at two-loop order. While this effect is significantly smaller than two-loop contributions from the hard thermal scale to the mass parameters, it still causes a sizable uncertainty. This motivates to increase the accuracy to O( 2 ) in future computations.

Discussion
The method presented in this article for phase transition thermodynamics avoids the triune poison of gauge dependence, slow convergence of perturbation theory and intractability. Our investigation follows [62,63,66,69] and implements thermal resummations related to the hard thermal scale of non-zero Matsubara modes using the dimensionally reduced 3d EFT. It further used the gauge-invariant -expansion within 3d perturbation theory to compute the thermodynamic quantities pertinent to electroweak phase transitions. In analogy to previous studies [62,63], we find alarming sensitivity to the renormalisation scale at one-loop level and we identified two-loop thermal masses to be the most crucial contributions to reduce this RG scale dependence. Their importance was previously highlighted in [62]. We note that these findings are in contrast to [80].
Thus, we suggest a minimal setup for perturbative accuracy, that still eliminates most of the undesired scale dependence. It combines NLO dimensional reduction including two-loop thermal masses and a one-loop effective potential in the 3d EFT with a simple expression in terms of the background field dependent mass eigenvalues. The difference to [62,64,69,116] is that the two-loop effective potential in 3d EFT is not included; the RG improvement related to hard thermal scale can still be acquired. 13 The described setup relies on the ability to construct the 3d EFT matching relations by dimensional reduction. Such a setup is, unfortunately, still rare in current BSM EWPT literature. To this end, we expect upcoming automated software [91] -designated for this problem -to increase the applicability of improved studies based on the 3d EFT. We further comment that ref. [29] has attempted to improve resummation of hard thermal loops and in particular to compute thermal masses beyond leading order based on "partial dressing" [117] but without using high-temperature expansion nor dimensional reduction to 3d EFT. Another computation of the two-loop thermal effective potential without high-temperature expansion appears in [113] (cf. also [118]).
The results of our gauge-invariant computation do not differ drastically numerically from conventional Landau gauge analyses. 14 Such a difference is completely overshadowed in a mere one-loop analysis by the uncertainty from the renormalisation scale. This, however, supports the common wisdom that while a gauge-invariant computation is theoretically important, it should not come at the expense of resummation and including relevant terms in the coupling expansion. In part, a similar conclusion was reached in [71]. In this regard, our computation significantly improves the previously suggested PRM scheme [65] by incorporating the required resummations consistently, while maintaining gauge invariance.
Finally, for the cxSM, the questions studied in [39,80] could benefit from the tools presented in this article. In particular, by focusing on what can be concluded about the possible thermal 13 This omission is deliberate in favor of technical simplicity. Computing the effective potential in the broken phase in terms of the background fields is significantly more complicated at two-than at one-loop level. Instead, thermal masses for scalars can be computed in the unbroken phase from hard mode contributions to two-point correlation functions, and such a two-loop computation is standard in dimensional reduction literature. Nonetheless, the two-loop effective potential would allow further RG improvement related to 3d EFT renormalisation scale.
14 Our analysis is performed in one single benchmark point and the comparison concerns only Landau gauge.
Therefore, we refrain from generic statements and acknowledge the theoretical importance of a gauge-invariant computation. In practice, most studies still ignore the issue of gauge dependence.
history of EWSB in this scenario, when both gauge invariance and RG improvement are installed. The presented tools can be used for wide scans of the model parameter space and to analyse physical implications. The latter can shed light on which regions of parameter space admit both a strong first-order phase transitions and dark matter candidates and what are the collider phenomenology signatures of these regions. In addition to the improved equilibrium thermodynamics computation of this article, future scans would benefit from a two-loop effective potential in 3d perturbation theory, as well as from the one-loop improved (zero-temperature) relations between MS parameters and physical input parameters. Both of these improvements are available in the real singlet-extended Standard Model (xSM) collected in appendices A and B of [64] and could be generalised to the cxSM. A further computation of the bubble nucleation rate is required for studying the gravitational wave production in the cxSM. Therefore, one can combine the 3d EFT presented in this article with recent technology [72]. In general, the same approach applies to other BSM theories and in particular those with extended scalar sector.
doublet. Using the shorthand notation g 2 0 ≡ 4 √ 2G F M 2 W , we have To relate the MS parameters to the above input parameters, we use analytic, tree-level relations In many other parameter space points, the relations between MS parameters and input parameters do not have such simple analytic relations for non-vanishing a 1 , v S 0 and the mixing angle α. In these cases, one can solve the MS parameters numerically by inverting the mass eigenvalues (m + is identified with m H 2 and m − with m H 1 ), together with the tadpole conditions .5) and the equation for the mixing angle All above relations hold at tree-level and receive quantum corrections in zero-temperature perturbation theory. In a consistent power counting up to O(g 4 ), one would need to include one-loop corrections to relations of MS and physical parameters (cf. [44,58,64,109,113]) -we do not include them here. Due to this omission we do not have the proper initial conditions in our partial RG-improvement at O(g 4 ). However, our discussion remains qualitatively intact as these initial conditions merely shift the corresponding benchmark point in the MS parameter space. Once this point is carefully related to particle collider phenomenology constraints, the improved initial conditions become quantitatively relevant to be accounted for.

A.2. Renormalisation group equations
By parametrising the MS renormalization scale through t = ln µ 2 the one-loop β-functions, or renormalisation group equations, for the MS parameters read In the beta-functions for µ 2 h and λ h , we depict explicitly only the new complex singlet contributions; pure SM contributions are collected in e.g. [88].

A.3. Mass eigenvalues for background field method
Parametrising doublet and singlet fields in analogy to eq. (2.2), but replacing zero-temperature VEVs, v 0 and v S 0 by generic, real background fields v and s results in mass eigenvalues for the background field method where we used the shorthand notation These mass eigenvalues, in terms of parameters of the 4d parent theory (2.1), can be inverted to obtain the tree-level relations between MS parameters and physical masses in appendix A.1. When computing the effective potentials in the 3d EFT, eqs. (2.10) and (2.14), we denote background field dependent mass eigenvalues by an additional subscript 3. Then all parameters and background fields therein are those of the 3d EFT. In addition, when computing the effective potential in Fermi gauge, Goldstone mass eigenvalues m 2 χ are replaced by [120] where m 2 2,± are double-degenerate. The singlet does not contribute to mass eigenvalues for weak bosons, viz.
and the mass eigenvalue for the photon vanishes.

A.4. Matching relations for dimensional reduction
We define a shorthand notation where ζ s for ℜ (s) > 1 is the Riemann Zeta function. Matching relations between parent 4d theory and 3d EFT read Here, we explicitly suppress pure SM contributions since they can be found e.g. in [81] together with the matching relations for the 3d gauge couplings and parameters of the temporal sector. The tutorial of [81] also computes similar matching relations for a real singlet field. The tadpole coupling a 1,3 does not receive loop corrections since cubic couplings are absent. All matching relations eqs. (A.27)-(A.33) are gauge invariant. By computing the matching in Fermi gauge, we could demonstrate an exact cancellation of the gauge fixing parameters (cf. similarly [63,73,81]).

B. Renormalisation group improvement
At high temperature, the running with respect to the 4d RG scale µ, described by betafunctions β(g 2 ), is not altered compared to zero temperature. However, an overall RG improvement for the effective potential at high temperature is more subtle because thermal scale hierarchies reflect to the structure of the potential [60,62]. In the mapping between 4d and 3d EFT (cf. sec. A.4), the 4d RG scale µ cancels between LO running and explicit logarithms within the L b -terms at NLO in the matching relations of the 3d parameters. This corresponds to cancellation of the RG scale related to hard mode contributions in eq. (2.6).
In the 3d EFT perturbation theory constitutes another, independent, renormalisation scale µ 3 . Without higher dimensional operators that appear at O(g 6 ), the 3d EFT is superrenormalisable. Hence, counterterms in dimensional regularisation have a finite number of terms sufficient to cancel all UV divergences at any loop order. In the SM -and many of its extensions -only the mass parameters (or tadpole parameters) require renormalisation. The 3d couplings are RG invariant and the mass counterterms are exact already at two-loop order. This gives rise to the exact running of mass parameters in terms of µ 3 . Since the running of the mass appears at two-loop order, the cancellation of µ 3 in eq. (2.6) leap frogs over even and odd terms. The running of g 2 -terms cancels logarithms at g 4 -terms, and the running of g 3 -terms cancels logarithms at g 5 -terms and so forth.
In this context, the discussion in ref. [39] on the RG improvement of the cxSM is incomplete. Therein, the one-loop effective potential is divided into a tree-level piece, the zero-temperature Coleman-Weinberg (CW) potential, and the one-loop thermal function. Consequently, the parameters in the tree-level potential are replaced by the parameters solved from one-loop RGEs which correctly eliminates the explicit logarithmic RG scale dependence in the CW potential. Contrary to [39] the renormalization scale does enter the high-T effective potential. While it is true that the one-loop thermal function is not explicitly dependent on the renormalisation scale, there is still an implicit running of the parameters inside the thermal function. The latter, contributes at same order as running of tree-level potential. Crucially, should one implement high temperature expansion to this thermal function, it is straightforward to observe that its leading behaviour of the quadratic terms -that contribute to one-loop thermal masses -is of same order in formal power counting as tree-level terms, i.e. O(g 2 ).
The running of the one-loop thermal mass contributions is therefore an effect of O(g 4 ) and this running is cancelled by logarithmic terms for thermal masses that appear at two-loop order. It is indeed these contributions that we include in our analysis with a NLO dimensional reduction. We also demonstrated their crucial importance in our numerical analysis of sec. 3.

Epilogue. Standard Model with a light Higgs
As an illuminating example, we discuss how even in the pure SM the running of the one-loop thermal mass can cause an alarming leftover scale dependence if the full two-loop thermal mass is not computed. For illustration, we work with a toy model of a pure SM (i.e. we omit complex singlet contributions from the dimensional reduction and 3d EFT of earlier sections) and we vary the Higgs mass in m h = (50 . . . 130) GeV. For simplicity, we determine the critical temperature T c from the condition that the minima of the one-loop effective potential in 3d EFT are degenerate. We use the ratio φ c /T c as an estimate of the transition strength. By doing so, we provide an analogy to many BSM studies that analyse the electroweak phase transition simply in Landau gauge in terms of the gauge-dependent T c and φ c /T c . In fig. 5 we depict these quantities as a function of the Higgs mass.
The leftover 4d scale dependence is illustrated by the resulting bands from varying µ = (0.5 . . . 2.0) × πT and by showing the result with the two-loop (one-loop) thermal mass in dark grey (light grey). A larger band in the latter case signals a larger intrinsic uncertainty, and reinforces our key message of the importance of the two-loop thermal mass. In BSM theories that involve large portal couplings to the Higgs, this theoretical uncertainty related to a leftover scale dependency can be even worse. In the SM, the top quark contributions dominate over those related to gauge couplings and the Higgs self-interaction. We note, that the running of the top quark contribution in the one-loop thermal mass that causes the dominant effect on the broader band, is exactly the same contribution that causes a similar uncertainty for a SMEFT study [63], since this model has exactly the same one-loop thermal mass. The dominant uncertainty in SMEFT roots in the SM sector and is not related to the new higher dimensional sextic Higgs operator.
From fig. 5 (right), we observe that the transition strength increases for lower Higgs mass which corresponds to smaller Higgs self-coupling and a smaller dimensionless 3d quantity x ≡ λ h,3 /g 2 3 . This quantity x is the expansion parameter of 3d perturbation theory [59,60]. The smaller it is, the better 3d perturbation theory works. 15 However, for smaller x the 3d matching relations are also relatively more sensitive to the 4d renormalisation scale, and in the end the total uncertainty is larger. This provides a counter-example for the common folklore which states that strong phase transitions (with large φ c /T c ) are better described in perturbation theory, than weak transitions. On the other hand, fig. 5 makes it evident that perturbation theory is oblivious to the crossover character of the SM phase transition after an end-point around m h ∼ 70 GeV [7]. There is still a discontinuity in critical value φ c (indicated by its non-zero value). In fact, perturbation theory predicts a first-order phase transition and fails even at a qualitative level. This showcases the theoretical challenge to describe phase transition thermodynamics. 15 A similar study for a simpler setup of a real scalar theory [89] found the explicit 3d expansion parameter and inspected convergence and RG improvement within 3d perturbation theory up to three-loop order.