A Critical Look at the Electroweak Phase Transition

The electroweak phase transition broke the electroweak symmetry. Perturbative methods used to calculate observables related to this phase transition suffer from severe problems such as gauge dependence, infrared divergences, and a breakdown of perturbation theory. In this paper we develop robust perturbative tools for dealing with phase transitions. We argue that gauge and infrared problems are absent in a consistent power-counting. We calculate the finite temperature effective potential to two loops for general gauge-fixing parameters in a generic model. We demonstrate gauge invariance, and perform numerical calculations for the Standard Model in Fermi gauge.


Introduction
The electroweak symmetry appears exact in the early universe, but not in our current day and age. As the universe expands and cools down the Higgs field develops a vacuumexpectation-value (VeV)-breaking the symmetry. There is a phase transition from a symmetric to a broken phase.
Although well established in the Standard Model, the electroweak phase transition remains elusive. It is, as yet, unknown when and how the transition took place; if it was violent, or calm; first-order, or continuous. Continuous transitions are rather innocuous compared to their first-order cousins. Indeed, a first-order phase transition is a turbulent and violent affair that likely has far-reaching consequences for the universe's development. Such transitions are part and parcel for understanding the observed matter-antimatter asymmetry [1]. Furthermore, gravitational waves from a strong phase transition reverberate throughout the universe and might be picked up by next-generation experiments [2]. Describing these phenomena goes hand-in-hand with understanding phase transitions.
Both perturbative and lattice methods accomplish this task; both methods with their fair share of virtues and vices.
On the perturbative side, the name of the game is the effective action. Calculations are carried out with Feynman diagrams-both quantum and temperature effects are included in this approach. But problems loom around the corner.
It is unclear if perturbative methods can at all be trusted [3,4]. There are ample issues related to gauge invariance [5]; appearance of infrared divergences [6]; not to mention a breakdown of perturbation theory itself [7]. By charging headlong one might even erroneously conclude that first-order phase transitions disappear for some gauge choices. And the perturbative expansion is likewise dicey. It can, and does, break down. Consider a scalar theory with quartic coupling λ. Finite temperature calculations include diagrams d N , at loop order N , scaling as d N ∼ T (λT 2 ) N −1 . So loops are not suppressed for large temperatures [8]: d N /d N −1 ∼ λT 2 ∼ 1 for T 2 ∼ λ −1 . What's more, the phase transition occurs at these very temperatures-as we discuss in section 2.3. So perturbation theory slowly but surely breaks down. This is not surprising in itself. After all, phase transitions occur when loop corrections overpower tree-level terms, which calls the loop expansion into question.
Yet it is long known how to alleviate these problems [7][8][9]: a resummation is needed. As particles interact with the thermal bath they receive a thermal mass [7]-a reorganization of the perturbative expansion around this effective mass improves convergence. Although the need for resummation is established, there's no consensus on the implementation. There are a number of conflicting strategies [5,8,[10][11][12]. More often than not resummations are gauge dependent.
But perturbation theory is not the only avenue. Indeed, lattice calculations are quite good at describing phase transitions. Though not without their own share of issues. Lattice simulations are resource expensive-large parameter scans are, as yet, unfeasible. Still, lattice is well-suited at studying single models. The Standard Model's phase structure has indeed been investigated via lattice calculations [13][14][15][16].
Perturbative calculations are on the other hand computationally cheap; so they are felicitous for studying complicated models. In this paper we propose a gauge invariant resummation. We develop robust perturbative techniques for describing phase transitions. These techniques are gauge invariant and aren't stymied by IR-divergences. Sticking to a strict power counting scheme is integral; a proper perturbative expansion-with powers correctly accounted for-is inherently gauge invariant. Our method is an amalgamation of (i) the early work of Arnold and Espinosa [8], and (ii) the gauge invariant methods developed by Laine [17] and emphasized by Patel and Ramsey-Musolf [5]. We restrict ourselves to observables at the critical temperature in this paper, and leave tunneling for the future.

The powers of perturbation theory
Perturbative calculations are organized in powers of a small quantity. This might be a collection of couplings, or a ratio of energy scales. All terms must, at a given order, be included. The consequences of forgetting terms are dire-including gauge dependence and exasperating divergences. The same holds when calculating the effective potential. This section discusses subtleties and dangers of perturbative expansions. We introduce a systematic way to treat the breakdown of the "naive" loop expansion. We apply these considerations to the effective potential and show how and why a proper power counting is useful. The issues and concepts we discuss are well known, but we introduce our own notation.

Power and loop counting
In order to illustrate how a perturbative expansion might break down, and how it might be fixed, some terminology is in order. For example, in a standard loop expansion an observable A is typically expanded as A = A 0 + κA 1 + κ 2 A 2 + . . . (2.1) with κ denoting the number of loops. Yet all is not fine and dandy. For if a coupling is large, the expansion might not be justified at all. And there are situations where calculations-even in weakly coupled theories-are not organized in loop powers.
So we better make a clear distinction between power-and loop-counting. To that end, introduce a new power counting parameter that better represents the actual sizes of terms: ħ h. As a side-note, ħ h is not related to the reduced Planck constant which goes by the same symbol. We choose ħ h as the power counting parameter only to be congruous with earlier papers [5,17].
If the loop expansion is applicable, ħ h and κ are equivalent: However, there might be terms in A n scaling with negative powers of ħ h. Consider a toy example, where A n = a n + b n /ħ h n−1 if n ≥ 2. The expansion is and diagrams from all loop orders are intertwined. If this is the case, a resummation is appropriate. These ideas can be made lucid through a few examples.

Gauging the problem
The effective potential is in perturbation theory calculated order-by-order according to some power counting scheme, with ħ h denoting the aforementioned power: The idea is to find the global minimum φ min , which then gives the physical energy density: V min ≡ V (φ min ). The "standard" approach finds φ min by minimizing V (φ) numerically. But this procedure is problematic and gives a gauge dependent V min -which doesn't make sense for a physical observable. The effective potential at an arbitrary field value is not a physical observable, which the Nielsen identity makes glaringly clear [18], (2.5) ξ is here a gauge-fixing parameter and C(φ, ξ) is a calculable function known as the Nielsen coefficient. This is a non-perturbative statement describing the effective potential's gauge dependence. The equation suggests that V min is gauge invariant, but that φ min necessarily depends on ξ. There needs to be a delicate cancellation between the gauge dependence of φ min and V for V min to be gauge invariant. Why is it then not sufficient to minimize the potential numerically? The devil is in the details of perturbation theory, and a fiery analogy might be appropriate. Consider a particle's pole mass as calculated in perturbation theory, This is an implicit equation for m 2 P . Solve it by further expanding m 2 P on the right-hand side, according to which gives the well-known result Comparing the two equations (2.7) and (2.8) order-by-order in ħ h, we deduce m 2 0 = m 2 and m 2 1 = Π 1 (m 2 ). Likewise, as emphasized in [5], φ min must in turn be found order-byorder in ħ h, . . .
The minimum can then be plugged into the effective potential to give the physical and gauge independent energy density Notice how all terms are expressed at φ 0 . This is expected since the expansion is organized around the leading order value. Although consistent power-counting schemes are gauge independent, the appropriate counting is not determined a priori.

Coleman-Weinberg
A well known application of the effective potential is due to S. Coleman and E. Weinberg, where they establish the mechanism of quantum-generated spontaneous symmetry breaking [19]. They considered scalar electrodynamics without a scalar mass term: (2.14) There is no symmetry breaking at tree-level. But the symmetry can be still be broken by quantum effects-the Coleman-Weinberg mechanism.
Explicitly, expand the potential as usual sub-leading corrections are given by scalar and photon loops How can the symmetry be broken by quantum corrections, which, after all, are suppressed in ħ h? For this to happen, the 1-loop terms must compete with the tree-level terms. This indicates that λ must be small, λ ∼ e 4 , compared to the standard (loop) counting λ ∼ e 2 . This is accounted for by systematically counting lambda as ħ h: λ → ħ hλ, implying (2.17) Spontaneous symmetry breaking is now possible. Driven by a quantum effects. But the situation is peculiar at higher orders. In Fermi gauge there are diagrams d N at N loops scaling as d N ∼ e 4N /λ N . The new power-counting λ → ħ hλ ∼ ħ he 4 indicates that these terms are all of the same order. All of these terms must be included-they must be resummed. The authors in [6] showed how this resummation in concert with the ħ hexpansion gives a gauge invariant result. The resummation grabs the relevant terms from each loop order and organizes them in a gauge-invariant manner.
(2. 19) This paper is concerned with temperatures much larger than the masses: T 2 X . In this case the n = 0 zero-mode is distinctly different from n = 0 finite modes-the zero-mode propagator does not depend on temperature. Note that zero-modes only appear in boson propagators.
It is well established that the loop expansion breaks down in finite temperature field theory for high enough temperatures. These temperatures are large enough to invalidate the loop counting scheme. The worst eggs are the diagrams known as daisies [7]. They are made up of a soft momentum (p T ) inner loop, strung together with hard (p ∼ T ) selfenergy insertions. Inner loops contribute T and each hard self-energy contributes λT 2 , for some coupling λ. An N -loop daisy d N then scales as d N ∼ T (λT 2 ) N −1 , and isn't suppressed compared to the (N − 1)-loop daisy, to wit (2.20) So perturbation theory breaks down for temperatures of order T 2 ∼ 1/λ. The well-known resolution to this problem is to perform a resummation in which bosons acquire a thermal mass. This removes all the problematic terms and replaces them with a single resummed term. Or rather, gathering up all daisies effectively resums the mass. Now consider the implications for the finite temperature effective potential. Leading temperature corrections show up at one loop. We are interested in large temperatures, so take a high-temperature expansion for granted: T M Ψ (φ) for all fields Ψ. In the high temperature expansion at one loop there are terms that contribute as T 2 , T, T 0 , T −2 , . . ., which we denote as¹ Note that V min is evaluated in the same way as in the zero temperature ħ h-expansion. First calculate φ min perturbatively and then evaluate the potential at φ min . The difference is that φ min depends on the temperature: To untangle the notation a bit, consider the T 2 correction at ħ h and ħ h 2 : This expression is gauge invariant order-by-order in ħ h-as we have confirmed to two loops.
1 Here and in the following we always discard terms that are independent of φ.
For a general potential the naive leading-order contributions are where both loop counting, with κ, and naive power counting, with ħ h, are included. The situation is disparate at high temperatures. The leading behaviour is set by the classical potential V 0 , and the (largest) loop term is given by ħ hT 2 V 2 1 . The loop term can only compete with the classical potential for temperatures of order T ∼ 1/ ħ h, which begs for a reshuffling of the perturbative expansion-analogously to the reshuffle in the Coleman-Weinberg model discussed in section 2.3.1.
Making this power-counting manifest by rescaling T → T / ħ h, the new expansion is Higher loop terms are now as important as lower loop ones-the harbinger of a resummation.

Thermal resummations and power counting
Close to the phase transition temperature the potential takes a form akin to equation (2.25). The minimum background energy is found by minimizing the potential order-by-order in ħ h: φ min = φ 0 (T ) + ħ h 1/2 φ 1/2 + ħ hφ 1 + . . . The minimization conditions are O ħ h 0 : (2.27) . . . The short-hand notation ∂ ≡ ∂ φ is used extensively to avoid clutter. Note that the leading order VeV φ 0 (T ) is temperature dependent, and terms starting at ħ h 1/2 get contributions from all loop orders.
The energy is This result should be gauge invariant if the power counting is consistent. As it stands it is only possible to check gauge invariance if an infinite number of diagrams are included. Serendipitously enough, it's possible to resum all terms. It is instructive to discern why a resummation is necessary in the first place. The new counting T ∼ 1/ ħ h implies that the leading-order result is an amalgamation of 1-loop terms with tree-level ones. So the inverse propagator of a particle with squared mass x is enhanced, and loop-corrections are of the same order as tree-level masses.² This is akin to the familiar hard thermal loop resummation [7], but with minor differences. First of all, nothing has been said about which propagator lines should be resummed. There have been arguments in the past [8,20] to only resum zero-modes; because the worst divergences are removed, and diagrams are not double-counted. We take a different approach, and let gauge invariance guide the way. We have confirmed that no linear-in-φ term is generated in Abelian Higgs and the Standard Model, which is a nice consistency check [20].
There's an infinite number of terms at NLO, These are the most divergent pieces in the daisy diagrams. And they all contribute to the resummation of the zero-mode: Scalars and 3D-longitudinal gauge boson have been resummed in V 1 1 . Recall that subleading terms must be evaluated at the temperature dependent minimum, As discussed in [5], at one loop all gauge dependence manifests itself in the masses of the Goldstone bosons. So the gauge dependence only cancels if all Goldstone masses are identically zero. Thus in a theory with a standard loop counting the 1-loop potential would be evaluated in the temperature-independent tree-level minimum-where resummed Goldstone masses are non-zero. When Goldstone masses are resummed according to the process above, they are zero in the new minimum. And hence the NLO correction to V min is gauge invariant. This is part and parcel of the ħ h-expansion [5].
Moving on to NNLO, some novel patterns appear. Start by considering the temperature independent 1-loop term ħ hV 0 1 . In the Arnold-Espinosa approach [8] this term is not resummed. Yet there are good reasons to resum it. First, there are terms at ħ h coming from all loop orders. Second, V 0 1 isn't gauge invariant without a resummation. The reason is the same as for the ħ h 1/2 term.
But there is another reason for resumming this term. Two-loop terms of the form T 2 V 2 2 have two different origins. The first is purely due to zero-modes (soft momenta) and cannot be removed by resumming. The second part comes from a mix of zero-and finitemodes. Before terms were reshuffled, the ħ h-expansion contained terms ∼ T 2 φ 0 1 φ 2 1 ∂ 2 V 0 . These terms are of the same form: a finite-mode contribution φ 0 1 , and a zero-mode con-tribution φ 2 1 . But these terms are washed away-the temperature scaling pushed them higher in the expansion.
It turns out that resuming V 0 1 is equivalent to including the aforementioned ħ h terms lost by the scaling. To wit, resuming a mass X in V 0 1 demands a subtraction to avoid overcounting: with similar subtractions at higher orders. To sum it up, V 1 0 φ 0 (T ) is gauge invariant, and so are the remaining 2-loop terms after subtracting diagrams.
In this way all the gauge dependence of T 2 V 2 2 is cancelled in two steps. The resummation of V 0 1 removes the first chunk. And resumming at two loops (T 2 V 2 2 + . . . → V 2 2 ), together with the ħ h expansion, removes the last bit since Goldstone masses vanish at φ 0 (T ).
To be clear, we advocate that the scalar masses should always be resummed, beyond their contribution to the leading order potential. This is demanded by gauge invariance. Gauge bosons are another matter, because only 3D-longitudinal zero-modes have a large self-energy. Hence only zero-modes of vector bosons should be resummed. We give an extended discussion about how to resum vector boson masses at higher orders in appendix C.

Phase transitions
Whereas the previous section delineated how the perturbative expansion of the effective potential is reshuffled with the scaling T ∼ 1/ ħ h, this section applies these results to phase transitions, both first-and second-order.
To make the discussion lucid, focus on the generic potential

Second-order transition
Consider first a second-order transition. With the scaling T ∼ 1/ ħ h the energy is .
The leading-order term V 0 + T 2 V 2 1 determines the temperature dependent VeV φ 0 (T ). Terms in T 2 V 2 1 are gauge invariant and are of the form ∼ e 2 φ 2 T 2 for some coupling e [5]. So all that changes for finite T is m 2 → m 2 eff (T ). The transition occurs at the temperature where m 2 eff (T ) changes sign: m 2 eff (T 2nd ) = 0. This is a second-order transition.

First-order transition
Let's for a moment forget everything about proper power-counting and just try to naively describe a first-order transition, where the minimum abruptly changes from non-zero to zero for some temperature T c . This requires a barrier to develop between the two minima.
To be concrete, consider a high temperature expansion in the Abelian Higgs model. For high temperatures the potential is approximately Following [8], these various terms have to balance each other for a barrier to develop. The balance occurs if So does this scaling always work? No. It depends on the couplings: vector bosons' thermal masses, ∼ e 2 T 2 , dominate tree-level ones if λ ∼ e 2 , which would break any semblance of a power-counting.³ A counting like λ ∼ e 4 -as in the Coleman-Weinberg model-is likewise dicey. To wit this counting implies eφ ∼ T which invalidates the high-temperature expansion. To let λ scale as higher powers of e will only worsen the problem, and lower powers than 2 will similarly break the perturbative expansion. This leaves only one option [8], So we should really be counting powers of e, and be fastidious about the power-counting.
In the end the first-order scaling is a hybrid between a Coleman-Weinberg-like scaling (pushes terms up in order) and the second-order scaling (drags terms down to lower orders). There will be infinite towers of diagrams at each order in the perturbative expansion, just as for the second-order scaling. Though note that scalar masses now scale differently. For example, the resummed Goldstone mass scales as G ∼ m 2 eff (T ) ∼ ħ h 1/2 . This implies that previously sub-leading Goldstone self-energy terms of order Tħ h 1/2 must now be resummed. So resummed scalar masses are X = X + T 2 Π 2 X + T Π 1 X , where only leading order terms are included in Π 1 X . This is quite natural since V LO includes order T and T 2 terms; inherited by scalars through This does not apply to gauge bosons since their masses scale as before.

First-order counting and gauge dependence
The above discussion disregarded everything that had to do with gauge symmetry and further complications from the power counting. So it may not be surprising that a naive application of this method is gauge dependent. The effect is particularly transparent in R ξ gauges.
The R ξ effective potential is schematically where the Goldstone's zero-mode has been resummed. The development of a barrier required for a first-order transition is driven by terms proportional to e 3 T φ 3 . Note that these terms vanish for ξ = 3 2/3 . So the gauge dependence is no paltry effect. Not only does the potential depend on ξ, the very nature of the phase-transition is extremely sensitive of ξ.
The situation is alleviated with a proper power-counting. Consider the first-order transition scaling λ ∼ e 3 , m 2 eff (T ) ∼ e 3 T 2 , T ∼ 1 e . A new minimum develops when the quartic term competes with the mass term: φ ∼ T . Now, the Goldstone mass is of order G ∼ e 3 T 2 , while the photon mass is of order e 2 φ 2 ∼ e 2 T 2 . This means that the gauge dependent terms (to leading order) cancel, leaving Gauge dependent terms are subleading. What's more, gauge dependent terms are evaluated at φ 0 (T ), and by definition vanish after a resummation: This is another sign that first-order transitions can only be described perturbatively if e 2 λ.

Details of the perturbative expansion
Due to its numerous appearances, it is felicitous to use e instead of ħ h for counting powers. So e serves bilaterally as a power and a constant-a powerful constant indeed. Gauge bosons scale as Z ∼ e 0 , and scalars as H, G ∼ e. In the Standard Model for example e ∼ α W ∼ 0.1. The VeV scaling (φ ∼ T ∼ e −1 ) implies that the leading-order potential scales as V 0 ∼ λφ 4 ∼ e −1 . Next-to-leading order terms come from T 2 V 2 2 and V 0 1 (with scalars and powers of lambda pushed to higher orders); these terms scale as e 0 . Cracking on, NNLO is solely due to scalar T V 1 1 terms.⁴ N 3 LO goes as e and contains terms from T V 1 2 , T 2 V 2 2 , and T 3 V 3 3 . The potential and VeV are Where φ min is calculated order-by-order in e. Mark that a derivative with respect to φ adds a factor of e: ∂ ∼ e. So 4 Technically there are terms from T 2 V 2 2 , but these all cancel against resummation subtractions. implying . . .
Finally, the extremum energy is Schematically, a gauge boson Z and its 3D-longitudinal mode Z L , and scalars H, G, contribute to the different orders of the potential as where κ denotes loops. 2-loop calculations suffice to calculate the potential to NNLO. Now for the gauge dependence. Some features are quite clear up to NNLO. Scalar masses are determined from V L O , and all terms are evaluated at φ LO : Goldstone masses are zero, which removes most of the gauge dependence. Yet the expansion of the potential, . . , is in Fermi gauge only correct when ξ = 0. The discrepancy comes from ξ dependent terms, formally starting at e −1/4 . These terms all scale with some negative G power. When λ ∼ e 2 they are removed by the ħ h expansion. But when λ ∼ e 3 they are cancelled by a combination of resummation subtractions and the ħ h expansion. Since these ξ dependent terms always cancel among themselves, we've left them out of the expansion of V (φ). Although, these terms are relevant when explicitly checking gauge invariance. We also want to caution the reader that completely new divergences, compared to the second-order scaling, might appear when working with a finite ξ. For example, terms ∼ log G appear at intermediate steps at O e 0 ; though, these terms cancel in the end, albeit in a subtle way.⁵ The next check would come at N 3 LO and requires knowing the effective potential to three loops in a general gauge. It was only recently that the 2-loop effective potential was found for a general gauge [21]. So we'll cross that bridge when we come to it. In section 4 we calculate observables in the Abelian Higgs model and in the Standard Model with the first-order scaling outlined above, and point out possible complications and pitfalls.

A phase transition between two phases A and B occurs when
Or V A = V B for short. Since both V A and V B are gauge invariant by themselves, T c is guaranteed to be gauge invariant. There are two schools of thought on how to find the critical temperature. First, draw the phase-diagram for V A and V B , and change T until the two energies match up. This is the method proposed in [5], and it has the advantage that higher order corrections are easy to include. Though this method is gauge invariant, perturbative orders are again muddled.
An alternative is to find T c order-by-order in ħ h, as investigated by Laine [17]. But he noticed that an ħ h expansion for T c breaks down at ħ h 2 , and so the idea has long been dismissed. Yet this breakdown does not occur for all power-counting schemes.
Consider first the second-order scaling, . (3.23) The leading order energy vanishes in the symmetric phase: V A LO = 0. For the broken phase the leading-order energy is proportional to the Higgs' temperature-dependent mass: ,T c = 0 at leading order gives H(T )| T 0 = 0-causing problems at two loops. The 2-loop potential contains terms of the form T 2 log H(T ); since we'll expand around T 0 in the ħ h expansion these terms diverge and do not cancel between phase A and B. The expansion seems useless. In our mind this cements that the scaling λ ∼ e 2 cannot describe a first-order phase transition, and that the critical temperature is then simply determined from the leading order potential-with higher order corrections incalculable in perturbation theory. Yet a first-order transition is different. The effective Higgs mass H(T ) is finite both in the symmetric and broken phase, and no divergences can (naively) appear. Our explicit calculations, reviewed in section 5, show that first-order transitions are free of these subtleties. The expansion of the critical temperature is of the form T c = e −1 T LO + eT NLO + e 3/2 T N N L O + O(e 2 ). Derivatives with respect to T scale as e 0 when acting on G, H or V LO , and as e when acting on anything else. Denote the potential difference as ∆V (φ) ≡ V (φ) − V (0), whose expansion is . . .
Note that ∂ T ∆V L O scales as e −1 , which is why T NLO scales as e.⁶ The additional suppression (T N L O /T L O ∼ e 2 ) explains why corrections to T LO tend to be rather small, as seen in section 5. With T c found it is possible to calculate various observables at the phase transition. For example, the barrier height is Where ψ LO is the location of the leading-order maximum defined by (3.27) Temperature derivatives of φ L O (T ) can be rewritten as ∂ T φ LO = −∂ T ∂ φ V LO /∂ 2 φ V LO and similarly for higher orders. So everything boils down to an effective φ NLO : The new φ N LO automatically takes care of all implicit derivatives. So when expanding around T c we'll always use the temperature corrected φ NLO .

Thermodynamical Observables
Finding the critical temperature is all well and good, but there are a myriad of other observables. For example, the sphaleron transition rate is approximately controlled by v c T c [5], where v c is taken to be the minimum at the critical temperature. Larger values of this ratio indicate that sphaleron proccesses are suppressed enough after the phase transition to not erase any generated matter-antimatter asymmetry. Such phase transitions are strongly first-order if v c /T c ≥ 1. If this measure is to be physically meaningful, it has to be gauge independent. But the minimum, numerically found or perturbatively expanded, is gauge dependent and IR-divergent, and so a rather bad observable. Though there is a related observable, the scalar square vacuum expectation value |Φ| 2 . The scalar square expectation value is gauge invariant in the minimum, and is given by where m 2 is the mass in the scalar potential [12].⁷ Although this quantity is scale dependent, the scale can be nailed down when comparing against lattice [22]. So we take |Φ| 2 /T c as a proxy for the sphaleron transition rate. Let's see how to find W ≡ 2 ∂ ∂ m 2 V (φ, T ) order-by-order. Take the potential as where all resummations are included as in the previous section. And equivalently (3.31) The expansion is We straight off the bat see a different story than for V min : φ NLO contributes already at NLO-this didn't happen for V min because ∂ V LO | φ LO = 0.
We have verified that W is gauge invariant and finite up to order NNLO. The phase transition strength is then and similarly for φ NNLO . On the other hand of the spectrum there's the latent heat, denoted L(φ, T ), and defined as L ≡ T ∂ T ∆V [17]:

Summary of procedure
With the considerations of the previous subsections out of the way, we can now summarize the recipe for calculating the vacuum energy at finite temperature. The first thing to consider is the size of the four-point coupling λ compared to the coupling of other bosons (we use a generic gauge coupling e to facilitate discussions). If λ ∼ e 2 a second-order transition takes place at leading order. The perturbative expansion of T c seems to break down for higher orders. The critical temperature T 2nd can be readily found from the leading-order potential. In our notation it can be written . (3.36) If λ ∼ e 3 , then a first-order phase transition is possible. If this is the case, finding the critical temperature requires performing a perturbative expansion-there are a number of steps. In the high-temperature expansion, a bosonic mode X contributes to the 1-loop effective potential as while a fermionic mode contributes with a similar T 2 and T 0 term but no T term.
1 The leading order potential is where the sum over X ranges over all particles, the sum over X over bosons whose masses scale as e 0 (i.e. not scalars) but that are not resummed, and the sum over X are for e 0 -scaling masses that are resummed (3D-longitudinal modes).
2 To find higher orders (NLO, NNLO, . . . ) use the following scaling rules. Perform the high-temperature expansion for each loop order. The scalar squared masses count as a factor of e, other squared masses as e 0 . Each factor of T, φ (outside masses) count as e −1 .
3 Excluding the T 2 V 2 1 term, all scalar masses should be resummed. The new masses are found from the leading-order potential, e.g. the generated extra terms, However, this requires singling out the terms that should not be resummed, and has to be done while resumming and in a specific order. In practice it is simpler to insert the resummed masses first and then also subtract off all double countings,⁸ (3.40) 5 Find φ min in the perturbative expansion by solving ∂ V = 0 order by order in e. Each φ derivative scales as e. 6 Evaluate V min perturbatively by expanding φ min in the expression for V . 7 Find T c by solving V min = V (0) either by varying T continuously or by performing a perturbative expansion (our recommendation). When expanding T c use the effective VeV to include implicit T derivatives: And T derivatives scale as e 0 when acting on m 2 eff , and as e when acting on anything else. 8 Other observables are found by similarly expanding around T c and φ min .

Examples
This section shows explicit calculations; we'll use Abelian Higgs and the Standard Model as guinea pigs.

The e ective potential
Our aim is to be accurate and have confidence in our accuracy. To test gauge dependence, and by extension our power-counting, we need to calculate the effective potential to 2loop order in a generic gauge. Thankfully, such calculations are tractable since the 2-loop effective potential for a generic model, and gauge, is known at zero temperature [21]. We extend these results to finite temperature and calculate some observables to NNLO accuracy. The calculation is structured so that everything depends solely on a number of master integrals. Unrenormalized master integrals are distinguished from renormalized ones by a boldface. For example, the unrenormalized 1-loop thermal integral, depending on a bosonic mass X , is f(X )-T dependence is left implicit. High-temperature expansions of all master integrals are given in appendix A.
The 1-loop effective potential is a mere functional determinant. And so contributions from different fields separate. Letting X represent bosons, and Y fermions, the generic result is where λ, g, a, c are combinations of couplings and masses, and f SSG is an integral function given in terms of master integrals and masses. At two loops there are are also contributions from counterterm insertions in the 1-loop potential. We opt to perform these insertions explicitly, rather than using the renormalized integral functions of [21].
We review a few integral functions in appendix B; these deviate slightly from [21] with the inclusion of new temperature pieces.
We perform several cross-checks to ensure the validity of our results. Both in Abelian Higgs and the Standard Model. We have tested gauge invariance in the standard loop counting and in the modified ones; renormalization group invariance; and the removal of the "dangerous" 2-loop T 3 φ term [20].

Model de nition
We use the same conventions as [21] to make comparisons easy. The Abelian Higgs model is defined by the Lagrangian The parameters in the scalar potential satisfy m 2 < 0, λ > 0, so that there is spontaneous symmetry breaking at tree-level. Focus on Fermi gauge. The gauge-fixing and ghost terms are Expand the scalar field Φ around its VeV φ as where H and G are real scalar fields. This gives the tree-level potential V 0 (φ) = (1/2)m 2 φ 2 + (1/4)λφ 4 . Using the convenient notation that the squared mass of a particle is denoted with the name of that field, the squared masses-as functions of φ-are H = m 2 + 3λφ 2 , (4.7) G = m 2 + λφ 2 , (4.8) Ghosts are massless in Fermi gauges. The Goldstone field, G, mixes with the gauge boson, Z; propagators of G, Z, and mixed G-Z are expressed in terms of the gauge-dependent masses Z ± , which fulfill (4.11) In addition, Z ± (φ 0 ) = 0, while Z + (0) = G| φ=0 , Z − (0) = 0. These relations ensure that the effective potential is gauge invariant when expanded around the tree-level extrema.

The perturbative expansion
The 1-loop potential is found by summing over all fields, Contributions from ghosts are ignored because ghost masses are φ independent. After renormalization the result is where un-bolded functions are finite. They are given in appendix A. The most important terms in V 1 are where φ independent terms are ignored. As emphasized in [5], the T 2 term in V 1 is gauge invariant.
If the scalar coupling λ scales as λ ∼ e 2 , a second-order phase transition takes place. The leading-order potential for this scaling is The critical temperature is reached when the φ 2 term vanishes, giving On the other hand if λ ∼ e 3 , a first-order phase transition can take place. Finding the leading-order potential (and higher orders) is then more involved; additional terms from V 1 become important: The leading-order potential for this scaling is Masses are of order Z, Z L ∼ e 0 , H ∼ e, Z ± ∼ e close to the minimum because φ ∼ e −1 . A resummation is needed. Both scalar masses are resummed in one sweep via H → The resummed potential describes a first-order phase transition and the critical temperature is known analytically to leading order.
Next, sub-leading corrections. Both V 0 1 and V 2 2 contribute at NLO. Scalar and longitudinal vector masses are resummed as in appendix C. Counterterm insertions also contribute: If a mass X is renormalized in dimensional regularisation by Z X , the finite counter-term contribution is The critical temperature is found in powers of e, T c = e −1 T LO +eT NLO + e 3/2 T NNLO +. . ., as in subsection 3.3, Only scalar T V 1 1 terms contribute at NNLO, and these terms are gauge dependent. But to find V min , V N N L O is evaluated at φ L O -where the Goldstone mass vanishes. So the result is gauge invariant. The corresponding contribution to the critical temperature is Higher orders require 3-loop calculations.

Model de nition
Free parameters are m 2 and λ in the scalar sector; g 3 , g, g (SU c (3), SU L (2), U Y (1)) in the gauge sector; the top-quark Yukawa coupling is y t in the Yukawa sector. We use Fermi gauge with the gauge-fixing parameters ξ γ , ξ Z , ξ W ; the effective potential is independent of the gauge-parameter of the gluon, ξ c , to two loops. In our numerical calculations ξ γ , ξ Z , ξ W ≡ ξ.
The field content for n G generations is Vectors: A, Z, W R , W I , g 8 , where g 8 denotes the color octet gluons.⁹ All bosons are real with the index R denoting the real part and I the imaginary part of the corresponding complex field. The Goldstone G 0 corresponds to the longitudinal mode of Z; G I to that of W R ; and G R to that of W I . Though the real and imaginary part of a field have the same squared mass, the R and I labels are handy for calculations [21].

The perturbative expansion
The barrier height is given by the coefficient of the φ 3 term. In Abelian Higgs this corresponded to the gauge coupling e. In the Standard Model, linear T terms come from W and Z bosons: W 3/2 = (g/2) 3 φ 3 , Z 3/2 = ( g 2 + g 2 /2) 3 φ 3 . Both coefficients are of similar numerical size; we can count the gauge couplings as g, g ∼ e.
The renormalized 1-loop potential is With leading terms For a second-order transition λ ∼ e 2 , the leading order potential is With critical temperature For a first order scaling (λ ∼ e 3 ) the leading order potential is The potential can not be minimized analytically, but the minimum is found numerically in a breeze. Higher-order corrections to the critical temperature are found as in Abelian Higgs.
Though there is a new complication beyond 1-loop: 3D-longitudinal modes of Z and A mix. See appendix C.1 for the details.

Numerical results in the Standard Model
In this section we report on the numerical results from applying our method to the Standard Model, and we compare it with the traditional method of numerically minimizing the potential. The calculation is performed using the high-temperaure expansions given in appendix A, and the organizational framework of [21]; we use Fermi gauge throughout.
We take the input parameters to be [23] When we vary the Higgs mass while keeping the VeV v fixed, we need to vary the potential parameters λ, m 2 in tandem, since

Traditional method
The traditional method finds the critical temperature T c by first calculating the effective potential to a given loop order, minimizing the potential numerically, and changing the temperature till the symmetric and broken energies coincide. We performed this calculation for the Standard Model, using the parameters laid out above. For the 1-loop potential we elected to perform the high-temperature expansion to order O(x 2 T 0 ) instead of using the numerical integrals directly, noting that the expansions are accurate in the temperature range we are considering. 2-loop sunset diagrams are truncated at O(x T 2 ). Scalar and 3D-longitudinal vector masses are all resummed with the leading-order-T 2 self-energy; the resummation targets zero-modes and the scalar O(x 2 T 0 ) contributions. We subtract the relevant diagrams from the 2-loop potential to prevent double-counting.
There are various observables to pick and choose from; some more gauge dependent than others. For example, figure 1a shows that the critical temperature depends weakly on the gauge fixing parameter ξ. This is not unexpected, because two-loop corrections are suppressed with e 2 compared to the leading-order T c -likewise with 1-loop gauge dependent terms. So all gauge dependence is suppressed by e 2 ∼ 1/100. This is in contrast to other observables where the gauge dependence is more prominent. For example, both the barrier height in figure 1b and the transition strength in 1c are quite gauge dependent,even for small ξ. Indeed, while the 1-loop barrier height is relatively insensitive to ξ-the 2-loop barrier height is not. And the same story with the latent heat in figure 1d.

Gauge-invariant method
For the gauge-invariant method we assume that the quartic coupling scales as ∼ e 3 . All calculations are done according to section 3.
It is evident from figure 2a that higher-order corrections are suppressed when determining T c . This is expected. Sub-leading corrections to T c are suppressed by a factor of e 2 ∼ 1/100 to the leading order result.
Higher-order corrections to other observables are more pronounced. There's quite an increase of the barrier height for example. Though, the extra terms in the ħ h expansion are putting in some work. The extra T NLO ∂ T ∆V LO term in (5.15) reduces the barrier height quite a bit-especially at large m H . However, note that radiative corrections to the barrier height are large both in the gauge-invariant and in the traditional method. This doesn't necessarily mean that perturbation theory is unreliable. Only that the leading-order barrier height is small. Indeed, figure 2b shows that the NNLO result is of the same order as NLO. An N 3 LO calculation would be a great cross-check on the convergence.
Finally note that all the results are unreliable for small m H . Because the λ ∼ e 3 scaling is not valid. So perturbation theory breaks down.

Comparison of traditional and gauge-invariant method
Note from figure 3a that T c is quite insensitive to ξ, and that the results of the gaugeinvariant method coincides with that of ξ = 0 (Landau gauge). This is however not the case for other observables. The barrier height is tremendously gauge dependent at T c , as seen in figure 3b. And even the Landau gauge result is an order of magnitude larger than the gauge-invariant result for certain Higgs masses. Which indicates that finite pieces, missed by the traditional method, are significant.
The phase transition strength in figure 3c and latent heat in 3d also showcase a ξ sensitivity. We conclude that all results are quite sensitive to the gauge parameter with the exception of T c -which for Landau gauge coincides with the gauge-invariant method.

Discussion
We showed in this paper how to include gauge-invariant resummations in the finite temperature effective potential. Beyond gauge invariant results, our method includes finite contributions that are missed by contemporary methods. We showed how first-order transitions, which appear highly gauge-dependent, are consistently described in this framework, and we used these methods to calculate a variety of observables-comparing our results to those of the standard method. Part and parcel of the method is the use of a consistent power-counting. Though the specific first-order scaling was first explored in [8], a gauge-invariant method have until now remained elusive. True, some gauge invariant calculations are known [5,17], but these are incomplete or focus on second-order transitions.
Others [24] have put bounds on the gauge dependence. These authors showed that gauge dependence is suppressed for small ξ. A pragmatic approach could then be to take ξ = 0 and ignore gauge dependence all-together. Yet gauge-dependence is but the forerunner of the real issue: an inconsistent power counting. It hardly matters that the result is weakly gauge dependent when there are missing gauge-independent terms of unknown size.
In section 5 we compared two different methods for calculating the critical temperature and the barrier height. The first is the vanilla gauge-dependent method, and the other is our gauge-independent method. We conclude that only some of the observables have small gauge dependence. Not all.
For comparison, the renormalization scale dependence might or might not be larger than the gauge dependence. But these uncertainties are fundamentally of different nature. The dependence on the renormalization scale expresses that our perturbative result is not perfect, and is reduced for each order included. Or by EFT techniques. The gauge dependence shows that terms are missing, and the problem can even worsen when further orders are included. In short, fictitious renormalization dependence is compatible with perturbative results while gauge dependence is definitely not.
There are several available avenues to continue the research in this paper. One is to calculate N 3 LO corrections. These include O(T 3 ) contributions from three loops, and O(T ) from two loops. The effective potential at zero temperature is known to three loops in Landau gauge [25], so the work required would involve translating the various integral functions and master integrals to finite temperature. This calculation would also require extending the high-temperature expansion of the thermal sunset master integral to O(T )-see [26] for a recent calculation.
Even though the step from NLO to NNLO did not induce a significant change in T c , there is reason to expect that the step from NNLO to N 3 LO will be bigger. The terms at NNLO correspond to a half-power of e, and they are sparse; the order N 3 LO corresponds to a full power of e, with sundry diagrams. This contribution could be bigger just from combinatorics. If this is the case, this awkward pattern might continue up the ranks, where every second order in perturbation theory contributes an insignificant amount.
Another avenue: because the overarching goal of this calculation is to study extensions of the Standard Model, it would be interesting to see this method applied to other models. For models with more complicated scalar potentials some care will have to be applied in comparing the sizes of couplings. When many different couplings are involved it might be more difficult to consider the different scaling laws needed to create a barrier. But possible. With scaling laws established, then comes the issue of performing the perturbative expansion. The leading order contribution is straightforwardly calculated once the vector bosons' thermal masses are known, and is in fact easier to calculate than the traditional way (using the full 1-loop potential for all particles).
At next-to-leading order requires a two-loop calculation, which is beyond the norm of phenomenology. Because the two-loop effective potential is known for a general model [21], and since we have extended it to finite temperature in this paper, this calculation can in principle be fully automated.
Finally, thermal resummations can typically be implemented very economically using high-temperature EFT methods. It would be interesting to explore whether the powercounting of this method can be realized in such a high-temperature EFT.

A The thermal master integrals
In this appendix we give the leading terms in the high-temperature expansions of the thermal master integrals, using dimensional regularisation and M S. For readability we suppress the implicit argument T to the functions.
For the 1-loop functions the expansions are known in closed form [7,9]; we just include the orders we need for the calculation. The 2-loop sunset integrals I(x, y, z), I F (x, y, z) are given to order T 2 .
For each integral we show its definition and its ε expansion. We then give the results for the leading orders in ε and T .
In dimensional regularisation it is separated according to in dimensional regularisation we split this integral up according to In the high-temperature expansion, these individual components are given by (A.14) A ε (x) = T 2 12 log A 24 Q 2 16π 2 T 2 + T x 4π log 4x Q 2 − 2 + O(x). (A.15)

A.2.2 Fermionic bubble
The fermionic bubble is A F (x) ≡ T n p 1 p 2 + x + (π(2n + 1)T ) 2 , (A. 16) in dimensional regularisation we split this integral up according to In the high-temperature expansion, these individual components are given by intermingle the modes. The integral functions used in [21] capture the Lorentz structure of the different diagrams; to perform this resummation we will need to project out the contributions of the different modes.
In that vein, consider the propagator of a massive vector boson, where we hid masses and gauge-fixing parameters in the Lorentz invariant functions A and B. We can focus on the projection operators: P µν (k) ≡ g µν − k µ k ν k 2 projects onto the three transverse modes; k µ k ν k 2 projects onto the longitudinal mode. To find the individual contributions from the transverse modes we introduce the 3D-longitudinal projector P µν L and the corresponding transverse projector P µν T , such that Using the properties of these projectors [27], we can derive the contributions from the 3Dlongitudinal modes to the various integral functions. We use the notation that an index V L corresponds to a 3D-longitudinal mode, and an index V T for the transverse modes. The non-zero integral functions that include at least one 3D-longitudinal mode are the following: In the above list we did not include permutations. As an example, there is also the integral function f V L V T (x, y). For each of the integral functions above there is also a corresponding one that only includes the transverse modes. We can find them by using the full result together with the formulas above. where the thermal self-energies are [7] Π W (T ) = 11 6 g 2 T 2 , (C. 19) This mass matrix is diagonalized by an angle θ -depending implicitly on T and φ. The explicit form of θ can be found using the mass matrix above, we neglect to show it here and simply note the limits Massive eigenstates are Z L , A L ; they have squared masses where we mean that Z L has a + sign in front of the square root, and A L has a − sign. These masses behave as expected, The final complication is that resummed 3D-longitudinal modes have different coupling