Perturbative effective field theory expansions for cosmological phase transitions

Guided by previous non-perturbative lattice simulations of a two-step electroweak phase transition, we reformulate the perturbative analysis of equilibrium thermodynamics for generic cosmological phase transitions in terms of effective field theory (EFT) expansions. Based on thermal scale hierarchies, we argue that the scale of many interesting phase transitions is in-between the soft and ultrasoft energy scales, which have been the focus of studies utilising high-temperature dimensional reduction. The corresponding EFT expansions provide a handle to control the perturbative expansion, and allow us to avoid spurious infrared divergences, imaginary parts, gauge dependence and renormalisation scale dependence that have plagued previous studies. As a direct application, we present a novel approach to two-step electroweak phase transitions, by constructing separate effective descriptions for two consecutive transitions. Our approach provides simple expressions for effective potentials separately in different phases, a numerically inexpensive method to determine thermodynamics, and significantly improves agreement with the non-perturbative lattice simulations.


Introduction
Gravitational waves from a cosmological phase transition could provide a window to directly observe the very early universe, preceding the birth of the cosmic microwave background as well as Big Bang nucleosynthesis. This would offer a probe of the fundamental constituents of matter and their interactions which is complementary to particle colliders.
In recent years, studies of the electroweak phase transition (EWPT) have sparked a lot of interest, motivated by the possibility of explaining the baryon asymmetry of the universe [1,2], and also of generating a stochastic gravitational wave (GW) background [3] observable by LISA-generation experiments [4][5][6]. In the Standard Model, electroweak symmetry breaking occurs via a smooth crossover [7,8], so a possible first-order electroweak phase transition requires the existence of new physics beyond the Standard Model (BSM). The search for BSM physics that could alter the thermal history of electroweak symmetry breaking provides a target and challenge for future collider experiments [9]. Of particular interest are multi-step phase transitions in the presence of multiple Higgs-like fields, where the transition to the EW phase can be preceded by another phase at a higher temperature [10]. In this work, we discuss a two-step EWPT [11].
In determining the thermodynamic properties of a BSM theory, thermal enhancements of infrared (IR) physics play an important role. In practice, this means that perturbative computations require thermal resummations [12][13][14][15]. In the imaginary-time formalism of high-temperature quantum field theory [16], the most elegant solution to organise these resummations is by means of effective field theory [17][18][19]. In such a computation, a dimensionally reduced effective theory is constructed for the IR sensitive zero Matsubara modes, while all non-zero Matsubara modes are integrated out and their effects are captured in the effective parameters of the EFT. Physically, this accounts for thermal screening, whereby the hard thermal scale modifies the dynamics of the softer IR physics that drives the EWPT.
The most infrared modes of the magnetic gauge bosons become strongly coupled at high temperatures [20], leading to non-perturbative effects on the thermodynamics which require use of Monte-Carlo lattice simulations [21][22][23]. For decades this non-perturbative physics at the ultrasoft scale has caused worry, calling into question the applicability and validity of perturbative determinations of thermodynamics. In the work at hand, we argue that these worries have been somewhat misplaced and argue that in fact first-order phase transitions take place above the ultrasoft scale. Indeed, it has long been known [15,24] that there is a scale in between the soft and ultrasoft scales, the latter of which have been the focus of studies utilising high-temperature dimensional reduction. A proper treatment of physics at this in-between supersoft scale requires a chain of EFTs that we construct in this work. In this approach, these non-perturbative effects are typically subleading for first-order phase transitions, so that both the leading order thermodynamics and several corrections can be obtained by a purely perturbative expansion in powers of a small expansion parameter. More recently, supersoft-scale EFTs have been studied in [25][26][27][28][29][30] (see also [31][32][33][34]).
By direct comparison to previous non-perturbative lattice simulations, we demonstrate that EFTs at the supersoft scale describe thermodynamics with striking accuracy. Further information on the validity of a perturbation approach can be extracted from the expansion itself, by considering renormalisation group (RG) invariance, gauge invariance and the convergence of successive terms. For the former, it has been shown in [35,36] that perturbative computations below two-loop order suffer from large intrinsic uncertainties in terms of sensitivity to RG scales, which reflect the magnitude of missing perturbative corrections. Throughout the paper, we denote RG scales as Λ 4 and Λ 3 , in the full four-dimensional parent theory and in the dimensionally reduced EFT respectively.
The remainder of this article is organised as follows. In Section 1.1 we motivate our analysis by reviewing previous lattice results for a two-step phase transition in the real-triplet extended Standard Model, as well as shortcomings of previous perturbative analyses. In Section 2 we discuss thermal scale hierarchies for generic first-order transitions. Based on these hierarchies, in Section 3 we introduce corresponding effective field theories, paying particular attention to a scale in between the soft and ultrasoft scales, that we dub the supersoft scale. In Section 4 we present what we refer to as strict perturbative expansions for the thermodynamics of phase transitions discussed in the preceding section. In Section 5 we provide a concrete application to the real-triplet extended SM, presenting numerical results while relegating a number of technical details throughout the article to appendices. In Section 6 we summarise and discuss our results.

Motivation
For a two-step electroweak phase transition, the current state-of-the-art determination of the equilibrium thermodynamics is provided by the non-perturbative lattice simulations of Ref. [37], concretely performed for a real-triplet extended SM [38]. In this work, we keep our discussion generic, and applicable to a wide variety of cosmological phase transitions, yet for the numerical analysis turn to the real-triplet extended SM. In this model, the scalar sector comprises the Higgs doublet ϕ and a real triplet scalar Σ a , where a = 1, 2, 3 is an SU(2) adjoint index. We follow the conventions of [37] and define the scalar part of the Lagrangian in 4d Euclidean space as 1 where the definitions for covariant derivatives are standard, and can be found in [39]. This model admits a two-step EW phase transition, where the system undergoes a first phase transition to the triplet phase at some high temperature, after which the system undergoes a second transition to the EW phase. Phase transitions in this model were first studied in perturbation theory in Ref. [38] and Ref. [39] performed the dimensional reduction from the hard to the soft scale 3d EFT for this model. Non-perturbative lattice simulations of the 3d EFT were presented in [37], together with the perturbative computation of the two-loop thermal effective potential. In our numerical analysis, we study the two benchmark points of [37]. These points are defined as  1 We make an exception for the mass parameters µ 2 ϕ and µ 2 Σ for which we use the opposite sign compared to [37,39] [37]. Circular and triangular markers depict lattice results, while solid lines show the two-loop perturbative counterpart following the approach of Ref. [37]. We have added coloured bands which show the RG-scale dependence of the perturbative calculation, as the scale within the EFT is varied over the set Λ 3 ∈ {0.5T, T, 2T }. In this work at hand, we will fix the relatively poor agreement between the perturbative and lattice results apparent here, including the spurious divergence in the perturbative results for the triplet condensate.
where M Σ is the physical triplet pole mass, and a 2 and b 4 are MS parameters at the input renormalisation scale Λ 4 = M Z equal to the Z-boson pole mass. Both of these points exhibit a two-step EW phase transition. According to the lattice study of [37], in BM1 the first (higher temperature) transition to the triplet phase is a crossover, whereas in BM2 the first transition is first order. The second transition is of first order in both benchmark points.
In Ref. [37] a comparison to a state-of-the-art perturbative calculation was provided, utilising the two-loop order effective potential, computed within the dimensionally reduced 3d EFT. The computation utilised the ℏ-expansion of the effective potential [40] -where the potential is perturbatively expanded around its leading order minima, order-by-orderthereby ensuring order-by-order gauge invariance [40][41][42]. However, it suffered from IR divergences related to the determination of the critical temperature of the first transition to the triplet phase. Such divergences were reported already in Ref. [42]: in the ℏ-expansion of the effective potential, the leading order potential does not have a first, but a second-order phase transition, and the critical temperature at leading order is identified with the temperature where the effective mass parameter of the scalar undergoing the transition vanishes. At two-loop order, i.e. O(ℏ 2 ), such a vanishing mass parameter inflicts an IR divergence on the scalar condensate (defined below). This is illustrated in Fig. 1, adapted from Fig. 2 of Ref. [37]. The condensates shown there are equal to derivatives of the free energy density (or equivalently the effective potential) with respect to 3d EFT parameters, . (1.5) These relations follow simply from the path integral definition of the free-energy density, or effective potential [21]. A discontinuity in these condensates signals a first-order phase transition. This is because the first temperature derivative of the pressure p = −T V eff , or free energy density, can be written in terms of the condensates, where ∆ denotes the difference between two phases, the prime denotes a temperature derivative, the η-functions are defined as η(κ i ) ≡ T dκ i /dT , and κ i runs over all 3d EFT parameters. The η-functions depend only on the ultraviolet (UV) thermal scale, and hence are smooth, continuous functions of temperature. It is the condensates that have discontinuities at phase transitions (or kinks for higher-order transitions).
In perturbation theory, jumps in the condensates are related to jumps in the position of the global minimum of the effective potential as a function of temperature. 2 At leading order, the square roots of the quadratic condensates, Eqs. (1.4) and (1.5), agree with the minima of the effective potential. However, this relationship breaks down beyond leading order, and while the condensates are manifestly gauge invariant, this is not the case for the minima of the effective potential [41,43]. A further benefit of the condensates is that they can be computed directly on the lattice as volume-averaged expectation values [21], unlike the minima of the effective potential [44]. As defined above, based on derivatives with respect to MS parameters, the condensates are RG dependent, inheriting their RG dependence from that of the MS parameters themselves. However, this RG dependence is simple and known exactly, due to the superrenormalisability of the 3d EFT. It can be subtracted off to define the following RG invariant combination 8) and likewise for the triplet condensate with β ⟨TrΣ 2 ⟩ = 12g 2 3 /(4π) 2 . From now on we will omit the RG subscript, but when we refer to condensates, we mean the above RG-invariant condensates. Note that these beta functions are exact, and independent of the phase, due to the superrenormalisability of the 3d EFT [21]. On the lattice, this combination is exactly RG invariant, whereas in any finite order perturbative calculation some RG dependence will remain, as a consequence of missing higher loop terms. By varying the RG scale, we gain some estimate of the size of these missing higher loop terms.
In this work, we resolve the failures of perturbation theory visible in Fig. 1 by providing a consistent setup for the radiatively-generated first transition, as inspired by Ref. [30]. In this setup, a physically correct picture of a radiatively-generated transition is provided by consistent power counting, where the barrier separating symmetric and broken phases exists already at leading order, due to integrating out heavy degrees of freedom that generate the barrier. For the second transition, the hierarchies of scale can change, necessitating the construction of separate EFTs for the two transitions, occurring at distinctive thermal scales. This novel construction can be used to avoid the spurious IR divergence of the triplet condensate at symmetry breaking, and provide sound predictions for the critical temperature and strength of the first transition. We find that a complete and gauge-invariant resolution of the failures of perturbation theory also requires perturbatively expanding the critical temperature, following ideas presented in Refs. [33,42]. In order to scrutinize the accuracy of our purely perturbative computation, we compare our results to the non-perturbative lattice simulations of Ref. [37]. In addition, we provide a thorough comparison to some alternative perturbative methodssuch as direct, gauge-dependent minimisation of the effective potential -and discuss their reliability and accuracy, despite their obvious theoretical blemishes.

Thermal scale hierarchies
In weakly coupled quantum field theories, scale hierarchies are a necessary prerequisite for a thermal phase transition. The argument goes as follows. Thermal effects arise through loop diagrams, which are subleading in the vanilla loop expansion. Yet, for there to be a phase transition, these thermal effects must change the effective dynamics of the transitioning field at leading order (LO). This requires a breakdown of the vanilla loop expansion, because the subleading order must match the leading order. 3 Finally, for equilibrium physics, which is time independent and hence absent a light cone, the only kinematic enhancements possible are simple scale hierarchies. In fact, as we will see, at phase transitions in weakly coupled theories there are typically multiple scale hierarchies.
Scale hierarchies wreak havoc with the loop expansion. Any large ratio of UV to IR energy scales Λ UV /Λ IR can multiply loop corrections, enhancing them relative to their naive 3 A partial way out of this argument is if a model is close to a phase transition already at zero temperature, then only small thermal corrections are needed to undergo the transition. However, this setup requires a hierarchy of scales already at tree-level, whereby the potential difference between the minima is parametrically small compared with the curvature of the potential. loop counting. EFT provides a systematic means to account for such enhancements. To construct a reliable perturbative expansion to describe a given energy scale Λ, one must first integrate out all energy scales which are parametrically larger. In studying thermal phase transitions, the first step therefore is to identify the energy scale of the transitioning field.
In constructing the EFT for the transitioning field, one integrates out heavy degrees of freedom step by step. Each heavy degree of freedom that is integrated out modifies the effective infrared dynamics, including the effective mass of the transitioning field. How then do we know when to stop integrating modes out? After constructing the EFT for energy scales Λ and below, if the mass of the transitioning field remains of order Λ through the transition, then one can be sure that all contributions which are enhanced by a ratio of scales have been captured. On the other hand, if there is an apparent second-order phase transition, then the effective mass of the transitioning field goes to zero, and any other massive degrees of freedom will become relatively heavier than the transitioning field. Thus, new hierarchies of scale arise, and the fields which remain of mass Λ through the transition must be integrated out.
The same conclusions can be reached from a rather different perspective. In general, a successful perturbative expansion requires a LO approximation which is relatively close to the complete result. If the LO approximation is qualitatively different from the complete result, then perturbation theory will fail, and may exhibit all manner of weird and wonderful pathologies. In the study of phase transitions, such pathologies arise when the LO approximation fails to get the order of the phase transition right.
At nonzero temperature, infrared modes (with energy E ≪ T ) of bosonic fields become highly occupied, and their collective effective coupling is enhanced. For strong first-order phase transitions, the transitioning field remains gapped, and weak-coupling (i.e. mean-field) expansions can work rather well, as long as the effective coupling is small. For weaker transitions, the bosonic field undergoing the transition is relatively lighter, so the convergence of the expansion is slower, until for transitions of second order or higher, the effective coupling is large, and there is therefore no weak-coupling expansion. One must then resort to other approaches, such as bootstrapping [45], lattice Monte-Carlo [21,46], weak-strong dualities [47], or the ϵ expansion [48,49].
Combining the observations of the previous two paragraphs, we reach the following conclusion: when studying phase transitions using a weak-coupling expansion, one should always start with a LO approximation in which the transition is of first order. This is the best that one can do with perturbation theory. If the transition is indeed of first order, then the perturbative expansion will converge well. On the other hand, if the transition is really of second order or higher, then perturbation theory will fail, but there is anyway no way around it.
In this section, we present a generic setup for the perturbative analysis of thermal phase transitions. We will use EFT to construct LO approximations in which a given phase transition is of first order. This will provide us with the best possible starting position for perturbation theory and gives results which are gauge invariant, real and free from spurious infrared divergences. It also improves agreement with the lattice, even when the transition is not of first order.

Thermal scale hierarchies
The starting point of our computation is the EFT picture for thermal phase transitions [17,18,50,51]. This starts from the assumption that we are at high temperatures compared to relevant mass scales, T ≫ m, and is based on the following chain of scale hierarchies in terms of a weak coupling g ≪ 1, and the temperature T . The factors of π arise from Matsubara modes (πT ) and loop integrals ( g 4π ). However, from here on we shall omit the factors of 4 related to loop integrals, as they are often compensated by group theory factors in Feynman diagrams multiplying the loop integral. The hard, semisoft, soft and supersoft scales are perturbative; for these scales the effective expansion parameters are small, ε eff ≪ 1. Indeed, there are separate expansion parameters for each energy scale, with ε hard ∼ (g/π) 2 , and the expansion parameters for softer scales are larger, indicating slower convergence. We dedicate the next section to discuss how EFT expansions arise for the semisoft, soft and supersoft scales. Energy scales higher than the hard scale are exponentially (Boltzmann) suppressed. At the other extreme, energy scales at or below the ultrasoft scale are nonperturbative, as the effective expansion parameter therein is of order unity [20]. However, as we argue below, the dynamics of strong first-order thermal phase transitions generally takes place at either the soft or supersoft scales. Only for very weak transitions can the dynamics take place at the ultrasoft scale. In principle, other energy scales between the hard and ultrasoft scales may arise, though we have not encountered them.
In applying the above power counting to a given model, the parameter g should be chosen such that ε hard ∼ g 2 /(π) 2 determines the convergence of the loop expansion for the hard scale. Thus, for models with a single dimensionless coupling, g can be identified with a 3-point coupling, or g 2 with a 4-point coupling. For models with multiple couplings g should be identified with the largest relevant coupling, as this is what limits the convergence of the loop expansion.
Hard scale The temperature sets the most UV scale for thermal fluctuations, as energies above this are Boltzmann suppressed. At high temperatures, these hard scale fluctuations dominate the free energy density. For equilibrium physics, Matsubara's imaginary time formalism reveals that nπT sets the energy scale of thermal-scale fluctuations, where n is an even integer for bosonic fields, or an odd integer for fermionic fields. All modes except the (bosonic) zero Matsubara mode n = 0 therefore have energies of at least the hard scale. Fields with masses above the hard scale can be integrated out as at zero temperature [52][53][54]. For the hard scale fluctuations, each successive loop is suppressed by ε hard ∼ g 2 /(π) 2 compared to previous one.

Soft scale
The effective dynamics of softer modes is screened by hard scale fluctuations. At one-loop order this screening induces an effective mass of order gT . Bosonic zero Matsubara modes are therefore generically of the soft scale, unless there is some mechanism for the partial or full cancellation of one-loop screening. The EFT construction between the hard scale and the soft scale is the well-known high-temperature dimensional reduction [17,18]. Technically, this amounts to integrating out nonzero Matsubara modes with masses of order πT , and constructing the EFT for the three-dimensional zero Matsubara modes of all lighter bosonic fields. The temporal components of gauge fields acquire a thermal Debye mass due to the heat bath breaking Lorentz invariance [55,56]. These modes always live at the soft scale. Their squared Debye masses are solely generated by screening of the hard scale, so are a sum of positive definite terms each of order (gT ) 2 . Heavy bosonic zero-modes, with masses comparable to the hard scale πT , are integrated out along with the nonzero Matsubara modes [31,34,39,52,53]. Within the soft scale EFT, each successive loop is suppressed by ε soft ∼ g/π compared to previous one, at least when the interaction corresponding to g is present at the soft scale.
Supersoft scale The thermal effective masses of Lorentz scalar fields can be parametrically smaller than the soft scale. If the quadratic mass parameter of a scalar field is negative at zero temperature µ 2 < 0, then there can be cancellations between this and the positive hard thermal contributions to the effective mass, where the subscript 3 is used to denote the effective mass of the field in the 3d EFT and c > 0 is an O(1) numerical coefficient. 4 For broad classes of transitions, the transitioning field becomes lighter than the soft scale at the critical temperature. The dominant subleading corrections to its effective squared mass come from integrating out soft-scale fields, and are of order O(g 3 T 2 /π). Thus, barring additional parametric cancellations, the mass of the transitioning field lives at the supersoft scale g 3/2 T / √ π. This is the case for symmetrybreaking first-order phase transitions [15,25]. The supersoft scale is in fact even more widely applicable to thermal first-order phase transitions, as we find below.
The construction of an EFT for the supersoft scale was introduced in Ref. [26]. The transition of a supersoft field can modify the masses of soft-scale fields at leading order. In this case, the effective Lagrangian at the supersoft-scale will have non-polynomial dependence on the transitioning field. It is nevertheless local, as the hierarchy of scales ensures the derivative expansion holds. This is akin to the EFT of inflation [58]. Within the supersoft scale EFTs we consider, the supersoft scale fields have parametrically small couplings at zero temperature, and each successive loop is suppressed by ε super ∼ (g/π) 3/2 compared to previous one.
Semisoft scale This scale lies between the hard and soft scales. It can arise naturally for very strong first-order phase transitions: if the jump in a scalar background field becomes as large as √ πT / √ g, then it can impart a mass of order √ πgT on other fields through the Higgs mechanism. Below we find that this situation occurs in Z 2 -symmetric multi-field models where there are two successive first-order phase transitions, and where there is sufficient supercooling between them. For such setups, when integrating out the semisoft scale fluctuations, each new perturbative order is suppressed by ε semi ∼ g/π compared to previous one, though multiple orders in this expansion appear at each loop order (suitably resummed).
Ultrasoft scale This scale and below are nonperturbative. In gauge theories, the spacelike Ward identities ensure that the spatial components of gauge bosons do not receive a thermal mass correction within perturbation theory. In the absence of a Higgs mechanism, they therefore remain massless until the ultrasoft scale, where they receive a nonperturbative thermal mass [20]. Lorentz scalar fields can also become ultrasoft in the near vicinity of a second-order phase transition. However, unlike for gauge fields, fine tuning is typically required for a scalar to be this light. The contribution of the ultrasoft scale to the free energy density is of order T (g 2 T /π) 3 , and hence subdominant to the contributions of the higher energy scales.

Effective field theory expansions
We begin by discussing the thermal effective potential, or the background-field-dependent free energy density of the plasma. From this the phase transitions of a model can be determined. We discuss the computation of the effective potential for a transition taking place at the soft or supersoft scale. The scale inducing such a transition must be heavier than the transitioning field, hence it can be the hard, soft or semisoft scale. A starting point of our discussion is the tree-level potential for the soft scale 3d EFT, in terms of real background fields v 5 By tree-level we mean that no loop diagrams from within the soft scale EFT are included. However loop diagrams from the hard scale are included in V soft tree , and are captured in the parameters of the EFT. This is reflected in the expansion in ε hard on the right hand side, where we will assume that the first two orders have been calculated. The calculation of the NLO term was pioneered in Refs. [17][18][19] and is a mainstay of high-temperature dimensional reduction. It has now been automated for generic models [59].
Based on the scale hierarchy between inducing and transitioning scales, we present EFT expansions of the effective potential. In such expansions, perturbation theory is organised in 5 Here, we indicate the size of the potential in terms of the mass term, quadratic in v. In the vicinity of a phase transition, terms with other powers of v are of comparable size. Note that the mass dimension of v is 1/2, following from the canonical normalisation of a scalar field in 3d. The relation to the corresponding 4d field is v 2 ∼ v 2 4d /T . terms of power counting with respect to dimensionless quantities within the 3d description. In Sec. 4, we then discuss the computation of thermodynamics using these EFT expansions. The validity of our perturbative EFT expansions will depend on the magnitude of the background field v 2 . For very weak transitions, a field with effective coupling g 2 3 ∼ g 2 T and mass ∼ g 3 |v| will become nonperturbative unless v 2 ≫ g 2 T /π 2 . On the other hand, for very strong transitions, the thermal scale hierarchy will break down altogether unless v 2 ≪ π 2 T /g 2 . These conditions ensure that a particle of mass g 3 |v| is much heavier than the ultrasoft scale, and much lighter than the thermal scale πT . Together they determine the relatively wide range of transition strengths which we can describe perturbatively, In the context of electroweak baryogenesis, the geometric midpoint v 2 ∼ T is favoured in simple BSM models which produce the observed baryon asymmetry [1,60]. The midpoint also often makes a convenient choice for power counting, but in what follows we will discuss both weaker and stronger transitions.

Transition for a soft field
As argued above, we wish to construct a perturbative expansion which predicts a first-order phase transition at leading order, which in this case means at tree-level within the soft-scale EFT. This implies a certain structure for the tree-level soft-scale potential: there should be at least two coexisting local minima, separated by a potential barrier. In the absence of such a tree-level barrier, we argue that there are no soft-scale phase transitions that can be described reliably within perturbation theory. Assuming the theory is weakly coupled at zero temperature, we can then conclude that either there are no transitions at all, or the transition takes place further into the IR. The perturbative expansion of the soft scale effective potential can be expressed in terms of a formal expansion in ε soft which is a dimensionless ratio of EFT parameters. Such an expansion parameter inherits its scaling ε soft ∼ g/π from the original theory, but can be treated as an independent expansion parameter in the following sense. Within the EFT, perturbation theory can be organised as an expansion in ε soft , and such a computation can be used to find critical values for 3d parameters, and the condensates as functions of the 3d parameters. Then, one wants to relate these to the temperature and the original parameters of the parent 4d theory, and this is done in an expansion in ε hard ∼ g 2 /π 2 in dimensional reduction. Indeed, there are two different expansions, one related to UV physics at the hard scale, and another to IR physics at the soft scale. The effective potential at the soft-scale admits the formal expansion where we have introduced V soft 0 ≡ V soft LO to simplify notation, and denote higher order corrections with increasing numeral in the subscript. Here the scaling of the LO potential is indicated in terms of the original weak expansion parameter of the parent theory, and we have assumed v ∼ T for simplicity. The magnitudes of higher order corrections are given with respect to LO in terms of ε soft ∼ g/π.
Diagrammatically, the expansion in ε soft aligns with the loop expansion within the EFT. The computation of the effective potential up to two loops is straightforward, and has recently been automated for general models [59]. It is illustrated in Fig. 2. Here dashed lines represent all scalars, wiggly lines gauge fields and dotted lines ghosts. The next term, of order O(ε 3 soft V soft 0 ), arises at three-loops. This is the last order which is computable in perturbation theory in theories with non-Abelian gauge fields [20]. For later convenience, we denote next-to-next-to leading order as N2LO and higher orders with increasing numeral.
Possible realisations of a soft-scale EFT showing a first-order phase transition at tree-level include a real scalar with cubic and quartic interactions [61], a real scalar with quartic and sextic interactions [35], and multi-scalar models where there is a transition between different broken phases [36]. The tree-level potential of the cubic-quartic model, here written in terms of a real background v, reads where the linear term has been removed by a shift v → v + const. 6 For κ 2 3 > 4λ 3 m 2 3 , this potential admits two minima separated by a maximum, and these two minima have the same height when κ 2 3 = (9/2)λ 3 m 2 3 . Starting with the potential of Eq. (3.4) as our LO approximation, let us consider loop corrections. These arise both from the soft scale (within the EFT) and from the hard scale (the construction of the EFT). The loop expansion parameters within the soft scale EFT are The powers of κ 3 and λ 3 follow from standard graph-theoretic identities [62], the inverse powers of mass arise from loop integrals and can be determined from dimensional analysis, and the factors of (4π) follow from the angular integrals arising in 3d loop integrals. These should be compared with the loop expansion parameters arising within the corresponding 4d theory (with analogous parameters dropping subscripts 3), which determine corrections arising from the hard scale.
In the vicinity of the critical temperature, the joint requirements that there are two coexisting minima, and that their potential energies are approximately equal, imply that all three terms in the potential are of the same order, so that κ 2 3 ∼ λ 3 m 2 3 , and hence that the convergence of the loop expansion within the EFT is determined by a single expansion parameter Using λ 3 ∼ g 2 T and m 3 ∼ gT together implies ε soft ∼ g/(4π). The loop expansion within the EFT therefore converges more slowly than the loop expansion used in constructing the EFT. In addition, the expansion within the EFT diverges for a second order phase transition m 3 → 0, though in the approach to this point the field becomes lighter than the soft scale.
For the Z 2 -symmetric model with quartic and sextic interactions, the potential reads where λ 3 < 0. The general conclusions about the perturbative expansion in the cubic-quartic model carry over to this case. In the vicinity of the critical temperature the loop expansion parameter within this EFT is [28,35,63]. For λ 3 ∼ g 2 T , the soft expansion parameter is again of order ε soft ∼ g/(4π).
Our third example is provided by a phase transition with two scalar fields participating in the transition, such that the transition happens between the different broken minima of the potential. For simplicity, we assume here a Z 2 -symmetric model where scalars are charged under a gauge group with a gauge coupling g. Concretely we discuss the following tree-level potential with two background fields x and y V soft tree,x,y = For power counting, we assume that the scalar masses lie at the soft scale, µ 2 x,3 , µ 2 y,3 ∼ (gT ) 2 , as well as possible gauge field and Debye masses m W , m D ∼ gT , and that the couplings are all equally perturbative λ x,3 , λ y,3 , λ xy,3 ∼ g 2 T . The loop expansion parameter within the EFT is then ε soft ∼ λ x,3 /(4π|µ x,3 |) ∼ g/(4π).
Depending on the signs of the mass terms, there is a symmetric minimum where both background fields vanish (x, y) = (0, 0), and broken minima at (x 0 , 0) and (0, y 0 ), where For certain choices of parameters, the symmetry breaking pattern (0, 0) → (x 0 , 0) → (0, y 0 ) is possible. The second step can be described reliably by this soft-scale EFT as, in field space, the two broken minima are separated by a barrier, so the transition is of first-order within the EFT.

Transition for a supersoft field
We continue by discussing the generic setup for a supersoft scale EFT of a single scalar field with a tree-level potential barrier. By tree-level we mean that the potential does not include any loop corrections from the supersoft scale, yet note that it can still be non-polynomial due to contributions from the soft scale fields that have been integrated out. Indeed, the barrier is typically generated by the soft scale. Our following discussion is inspired by [15,[25][26][27][64][65][66][67] that discuss one-step phase transitions with radiative barriers. For a two-step EWPT, this setup describes the first transition to an intermediate phase before the transition to the EW phase.
Our discussion here is schematic: we assume a single supersoft scalar field ϕ, with real background field v, that couples to a gauge field, with 3d gauge coupling g 3 . The one-loop effective potential at the soft scale reads [68] The first two terms are tree-level terms in the soft scale EFT, and the last term is a one-loop contribution where m 2 V ≃ g 2 3 v 2 /4 is the gauge boson mass eigenvalue and M represents the scalar mass eigenvalues Here m D is the Debye mass for the temporal component of a gauge field, and h 3 its coupling to ϕ.
For there to be a first-order transition, the potential should have more than one minimum separated by a barrier. At tree-level in the soft scale EFT this is not possible, since there is only one minimum at any given temperature: when the mass parameter µ 2 3 is positive, the only minimum is at v sym = 0, and when it is negative the minimum is instead at v broken = −µ 2 3 /λ 3 . At µ 2 3 = 0 there is a second-order transition. However, the one-loop gauge boson contribution is cubic ∝ |v| 3 and can provide a barrier between minima. Similarly, contributions from the temporal components of gauge fields contribute to the barrier -as well as other soft scalar fields in the case of multiple scalars -but in order to simplify the presentation we do not discuss their contribution further in this section, but assume that the 3d gauge field is the only soft field. For a transition to be of first order, the one-loop gauge boson term should be of the same parametric order as the tree-level terms of the potential, i.e.
Additional loop corrections from gauge bosons are suppressed relatively by g 2 3 /(πg 3 |v|), and those from the scalar undergoing the transition are suppressed by λ 3 /(πM 3 ). To keep track of perturbative corrections, we introduce the dimensionless power counting parameter ε soft ∼ g 3 /(π|v|), which counts soft-scale loops and satisfies 0 < ε soft ≪ 1 when perturbative corrections are small. Expressing all three parameters in units of v, Eq. (3.11) then implies that Here, and in what follows, we use ∼ to denote that two quantities have the asymptotic scaling as ε soft → 0 + . Note that this power counting is equivalent to that of Refs. [15,30] when ε soft ∼ g/(4π), and often the ratio ε soft ≡ λ 3 /g 2 3 is used as an expansion parameter and denoted by x in previous literature, e.g. [23]. Note in particular that with this counting, the one-loop scalar terms are of order ≃ M 3 /π ∼ µ 3 3 /π ∼ π 2 ε 9/2 soft v 3 and hence do not contribute at leading order to the potential Here we have indicated that this is the LO, or tree-level, potential of the supersoft scale EFT. The mass of the transitioning field is i.e. at the supersoft scale. Here Π ≡ − 3 8π g 3 3 |v| is the resummed contribution from the soft gauge field with mass m 2 The resummation arises from integrating out the soft fields, and it decorates the supersoft scalar propagator with one-loop insertions of the soft field. The broken minimum of the LO potential reads (

3.15)
This minimum is separated from the symmetric minimum at v sym = 0 by a barrier for temperatures such that the effective mass term lies in the range 0 < µ 2 3 < 9g 6 3 /(1024π 2 λ 3 ). We comment that should we include contributions of other soft fields, the expression for the broken minimum becomes readily much more complicated analytically. In addition, we point out that the supersoft EFT is constructed in the broken phase, or for sufficiently large background fields, where the gauge field is indeed soft and can be integrated out. Next, we consider higher orders in ε soft . The NLO corrections to the effective potential, suppressed by one power of ε soft , are given by two-loop digrams of purely soft modes. This can be formally performed by matching the effective potentials of the soft and supersoft EFTs, treating supersoft masses and momenta in strict perturbation theory, cf. [26,27,54,58]. In this approach, the propagator of the supersoft field is treated as massless for the matching computation [18], so that pure supersoft diagrams vanish identically, and only soft-scale contributions from mixed scalar/gauge diagrams are included. The outcome of this formal procedure is the potential for the supersoft EFT, and the procedure is illustrated in Fig. 3. In Appendix A.1, we present explicit computations up to N2LO for some relevant example models.
In the construction of the effective potential for the supersoft scale, there are multiple expansions. First, there is the expansion related to integrating out the hard scale, which we will here take for granted. Second, there is the expansion related to integrating out the soft scale, and finally there is the loop expansion within the supersoft scale EFT. The parameters for these latter two expansions are related as ε super ∼ ε 3/2 soft , so that ε super ∼ (g/π) 3/2 . Oneloop diagrams within the supersoft theory contribute at O(ε super ) relative to LO and have the simple form in terms of the resummed mass M depicted with the double line in Fig. 3. Note that at the order we work, and because we are only interested in observables for homogeneous background fields, we do not need to include the effect of the momentum dependent field normalisation in matching, c.f. e.g. [27]. The N3LO contributions to the potential, or O(ε 2 soft ) relative to LO, are given by threeloop soft-scale diagrams. This is followed by terms of order ε soft ε supersoft ∼ ε 5/2 soft relative to LO, which are given by the resummation of the NLO corrections to the mass within the supersoft one-loop diagram [54]. We do not compute either of these contributions, leaving them for future work. N4LO is the highest order that can be computed perturbatively, since at the next order one encounters Linde's Infrared Problem [20], where all loop topologies of the ultrasoft scale contribute at the same order in powers of couplings. These contributions are suppressed relative to the LO potential by ε 3 soft ∼ (g/π) 3 , in terms of the weak coupling of the original theory.

A two step phase transition
Let us return to the example model with two background fields (x, y), and the potential of Eq. (3.8). We consider the interesting case where there is a two step transition with the first step (0, 0) → (x 0 , 0) taking place at the critical temperature T c,1 followed by a second step (x 0 , 0) → (0, y 0 ) at a lower temperature T c,2 . Here x 0 and y 0 are generic nonzero background expectation values for x and y.
The thermodynamics of such two-step phase transitions depends on the relative magnitudes of the couplings and masses. As argued above, for perturbation theory to work, we need to find EFTs for the transitioning fields in which these transitions appear first order. We find there are (at least) two natural options for the power counting relations, which we outline below.

The first step
The first transition appears to be of second order within the soft-scale EFT. So, the transitioning field x becomes parametrically lighter than the soft scale.
Let us start by assuming that all couplings are equally perturbative λ x,3 ∼ λ y,3 ∼ λ xy,3 ∼ g 2 T , following the discussion after Eq. (3.8). Then, integrating out the soft scale fields around this transition, one finds that the largest possible discontinuity in the background field is x 2 0 ∼ g 2 T /(4π). At this point the effective mass of the transitioning field is at the non-perturbative ultrasoft scale, and the transition is either very weak or a crossover and perturbation theory is not viable. This conclusion crucially relies on the absence of hierarchies between the couplings. The first step (0, 0) → (x 0 , 0) can be strongly first order if λ x,3 is parametrically smaller than some other couplings to the x field. For example, if the power counting for the portal coupling is unchanged λ xy,3 ∼ g 2 T but the self coupling scales as λ x,3 ∼ g 3 T /π, then the x field is supersoft at the first transition, and the analysis of Sec. 3.2 applies directly. For the first transition, the leading order potential then reads where m 2 y ≃ µ 2 y,3 + 1 2 λ xy,3 x 2 , and we have included contributions from a vector boson with mass m 2 V = g 2 3 x 2 /4. We assume that only the x field becomes supersoft as it transitions, with other fields remaining soft. The balance of the three terms in the potential then implies that x 2 ∼ T . In the supersoft scale EFT, the mass of the transitioning field is resummed We emphasize that the effective couplings of the soft scale EFT do not vary much with respect to temperature, apart from their dimensional scaling, i.e. ratios λ x,3 /T , λ y,3 /T and g 2 3 /T are approximately constant, whereas scalar mass parameters, µ 2 x,3 and µ 2 y,3 , can vary non-trivially with temperature, in particular they can go through zero. Fig. 4 shows a schematic plot of the temperature dependence of the effective masses in the vicinity of the first transition. The left panel shows a generic scalar mass parameter, which grows with temperature. Shown in the right panel, however, when contributions from the background field x are included, the scalar masses become decreasing functions of temperature below T < T R 0 for which the x-minimum exists. The right panel depicts a generic soft mass M and the resummed x-field mass in the x-phase of the supersoft EFT. Additionally in the left panel the Debye mass is shown, together with contributions from the background field. The Debye mass is typically significantly heavier than scalar masses. If the transition is weak, so that the correction to the Debye mass due to the background field is relatively small, it should be possible to integrate out the corresponding field following Ref. [17].

Two consecutive supersoft scale transitions
Since 3d effective couplings do not vary significantly with temperature, we can assume they satisfy the same formal power counting relations for the second step (x 0 , 0) → (0, y 0 ), as for the first. Furthermore, we also assume that λ y,3 ∼ λ x,3 . 8 Given this, what is the scale for the second transition? 8 In practice -for strong electroweak phase transitions in models with relatively heavy BSM fields, and with a generic gauge coupling g4-the portal coupling λxy is often the largest coupling, and there is a hierarchy λy ∼ λx < g 2 4 < λxy [37,[69][70][71]. This suggests to organise all power countings with respect to λxy instead of the gauge coupling. for both x-and y-fields), whereas the Debye mass is noticeably larger (note that the Debye mass is enhanced by group theoretic factors); (ii) for the final mass eigenvalues, there is a clear hierarchy between the light transitioning field and other fields. Note that such hierarchies might not be so clearly manifest at any given parameter point, yet this plot is inspired by our numerical application in Sec. 5.
The first possibility is that there is no significant supercooling between the critical temperatures, and while the masses of both x-and y-fields increase with increasing background field, they are still supersoft at T c,2 , in the x-and y-phases respectively (note that the x-field is soft in the y-phase, and vice versa). Alternatively, given enough supercooling down from T c,1 , the masses of transitioning fields could grow to become soft ∼ (gT ) 2 .
First, let us assume that the second transition to the y-phase occurs without much supercooling. Then the x-phase free-energy should still be computed in the supersoft scale EFT which described the first step, Eq. (3.17). What about the y-phase free energy?
Given the assumptions that λ y,3 ∼ λ x,3 ≪ λ xy,3 ∼ g 2 3 , the mass hierarchies in the y-phase are mirror those in the x-phase. The y-field is supersoft in the y-phase, and all else is soft. At leading order, in the vicinity of T c,2 the y-phase potential therefore reads where m 2 V = g 2 3 y 2 /4 ∼ (gT ) 2 and m 2 x ≃ µ 2 x,3 + 1 2 λ xy,3 y 2 ∼ (gT ) 2 , i.e. the gauge field and x-field are soft in the y-phase. The critical temperature T c,2 is determined from the condition that

Transition for a soft field, induced by semisoft scale
On the other hand, if there is enough supercooling between T c,1 and T c,2 , the masses of the transitioning fields can grow to the soft scale at the second transition. This is illustrated schematically in Fig. 5. In this case, background fields at minima are parametrically larger, and the leading order potential agrees with the standard tree-level potential. In the x phase, this is 20) which implies that at the broken minimum x 2 ∼ µ 2 x,3 λ x,3 ∼ π g T . Similarly, for the y-phase the background field satisfies y 2 ∼ µ 2 y,3 λ y,3 ∼ π g T . This is schematically illustrated in Fig. 6. In the x-phase, the large background field of x enhances and dominates the mass of y, and vice versa for the mass of x in the y-phase, This new scale hierarchy allows us to construct a soft scale EFT, where we integrate out the semisoft scale. The minima in the broken phase of these EFTs scale as x 2 T ∼ π g ≫ 1 (and likewise for the y-field) and describe an extremely strong transition. The critical temperature T c,2 is determined from the condition that the free-energies of both phases are equal, where each is computed within a separate soft-scale EFT, with fields at the semisoft scale integrated out. Note, that in the field space of two fields, there is a barrier at leading order that separates the phases.
At leading order the size of the field-dependent effective potential can be read off from Eq. (3.1) with v 2 ∼ πT /g, leading to V soft eff ∼ gπT 3 . Counting powers of couplings, and ratios of scales, the perturbative expansion in each phase takes the form where ε semi ∼ g/π and we have used the shorthand notation V 0 ≡ V soft LO and truncated the series to fifth order, i.e. N5LO. A diagrammatic rundown of these corrections is depicted in Fig. 7. Computation beyond that requires soft mass insertions at two-loop order, which would provide the result at N6LO, and 3-loop semi-soft scale diagrams [72] are required for N7LO. The computation of the effective potential splits into two: first, matching from the semisoft to soft scale, and then soft scale contributions. In the matching, all loop momenta are formally semisoft [54], and one expands propagators in soft mass parameters, before integration. For example, in the y-phase the one-loop x-field bubble diagram is expanded as 9 where M 2 (y) ≡ 1 2 λ xy,3 y 2 is the semisoft contribution to the mass. Utilising such expansion in the soft mass, in To compute the full N6LO piece would require performing similar expansions at two-loop level. For example, the pure scalar sunset integral is , (3.25) where m 2 y is soft and m 2 x = µ 2 x,3 +M 2 (y) includes both soft and semisoft contributions. Again, for the matching one first expands the integrand in the soft-scale quantities m 2 y and µ 2 x,3 , and then evaluates the integrals with semisoft momenta in the loops. For this, we can make use of the generic two-loop result given in Eq. (C.81) of Ref. [73] (c.f. also [74] and references therein) and define .

(3.26)
Expanding the integrand of Eq. (3.25) in soft masses m 2 y and µ 2 x,3 , and expressing the result in terms of S αβδ , we obtain Here the last term shown describes soft mass insertions, that contribute at N6LO. Other twoloop diagrams could be treated in an analogous manner, first expanding the integrand with respect to soft masses and only then computing the resulting integrals. However, a similar treatment with integrals involving gauge field propagators is somewhat more laborious, and we have decided to truncate our computation to N5LO in this work at hand, and leave higher orders for future. We give a more detailed description of the computation depicted in Fig. 7 in Appendix A.2, working through a concrete example for the SM augmented with a real triplet. Finally, in the computation of the one-loop correction to the effective potential in the soft scale EFT, the mass of the soft field needs to be resummed as Resummation by the V 1 contribution (but not V 2 or higher) is needed at N5LO. The effect of matching the momentum-dependent field normalisation will contribute first at N7LO for observables depending only on homogeneous backgrounds fields.

Strict perturbative expansions
In the previous section, we have discussed the construction of EFTs to describe first-order phase transitions, and the calculation of their effective potentials. With this in hand, there are a number of different calculational approaches which one could adopt to analyse the phase structure and thermodynamics.
The equilibrium thermodynamics of a model can be derived from the pressure -determined by the effective potential evaluated at its minima -as a function of temperature, and its derivatives with respect to temperature. The pressure can be written as [75] where p 0 is the coefficient of the unit operator in the construction of the 3d EFT [18] and V eff is the effective potential within the 3d EFT [68]. This coefficient of the unit operator is linked to the symmetric phase pressure as p sym (T ) = p 0 − T V eff (0), where the effective potential is evaluated at the origin [75,76]. The critical temperature T c is defined as the temperature where the pressure difference between two phases vanishes ∆p(T c ) = 0 and this translates to a condition that the effective potentials at different minima are degenerate. Gauge invariant condensates were already discussed in Sec. 1, and they can be computed as derivatives of V eff with respect to the parameters of the 3d EFT. The strength of the phase transition can be characterised in terms of released latent heat, which is related to the pressure as L = T ∆p ′ , where prime denotes a derivative with respect to temperature. All these quantities depend on differences between phases, and hence we do not need to compute p 0 . However, note that the speed of sound in each phase depends directly on p 0 , and this is relevant for the determination of the gravitational wave power spectrum, see [75,77]. Next, we turn to different calculational approaches. The most direct approach would be to just numerically minimise the effective potential. There are however a number of issues with this approach. First, the effective potential is generically complex, with imaginary parts arising away from the minima of the leading-order potential where squared masses can become negative [78]. Typically this issue is simply ignored by working only with the real part, or by replacing squared masses by their absolute magnitude. 10 Second, the effective potential is gauge dependent, and so are its minima, when computed directly.
An alternative approach is to adhere strictly to the confines of the perturbative expansion, and to perform a strict expansion in powers of ε eff . This approach is sometimes called the ℏ expansion, though in general the expansion parameter need not have anything to do with Planck's constant. Rather than directly minimising the full effective potential, one first minimises the leading-order effective potential, and then includes the corrections from higher orders perturbatively. This approach has the benefit of being exactly gauge invariant orderby-order [40,41]. It is also manifestly real. The difference between this approach and direct minimisation is due to a subset of higher order terms in the expansion of the minima (the tadpole expansion), which are resummed by direct minimisation.

Method Gauge invariant
Real direct ✗ ✗ mixed ✓ for sufficiently small ε eff strict ✓ ✓ Table 1: Basic theoretical properties of different perturbative methods described in the text. The column headings here refer to the properties of real physical quantities computed using these methods, such as the free energy or critical temperature. Note that, while one can always get a real result from a complex quantity simply by discarding the imaginary part by hand, we do not consider such a result to be a genuinely real prediction of a given method. As Goldstone squared masses go through zero at the LO broken phase, the direct method can yield spurious imaginary parts even when perturbative corrections are arbitrarily small. The mixed method yields real physical results when the expansion parameter ε eff is sufficiently small, but as we argue around Eq. (4.19) below, it can yield spurious imaginary parts when higher-order corrections exceed some finite bound.
Further possibilities arise for quantities which require additional intermediate steps in their computation from the effective potential, such as the critical temperature. In these cases, one can choose to make an additional strict expansion, or to mix the strict expansion of the potential with a direct approach at solving ∆p(T ) = 0 for the critical temperature. Unlike for the minima of the effective potential, in this case both possibilities are real and gauge invariant. For the critical temperature, there is however an important difference between these two approaches: If the critical temperature at some higher order is not within the range in which there is metastability at leading order, then the direct approach fails, while the strict perturbative approach continues to work. This issue is discussed further below. Table 1 summarises the theoretical properties of the different perturbative expansion schemes we have considered. In Section 5, we will further describe and test all these different approaches for a numerical example, comparing them to lattice Monte-Carlo data. For the remainder of this section, we formulate strict expansions for the effective potentials of the previous sections, as well as expansions for thermodynamic quantities of interest. In addition, we describe mixed approaches that combine direct and strict methods.

Strict expansions for a soft field
In the case of a soft-scale field undergoing a phase transition, the effective potential has an expansion in the effective expansion parameter ε soft , The minima of the potential can also be expanded as where the coefficients v i are determined by solving for the minima of the potential as an expansion in ε soft . For a single field, this results in [40] V where we have introduced ∂ v ≡ ∂ ∂v , and all terms on the right-hand side of the equalities are evaluated at v 0 . This last point is crucial for the desirable properties of this expansion, such as order-by-order reality and gauge invariance. In analogy, we can write the pressure as where the coefficients incorporate corrections from the expansion of minima p 0 = −T V 0 , . The generalisation to multiple scalar fields follows by upgrading the multiplications in Eqs. (4.4) and (4.5) to matrix multiplications. The background field becomes a vector v a with index a, and ∂ a ≡ ∂ ∂ v a is the corresponding gradient operator. The second derivative of the potential then becomes a matrix, and we find where (∂ a ∂ b V 0 ) −1 is the matrix inverse of the Hessian matrix ∂ a ∂ b V 0 , and we have used Einstein summation convention. As above, all terms on the right-hand sides of these equations are evaluated at the LO minimum v a 0 . For a 2-field model, with v = (x, y), one can invert the Hessian matrix explicitly, resulting in [37] where again the right-hand side is evaluated at the LO minima v a 0 = (x 0 , y 0 ). The LO minima v 0 describe different phases of the system, and the most likely phase in thermal equilibrium corresponds to the global minimum. Later on, we refer to Eqs. (4.4), (4.7) and (4.9) as the soft EFT expansion for the effective potential. Technically, this is a strict expansion of the effective potential around its LO minima, in a 3d EFT at the soft scale. Notably, this expansion is gauge invariant order by order [40,42], since it satisfies the Nielsen-Fukuda-Kugo identities within the EFT [40,41], and the construction of the EFT through dimensional reduction is gauge invariant [27,35].
The strict expansion strategy can be extended to determining the critical temperature [42], by writing where T 0 is solved from ∆p 0 (T 0 ) = 0, (4.11) or equivalently −∆V 0 (T 0 ) = 0. The next two orders in the expansion are 11 Similarly, expanding the latent heat where We emphasize that in consistent EFT expansions we do not encounter the problem reported in [42], whereby the condition ∆p 0 (T 0 ) = 0 leads to a vanishing mass parameter µ 2 3 (T 0 ) = 0 which then results in spurious IR divergences at two-loop order for T 2 . This problem has its roots in the fact that the leading order potential used in [42] is the one at tree-level at the soft scale, and does not have a barrier, but describes instead a second order phase transition. In the approach we have advocated, the assumption of a barrier in the LO potential is build-in to the EFT construction. In the case of an apparently second-order transition at the soft scale, the correct description for the transition is the supersoft scale EFT, cf. Secs. 3.2 and the discussion below in this section. Furthermore, EFT expansions are manifestly real: since all expressions are evaluated at LO minima, all squared mass eigenvalues are non-negative and hence no imaginary parts arise from the computation of higher order corrections.
The combination of the soft EFT expansion for the effective potential and strict expansions for T c and L is the strict method of Table 1. This method is computationally 11 We remind the reader that prime is used to denote the temperature derivative, V ′ 0 ≡ dV 0 dT . very efficient on an algorithm level, since once T 0 is solved numerically from the condition ∆p 0 (T 0 ) = 0, all corrections are simply evaluated at T 0 from expressions that are known analytically. 12 Expansions in ε soft for the condensates follow naturally from the EFT expansion for the effective potential. Concretely, for the scalar quadratic condensate and where for simplicity we assumed that the scalar is a Higgs doublet.
Note that since v 0 often depends on the 3d mass parameter µ 2 3 , one first evaluates the effective potential in each phase, and only then differentiates.
Alternatively, following a strategy of [37,79,80], one could determine the critical temperature by computing ∆p(T ) using the soft EFT expansion, and then numerically solving for its root ∆p(T c ) = 0. Such a direct determination of T c is indeed computationally efficient, as it does not require numerical minimisation of a complicated effective potential. This is the mixed method of Table 1: the effective potential is evaluated in a strict expansion around the leading order minima, while thermodynamic quantities are determined directly as functions of temperature. Indeed, the same method can be applied to the determination of the latent heat [37,80].
With the mixed method, or more generally if we consider evaluating quantities at temperatures other than T 0 , we come across the following problem: (here we denote the expansion parameter again by ε eff for generality) The LO minima v 0 exist in some range of temperatures T ∈ [T L 0 , T R 0 ]. The LO critical temperature lies in this range T 0 ∈ [T L 0 , T R 0 ]. In a power counting sense the width of the range is T R 0 − T L 0 = O(T 0 ), so for ε eff → 0 the full critical temperature T c = T 0 + ε eff T 1 + ε 2 eff T 2 + ... must also lie in this range. However, for a finite expansion parameter ε eff it is possible that T c = T 0 + ε eff T 1 + ε 2 eff T 2 + ... is outside the range. This leads to the problem that for any temperature dependent function F (T c ) an expansion F (T c ) = F 0 (T c ) + ε eff F 1 (T c ) + ... does not exist. Furthermore, even if the LO minima exist at T c , in the mixed method the range of existence of a given phase is fixed at LO, and does not change at higher orders.
The proposed solution is to consider physical quantities as functions of 19) to treat the difference as of leading order ∆T = O(T 0 ), and then to power expand everything in ε eff . The origin of the independent variable ∆T is fixed to the critical temperature orderby-order, so that the expansion cannot cause T − T c to change sign, or to grow too large in magnitude. This helps to extend the principles of the strict method to temperatures other than the critical temperature. Note that for a fixed T , the corresponding value of ∆T depends on the order to which we compute T c . It is the difference from the critical temperature in a given approximation. Now consider the expansion of some generic quantity F : Evaluating this at ∆T = 0 gives or, equivalently, So, the range of existence of phases is shifted at each order, by the amount that T c changes at that order. This essentially solves the problem of the static range. Note however that the width of the range does not change from order to order; a feature which still seems undesirable.
Strict expansions for a supersoft field In the case of a supersoft field, EFT expansions follow essentially the same logic. However, due to the different structure of the effective expansion, the resulting expressions are slightly simpler. The effective potential consists of two expansions Since ε super ∼ ε 3 2 soft , following Ref. [30] we can write formally where ε ∼ √ ε soft and V 1 = 0 identically. In this expansion, the soft and supersoft expansions are mixed together. In principle one could do everything in a fully EFT way, essentially resumming V supersoft eff = ε 0 super (V 0 + V 2 ) + ε super V 3 , i.e. both V 0 and V 2 are treated as LO within the supersoft EFT. This is in analogy to not mixing the hard and soft expansions in dimensional reduction. However, here we choose the former option and mix the expansions together, so the supersoft scale EFT should be understood in the sense of this mixed expansion. We note, that these two alternatives agree up to the order computed [34], yet resum different sets of formally higher order corrections.
We emphasize that the order ε 1 is not present in the effective potential, and this leads to multiple simplifications in formulae below. We formally expand the minima as (4.27) where v 1 = 0 since V 1 = 0. The expansion for the potential evaluated at the minimum reads This expression is particularly simple up to and including O(ε 3 ), since subleading corrections to the minimum start at v 2 = O(ε 2 v 0 ), and the condition ∂ v V 0 = 0 at v 0 ensures that this does not contribute to the potential until O(ε 4 V 0 (v 0 )). Later on, we refer to Eq. (4.28) as the supersoft EFT expansion for the effective potential. Technically, this is a strict expansion of the effective potential around its LO minima, in a 3d EFT at the supersoft scale. In analogy to case of the soft field in the previous section, this expansion is gauge invariant [27,30,32], and Eq. (4.28) can be utilised in the mixed or strict methods of Table 1.
For the strict method, we expand Again, we emphasize, that these EFT expansions are free from spurious IR divergencies reported in [37,42]. In the EFT expansion, a radiative barrier provided by the soft fields is included to the LO effective potential, and the condition for T 0 does not lead to a vanishing mass parameter. Hence, the spurious singularities at higher orders are avoided, and concretely T 2 and T 3 are finite. Similarly, the latent heat has an expansion By writing the pressure as Finally, the quadratic condensate has an expansion Strict expansions for a soft field, induced by semisoft scale At this point, after the previous discussions, the methodology of strict expansions should be clear. However, for the sake of completeness we repeat the corresponding discussion here. The only difference is the form of the expansion of the potential, which in the case of the semisoft-induced soft-scale EFT reads Compared to the previous sections, we have many more orders available, and this time it is the N2LO term that is missing (i.e. V 2 = 0) due to the nature of the matching between semisoft and soft scales. From this expansion all else follows by Taylor expansion. The minima, pressure, latent heat, field condensates and other thermodynamic quantities can all be expanded as First, one solves ∂ v V = 0 for the minima, with each successive order determined by a linear algebraic equation, av n + b = 0 where a and b are given in terms of lower orders. From this one can construct the pressure, the first few orders of which are where as always the expressions on the right hand sides are evaluated at v 0 and higher orders can be generated iteratively. Notably, unlike in the supersoft EFT case, where p 1 was zero due to V 1 being zero, in this case all orders for the pressure are nonzero despite V 2 being zero. Hence, solving ∆p(T c ) = 0, and computing L = T c ∆p ′ (T c ), the order-by-order results reproduce the same expressions at the first few orders that were already encountered in Eqs.  4) and (1.5). This completes the formal outline of our setup, and next we turn to applications.

Cosmological phase transitions
In this section, we will test the EFTs and perturbative expansions presented above, as applied to possible cosmological thermal histories. As introduced in Sec.1 we use the real-triplet extended SM as our concrete playground, in order to compare our perturbative EFT methods to the lattice results of Ref. [37]. We will also use the renormalisation scale dependence of our perturbative results to provide an intrinsic measure of their uncertainty. However, for the thermodynamics of this model, the lattice results are expected to be correct up to very small statistical uncertainties, so they are the ultimate arbiter.
We would like to comment on a subtlety in comparing to Ref. [37]. The lattice simulations of [37] were performed for a 3d EFT without the temporal components of gauge fields, their effects being captured by the parameters of the EFT. Such EFTs are commonly referred as "ultrasoft" scale EFTs, and have been studied using non-perturbative lattice simulations e.g. in [21,22,70,81,82]. However, the derivation of such theories does not require scalar masses to be ultrasoft. Indeed, the only necessary assumptions therein are that (i) m 2 3 ≪ m 2 D , i.e. scalar masses are much lighter than Debye masses, and (ii) where h 3 is a generic portal coupling between scalars and temporal gauge field components, and v is the Higgs background field. Based on the power counting arguments of the previous sections, the Higgs and triplet scalar fields are not expected to become ultrasoft, except in the near vicinity of a second-order phase transition. For a reasonably strong first-order phase transition, they are either of the soft or supersoft scales. As a consequence, the temporal components of the gauge fields should be treated as described in Sec. 2 in the construction of the EFT for the transitioning fields. In particular, assumption (ii) can break down when the background field v 2 becomes large. In this case, one should not integrate out temporal gauge field components as described in [17].
Nevertheless, in order to provide an apples-to-apples comparison with the lattice results of [37], in all the numerical computations of this work, we incorporate contributions from the temporal components of gauge fields as in [37,39]. Furthermore, in the lattice simulations of [37], dynamical effects of the U(1) subgroup in the 3d EFT were not included, hence we set g ′ 3 = 0 in perturbation theory as well, but note that we still keep effects of g ′ in the dimensional reduction matching relations, as in [37].
The effective potential for the soft-scale EFT of the real-triplet extended SM can be found in [37], up to N2LO, or two-loops. Explicit results for the corresponding supersoft-scale effective potential are given in Appendix A.1, and results for soft scale effective potentials including effects from the semisoft scale are collected in Appendix A.2.

One-step symmetry-breaking transition
The two benchmark points studied non-perturbatively in Ref. [37] both showed a succession of two phase transitions, with the pattern: symmetric to triplet to Higgs phase. 13 For now, 13 Here we use the phase phase transition somewhat loosely. Depending on the mass and couplings of the triplet scalar, the symmetry-breaking transition may be a smooth crossover, like the liquid-gas transition of  we will focus on the transition from the symmetric phase to the triplet phase.
As with other symmetry-breaking transitions, the symmetric to triplet transition appears to be of second order at tree level in the soft-scale EFT. Following the arguments of Sec. 2, this implies that in fact the transition takes place at lower energies, so the soft scale should be integrated out. Barring further cancellations, we thus expect the transition to take place at the supersoft scale. However, for completeness, and for comparison to the previous literature, we also consider the possibility that the transition takes place at the soft scale (though this will lead to IR divergences). For each EFT we consider each of the perturbative methods introduced in the previous section. The matrix of possibilities are summarised in Table 2.
We have computed the triplet condensates according to each of the matrix of possibilities shown in Table 2. Our results, together with the lattice results of Ref. [37] are shown in Fig. 8 for BM1, and in Fig. 9 for BM2. Exceptions are made for the strict method in the hard→soft EFT as well as the direct method in the soft→supersoft EFT, for which we do not show the results in Figs. 8 and 9. For the former approach, strict expansions fail as originally realised in [42]. At leading order in the soft theory the transition is of second order, and strict loop corrections do not modify the position of this transition, so that the mixed and strict methods in the hard→soft EFT are equivalent. While the use of direct minimisation in the soft→supersoft EFT has not appeared before in the literature, due to its gauge dependence there is no clear reason to prefer it to the other methods applied to the supersoft EFT, and we relegate the results to Appendix B.
In the following, we discuss in turn the panels of Figs. 8 and 9, focusing on the successes and failures of the different perturbative approaches. In all figures shaded bands correspond to the range of predictions from varying the 3d RG scale over the set Λ 3 ∈ {0.5T, T, 2T }, and therefore give an intrinsic measure of the theoretical uncertainty.
(a) Hard→soft EFT, direct method: At LO there is a 2nd order transition at T c ≈ 134. 5 GeV, which is the temperature where the triplet 3d mass vanishes, and is below the   lattice transition temperature. The NLO result shifts the critical temperature significantly above the lattice result, and the transition is of first order. Extending to N2LO yields much closer agreement with the lattice results for both T c and the values of the condensates, yet note that there are a number of data points missing, appearing as cuts in the otherwise continuous result. This is due to the possibility of our direct minimisation algorithm 14 failing. While this feature could be ameliorated by improving the 14 We used Mathematica's NMinimize function, adopting the differential evolution method and choosing tolerance parameters so that producing the N2LO data for one benchmark point took around one hour on a laptop, with temperature steps ∆T = 0.25 GeV.   Fig. 8, but for BM2. Note that the substantially larger portal coupling in this case leads to larger uncertainties at lower orders, and slower convergence. Nevertheless, the soft→supersoft strict method shows good agreement with the lattice at N2LO. minimisation algorithm, we have stuck to the aforementioned algorithm for the following reasons: the failure of direct minimisation at N2LO is a fairly common occurrence compared to using the same algorithm for the potential at lower orders. This is due to IR-sensitive logarithmic terms at two loops, that can result in spuriously large contributions to the potential at field values where the corresponding mass eigenvalues in the logarithm vanish, possibly preventing convergence to the actual minimum. Such an issue may be mitigated with a cost in performance time, thought this would be a limiting issue for model parameter space scans. To emphasize this aspect, we have used the same minimisation algorithm at all orders, despite the algorithm occasionally failing at N2LO, causing gaps in the corresponding plots in Figs. 8 and 9. Note that RG improvement kicks in for the first time at N2LO, due to the structure of running in these 3d EFTs [36,68]. The width of the error bands at LO and NLO are comparable, and in both cases significantly underestimates the theoretical error.
(b) Hard→soft EFT, mixed method: This method is known to fail when the LO potential does not have a barrier between the minima [37,42]. From Figs. 8b and 9b one can see two clear failures of this approach: T c is unchanged by higher orders, and at N2LO there is a spurious divergence at T c . The broken minimum exists only after the triplet 3d mass parameter (µ 2 Σ,3 ) becomes negative (when it is positive the value of the effective potential in the triplet phase is imaginary), at which point the triplet minimum immediately becomes the global one: the critical temperature is therefore erroneously identified -at all orders -with the condition that the triplet 3d mass parameter vanishes. The divergence at the critical temperature arises from a logarithm of the triplet mass parameter as it goes through zero. At higher orders in this expansion it is expected that stronger IR diverges will occur.
(c) Soft→supersoft EFT, mixed method: While this method is gauge-independent, one can see from Figs. 8c and 9c that it yields spurious divergences at NLO and N2LO for some RG scales. This problem arises at the edge of the range of temperatures where the LO result ceases to exist, resulting in ill-defined behaviour for the condensate. This problem was further discussed around Eq. (4.19) in Sec. 4, where the strict method was proposed as a general solution.
(d) Soft→supersoft EFT, strict method: Finally, this approach resolves all the theoretical problems encountered by the other approaches, and seemingly converges towards the lattice results with impressive accuracy. For the more weakly coupled BM1, the LO result in the supersoft EFT already agrees well with the lattice, and higher orders lead to small improvements, especially noticeable in the broken phase. However, it is for the more strongly coupled BM2 where this method clearly outshines the others, especially in the vicinity of the phase transition, where this method converges towards the lattice without spurious artefacts.
We highlight that despite the respective successes of different EFT expansions in predicting T c and the value of the triplet condensate as a function of temperature, in all approximations perturbation theory incorrectly predicts the character of the transition for BM1. The transition is a smooth crossover, as measured on the lattice, whereas perturbation theory predicts a weak first-order transition. In BM2 the transition is actually first order, and perturbation theory predicts it correctly.

Two-step transition
Finally, we turn to a two-step phase transition, for which we suggest a novel prescription in terms of two separate EFTs for consequent transitions. From the previous section, we know that for the first transition to the triplet phase we need to use the supersoft EFT. For the second transition, we have multiple options. In the plots of the critical temperature, the horizontal band is the lattice result together with its statistical uncertainty, and the bars show perturbative results at each order in the EFT expansion with uncertainty due to varying the RG scale.
First, there is a possibility that the second transition happens at the soft scale. As discussed in Sec. 3.3.3, such a transition can be induced by the semisoft scale in addition to the hard scale. We depict the results based on this approach in Fig. 10. 15 This figure depicts both condensates as well as the critical temperature for the triplet to Higgs transition, at each order in the expansion. The plots indicate some degree of convergence, yet even at the highest orders we have computed they do not provide striking agreement with the lattice results. This signals either that still higher order contributions should be included (especially since several  Figure 11: Convergence of supersoft strict approximations, where in both cases ∆T has been defined relative to the higher temperature transition. Note that the higher temperature transition for BM1 is a crossover, so T c here corresponds to a pseudo-critical temperature, defined as the peak in the susceptibility for the scalar condensates. RG improvements kick in only at N6LO) or that the assumption that the transitioning fields live at the soft scale is not correct, for the benchmark points in question.
Hence, we next test the possibility that the transitioning fields of the second transition live at the supersoft scale. This requires two separate EFTs for each of the triplet and Higgs phases. Comparison of the N2LO result to lower orders and convergence of the expansion is depicted side-by-side in the triptychs of Fig. 11. These triptychs of plots indicate reasonably good convergence, from LO results which only rough accord with the lattice data, each additional order yields closer agreement. We observe that at N2LO our novel perturbative computation of the scalar condensates using EFT expansions at the supersoft scale provides a striking correspondence with lattice results. While the convergence is clearer in BM1, it is also the case in BM2 which has a portal coupling more than twice as large, albeit the critical temperature for the second transition is further away from the lattice result. In analogy to Fig. 1 in Sec. 1.1 we summarise our analysis in Fig. 12 which recollects the N2LO result of Fig. 11. This plot indeed demonstrates the key result of this article: EFT expansions resolve all theoretical blemishes that have haunted perturbative predictions of the past, and while doing so, provide results that are not far away from those obtained on the lattice.

Discussion
In the past few decades, a rich patchwork of perspectives and insights have been developed regarding the reliability of perturbation theory to describe cosmological phase transitions. In one thread of the inquiry, a range of thermal hierarchies of scale were identified, and corresponding resummation schemes to correctly account for them. The early development of high-temperature dimensional reduction was based around the hard, soft and ultrasoft scales [17,19], yet another scale, the supersoft scale, was identified as central to first-order phase transitions [15]. In a separate thread of inquiry, a range of different perturbative methods were developed for the study of equilibrium thermodynamics, from direct minimisation of the thermal effective potential to strict ℏ-expansions. Concerns were raised that the direct minimisation method led to gauge-dependent results [42], while strict ℏ-expansions appeared to lead to IR divergences [42]. Concern about gauge dependence was then later revived in [79], further inspiring [25]. In recent years, significant progress towards resolving the aforementioned problems was made in Refs. [26-30, 32, 35, 36, 54, 61, 80, 83]. In this work at hand, we have unified, generalised and expanded most of this progress to a revised EFT framework for equilibrium thermodynamics that builds from the dimensionally reduced 3d EFTs [17,18,[21][22][23]68], but also consistently applies strict power-counting expansions in perturbation theory [30].
Concretely, we have simultaneously tested both these threads of inquiry, and have found a consistent resolution to all the concerns in terms of self-consistent perturbative EFT expansions. The results of recent lattice Monte-Carlo simulations at two benchmark points in the real-triplet extended SM [37] have formed the bedrock of these tests. This has allowed us to obtain an unambiguous measure of the error in different perturbative approaches.
Our results, summarised in Fig. 12, attest to the correctness of a particular perturbative approach, which is both theoretically consistent and numerically reliable. This approach is rather simple, and in hindsight obvious. It is the following: 1. Successively integrate out UV modes, starting from the hard scale, and working towards the IR, and stopping when one meets the mass scale of the transitioning fields.
2. The mass scale for the transitioning fields can be identified by power counting, applied to the tree-level potential of the EFT. If there is an apparent second-order phase transition at tree-level within this EFT, then more modes must be integrated out.
3. Carry out strict perturbative expansions within the EFT for the transitioning field, ensuring to remain within the region of validity of the EFT.  [37]. The supersoft EFTs are constructed separately in each broken phase, with the triplet field treated as supersoft within the triplet phase and the Higgs field treated as supersoft within the Higgs phase. This approach yields gauge-invariant results, in good agreement with the lattice, and a significant improvement over previous perturbative approaches (see Fig. 1).
underlying expressions for effective potentials and for thermodynamic quantities are astonishingly simple -excepting complicated, yet closed form, expressions for LO broken mimima in cases of radiatively generated barriers -and the striking agreement with lattice results highlights that the underlying physics is well captured in perturbation theory. 16 The good agreement of perturbation theory with the lattice shown in Fig. 12 should be contrasted with that of Fig. 1. The crucial difference is that in Fig. 12 the supersoft scale has been correctly identified as the energy scale of the transitioning fields. Our results align with [30,61], which also compared perturbation theory to the lattice, and came to similar conclusions. Yet, in this work at hand we have for the first time applied these developments to a BSM theory, where the studied phase transition pattern is more complicated and leads to a rich chain of EFT setups. Indeed, in the course of this study, a number of further technical manoeuvres have been identified. We have shown how, when one is interested in thermodynamic observables away from the critical temperature, it is advantageous to reexpress quantities in terms of deviations from the critical temperature ∆T = T − T c . This resolves a problem of the existence of the LO result required by a strict expansion, and 16 Yet we emphasize again that perturbation theory cannot separate weak first-order transitions from crossovers, or describe purely non-perturbative phenomena related to phase transitions, such as condensation of monopoles [84].
underlies the difference between our mixed and strict approaches. We have also discovered that a new scale between the hard and soft scales, here dubbed the semisoft scale, arises naturally in strong two-step first-order phase transitions.
The strict perturbative EFT expansions presented in this work can still be extended by computing the final perturbative orders that are available before crashing against the Linde problem of non-Abelian gauge theories at four-loop order. Computing these final orders requires three-loop vacuum diagrams, and provides an intriguing future challenge analogous to that achieved in hot QCD [19,85,86]. Yet another, different kind of challenge will be to incorporate the presented EFT expansions to parameter space scans of phenomenologically interesting models. Herein the challenge lies in the implementation: different EFTs may be required to study different parameter points, and even in a single parameter point there can be several EFTs in different temperature regimes. This issue has also been raised in the recent Ref. [34] and is further discussed in Appendix B. This reference indeed discusses many of the same ideas as detailed in our work at hand, yet our computation includes concrete applications to a BSM model, as well as comparison to lattice data.
Finally, while we have limited ourselves to the study of equilibrium thermodynamics, EFT expansions are expected to carry over to the perturbative study of other properties of first-order phase transitions, such as the bubble nucleation rate and the bubble wall speed, as well as the sphaleron rate. For the bubble nucleation rate, Refs. [26,27,29,33] have in fact already used the approach proposed here. On the other hand, transferring what we have learnt here in this article to tunnelling will be challenging. While it is possible to utilise different EFTs in different phases for computing the free-energy of homogeneous phases, for bubble nucleation this must be generalised to non-trivial paths in field space. The formulation of the EFT description of bubble nucleation for such a two-step transition, warrants dedicated future studies.

A EFT expansions with doublet and triplet fields
In this appendix we present explicit expressions for EFT expansions of the effective potential for the real-triplet extended Standard Model. We can immediately read the effective potential for the soft scale EFT from [37]. Therein, the presented "ℏ-expansion" matches the soft-scale strict expansion of Sec. 4. Our task is then to compute the effective potentials at the supersoft scale EFT, as well as the soft scale EFT in the presence of fields at the semi-soft scale.
For this computation, we can read from [37] all the two-loop diagrams we need, which we reorganise into EFT expansions. In this section, all masses m i are mass eigenvalues in the 3d EFT, and are not to be confused with physical pole masses. We use Landau gauge throughout.

A.1 Supersoft EFT
In this section, we use the following notation for the supersoft scale effective potentials in the EFT expansion where each order beyond V 0 at LO is suppressed by ε ∼ √ ε soft relative to the previous order; see the discussion around Eq.(4.26). Notably, V 1 is identically zero.
Supersoft Higgs doublet For simplicity, we start with the 3d EFT of the Standard Model with SU(2) and U(1) gauge fields. We parametrise the Higgs doublet in the 3d EFT as i.e. we compute the effective potential in terms of real background field v > 0. The leading order contribution reads where m W = 1 2 g 3 v and m Z = 1 2 g 2 3 + g ′ 3 2 v. The two-loop contribution from gauge fields and ghosts reads )v 2 ln(2) This can be obtained from the result of Ref. [37] by dropping triplet contributions and by setting scalar masses to zero inside two-loop diagrams. The resummed mass eigenvalues for the scalar fields -Higgs (h) and Goldstone bosons (G ± , z) -read where the Goldstone square mass eigenvalue is triply degenerate. These expressions can be derived starting from a general potential written in terms of the gauge-invariant bilinear ϕ † ϕ, and expanding to quadratic order in fluctuations, where prime denotes a derivative with respect to v, and we have shown only the relevant bilinear terms. In terms of resummed masses, the one-loop supersoft contribution reads SU(2)+Higgs The above expressions for the supersoft effective potential becomes much more compact when the U(1) gauge sector is decoupled, i.e. in the limit g ′ where we have substituted m W = 1 2 g 3 v. This result has been previously obtained in Ref. [30].
Supersoft doublet with a soft triplet Next, we include the effect of a soft triplet, i.e. the triplet remains at zero background field, and is integrated out together with the gauge fields. One-loop triplet contributions arise at leading order, so that and m 2 Σ = µ 2 Σ,3 + 1 2 a 2,3 v 2 . Lagrangian parameters are defined in Eq. (1.1). The triplet squared mass is triply degenerate, since the neutral and charged triplet have equal masses for vanishing triplet background field. Note, that here we have assumed that also µ 2 Σ,3 is soft: this is relevant for one-step transitions directly to the Higgs phase in the presence of a soft triplet that can enhance the transition strength, but also for the second step of a two-step phase transition. In practice this requires that there is enough supercooling between the two transitions, and the magnitude of µ 2 Σ,3 (which can and often will be negative) can increase parametrically from the supersoft to the soft scale after the first transition. In this case, the expression for the minimum of the LO potential of Eq. (A.10) becomes utterly complicated, yet it can still be found analytically. Below we comment on the case where µ 2 Σ,3 is still at the supersoft scale, albeit m 2 Σ is soft. Resummed scalar masses related to the supersoft doublet read V 3 is of same form as before in Eq. (A.8), but with the above modified supersoft masses. Next we have where the triplet contributions read If µ 2 Σ,3 is supersoft despite m 2 Σ being soft, we can account for the effect of the triplet 3d mass parameter as a perturbative mass insertion in the matching. In this case, the LO potential has a simple expression 2,3 . The broken minimum of this potential has the simple expression Evaluating expressions in this LO minimum in strict expansions, results in relatively simple, closed-form expressions for many quantities. The resummed Higgs and Goldstone masses read Triplet NLO contributions are obtained by replacing Eq. (A.15) as where M 2 Σ ≡ 1 2 a 2,3 v 2 , and we have added the triplet one-loop diagram with a single IR mass insertion, which contributes at O(ε soft V 0 ) to the soft-to-supersoft matching.
Supersoft triplet In this section, we assume that the triplet is supersoft, whereas the Higgs is soft. For simplicity, we start with the case of a sole triplet, and only after add contributions from the Higgs. Note that the triplet is not charged under the U(1) gauge group. We parametrise the triplet field as where the triplet background field is denoted by x > 0. The leading order effective potential reads where m W = g 3 x, i.e. the W-boson contribution is resummed together with a tree-level potential. The Z boson is massless in the triplet phase. Two-loop diagrams with W bosons and ghosts yield an extremely simple result in which we have used m W = g 3 x. Note that no logarithmic terms arise: in the sole triplet case the triplet mass parameter does not run at this order. 17 In the supersoft scale EFT triplet masses are those resummed by the soft gauge-field contributions. For neutral and charged triplets in the broken triplet phase we have Consequently, the one-loop triplet diagrams yield where m 2 ϕ = µ 2 h,3 + 1 2 a 2,3 x 2 is the quadruply degenerate mass squared eigenvalue for the Higgs field. Again, we have first assumed here that µ 2 h,3 is soft. Resummed masses for neutral and charged triplets get new contributions accordingly, (A.30) 17 In the 3d counterterm δµ 2 Σ,3 presented in [37], the g 4 3 contribution comes solely from the Higgs doublet loops.
The NLO potential is the one-loop contribution with semisoft masses N2LO vanishes, V 2 = 0, and the result at N3LO is given by one soft-mass insertion of the triplet to the one-loop bubble diagram . (A.37) All higher order terms include contributions from the soft EFT, and we highlight these contributions separately below, in addition to matching contributions from the semisoft scale. At N4LO we get Here the soft EFT contribution results simply from one-loop bubble diagrams with unresummed masses. The matching contribution comprises of two-loop diagrams without soft mass insertions. At N5LO there are contributions from the triplet one-loop bubble with two soft mass insertions, as well as a contribution from within the soft EFT, The soft EFT pieces at N5LO result from resummations of V 1 . That is using these resummed masses, inside one-loop bubble diagrams, and then re-expanding in ε. We higlight that our result for the soft EFT expansion of the effective potential is RG invariant at the order we truncate our computation.
Triplet phase The effective potential in the triplet phase has an expansion where ϵ ∼ g/π. Mass eigenvalues read m 2 W = g 2 3 x 2 ∼ ( √ gπ)T 2 , and a quadruply degenerate m 2 ϕ = µ 2 h,3 + M 2 ϕ , where the semisoft piece is M 2 ϕ ≡ 1 2 a 2,3 x 2 ∼ ( √ gπ)T 2 and µ 2 h,3 ∼ g 2 T 2 is soft. Soft masses are m 2 Σ 0 = µ 2 Σ,3 + 3b 4,3 x 2 and m 2 Σ ± = µ 2 Σ,3 + b 4,3 x 2 . In analogy to the computation for the Higgs phase result, we get the following expressions, At higher orders, in analogy to the computation in the Higgs phase, we get  at N4LO and N5LO respectively. In particular, the soft EFT pieces at N5LO result from resummations of V 1 using Again, our result is RG invariant at the order we truncate.

B Direct minimisation
In the past few decades, direct minimisation of the real part of the thermal effective potential in Landau gauge has solidified itself as the standard meta in studies of cosmological phase transitions. Most studies resort to a one-loop approximation, as for generic models this has  Figure 13: Convergence of the direct minimisation method within the supersoft EFT for both transitions. Already lower order results are fairly close to lattice results, and convergence is clear. Notably in BM2, the first transition is seemingly weaker at higher orders compared to LO, as indicated by the size of the discontinuity in the triplet condensate.
an explicit and relatively simple formula, however the downside is that it suffers from rather large theoretical uncertainties. When going beyond one-loop accuracy, the functional form of the effective potential becomes much more complicated and direct minimisation becomes numerically expensive. On the other hand, strict EFT expansions are numerically cheap to evaluate, once higher order corrections to the effective potential are known, as they can be obtained by straightforward Taylor expansions around leading order results. All higher order corrections are obtained simply by evaluating these expressions numerically. In order to provide a comparison of these approaches, in this appendix we present results obtained by direct minimisation of the real part of the Landau-gauge effective potential, both for the supersoft and soft scale EFTs. First, Fig. 13 showcases the supersoft EFT and the direct method for both transitions at both benchmark points. Despite residual gauge dependence and the need to discard the imaginary part of the potential in minimisation, the result for the value of the condensate in the broken phase and for the critical temperature align with the lattice data fairly reliably. Notably, already the leading order result is reasonably good, yet N2LO is even better,  Figure 14: Convergence of the direct minimisation method for the hard→soft EFT. Results at lower orders are seemingly far off, while at N2LO, or two-loops, the perturbative result agrees fairly well with the lattice data despite its theoretical hiccups related to imaginary parts and gauge dependence. However, note that around the higher critical temperatures there are a number of perturbative data points missing. This is due to the numerical minimisation algorithm failing with default tolerance and method arguments, and is a common issue arising when directly minimising the real part of effective potentials at higher loop orders.
indicating convergence. On the other hand, given the computational cost of direct minimisation, there should be no practical reason why not to upgrade this computation by a strict expansion, the results of which are shown in Fig. 11. Finally, we present triptychs of the convergence of direct minimisation for the hard→soft EFT in Fig. 14 for both BM1 and BM2 and for both transitions. We observe that while LO (tree-level) and NLO (one-loop) leave much to hope for, the result at N2LO (two-loop) agrees fairly well with the lattice results for both benchmark points, for both the value of the condensates and the critical temperatures. This observation provides some support for using the direct minimisation method within the hard→soft EFT in parameter-space scans of BSM theories, as long as the computation is performed at two-loop order. Despite the lack of theoretical robustness, this method allows one to scan wide regions of BSM theory parameter space in a single EFT, in contrast to EFT expansions which require delicate usage of chains of EFTs, potentially even in a single parameter point. Yet in practice, one major downside of the direct approach at two-loop order is the computational cost of minimising complicated multivariate functions. This motivates pursuing the automation of EFT expansions, which are numerically significantly less expensive and furthermore theoretically sound.