On the perturbative expansion at high temperature and implications for cosmological phase transitions

We revisit the perturbative expansion at high temperature and investigate its convergence by inspecting the renormalisation scale dependence of the effective potential. Although at zero temperature the renormalisation group improved effective potential is scale independent at one-loop, we show how this breaks down at high temperature, due to the misalignment of loop and coupling expansions. Following this, we show how one can recover renormalisation scale independence at high temperature, and that it requires computations at two-loop order. We demonstrate how this resolves some of the huge theoretical uncertainties in the gravitational wave signal of first-order phase transitions, though uncertainties remain stemming from the computation of the bubble nucleation rate.


Introduction
The discovery of gravitational waves by LIGO [1] suggests a novel possibility to observe the early universe. The comparatively weak interaction between gravity and matter means that gravitational waves will free-stream after their production, carrying with them a fingerprint of whatever process produced them, a fingerprint which is potentially observable today. In particular, a strong first-order phase transition in the early universe would give rise to a stochastic background of gravitational waves (for recent reviews see Refs. [2][3][4]), which may be visible by planned gravitational wave experiments such as LISA [5], DECIGO [6], BBO [7] and Taiji [8]. In fact recent evidence for a possible stochastic background of gravitational waves by the NANOGrav experiment [9] has been interpreted as the gravitational wave signal of a first-order phase transition [10][11][12][13].
On the other hand, the discovery of the Higgs boson [14,15] and high precision probes of its interactions at the LHC and planned future high-energy colliders play a central role in advancing our understanding of particle physics. In particular, determining the symmetrybreaking pattern of the electroweak sector is a major goal. A first-order electroweak phase transition may provide the necessary departure from equilibrium for the generation of the observed matter-antimatter asymmetry [16], yet this possibility requires new particles with masses TeV which are not too weakly coupled to the Higgs [17]. The determination of the pattern of electroweak symmetry breaking is therefore a clear target for near-future collider experiments.
From an observation of a stochastic background of gravitational waves, one can in principle learn about the underlying physics of the production process. Indeed, it has been mooted that this may provide an experimental probe of particle physics beyond the Standard Model (BSM) which may be complementary and competitive with collider searches, see for example Refs. [17,18] and references therein. However, for this to be feasible, it must be possible to make robust quantitative predictions of the gravitational wave spectrum produced by a given particle physics model. For first-order phase transitions, this in turn requires accurate predictions of the thermodynamics of the phase transition.
The majority of recent literature studying the thermodynamics of cosmological firstorder phase transitions -in a wide variety of models -resorts to one-loop computations that use the (daisy-resummed) thermal effective potential. Such studies are aided by the known explicit form of the potential (in terms of background field-dependent mass eigenvalues) in this approximation [16,[19][20][21][22][23][24], so allowing for relatively straightforward applications even to complicated BSM models. This simplicity has led to a wealth of phenomenological studies, for example  which have broadened our understanding of the possible scope of planned gravitational wave experiments to probe particle physics (as well as viability of electroweak baryogenesis). However, as recently emphasised in Refs. [56,57], such one-loop computations suffer from huge theoretical uncertainties due to thermal enhancements of infrared physics, significantly limiting the possibility of making quantitative conclusions. These amount to orders of magnitude uncertainty in the peak amplitude of the gravitational wave spectrum. 1 Unphysical dependence on the renormalisation scale (µ) was noted as one of the largest theoretical uncertainties in typical one-loop computations [56]. This is in stark contrast with typical one-loop computations at zero temperature, in which renormalisation scale dependence is usually small, at least for theories with small couplings. As a consequence, this issue has been the subject of confusion in the literature. Incomplete attempts at renormalisation group (RG) improvement of the thermal effective potential were made for example in Refs. [47,52,[59][60][61][62][63][64][65]. Note that renormalisation scale dependence can be regarded as a proxy revealing the size of missing higher-order terms; a strong renormalisation scale dependence reveals important missing terms.
In Fig. 1 we show an example of the renormalisation scale dependence of the gravitational wave spectrum from a first-order phase transition, in the real singlet-extended Standard Model (xSM). Shown are predictions for one benchmark parameter point based on two different approximations: the widely used one-loop thermal effective potential in (dotted) green, and a more complete computation in blue. In the one-loop calculation, both the peak amplitude of the gravitational wave spectrum and the LISA signal-to-noise ratio vary by four orders of magnitude as the renormalisation scale varies over just a factor of four, a huge theoretical uncertainty. This intrinsic uncertainty is much smaller for the more complete computation The coloured lines and associated theoretical uncertainty bands show the predicted gravitational wave spectrum for a single benchmark parameter point in this model (labelled BM1 in Sec. 5). The uncertainty bands reflect the variation as a result of varying the renormalisation scale over a factor of 4 for these two different approximations. For more details see Sec. 6. The LISA signal-to-noise ratio, assuming a 3 year mission profile and using PTPlot [3] for the computation, varies from 0.017 to 0.19 in the O(g 4 ) approximation, and from 2.3 × 10 −6 to 0.17 in the one-loop approximation.
that is organised in powers of couplings (rather than the loop expansion) and accounts for (equilibrium) contributions consistently up to O(g 4 ). Here g is a formal power counting parameter, that is identified with the SU(2) gauge coupling in the case of the xSM. For a detailed discussion of this analysis, see Sec. 6. As we show in this article, for Z 2 -symmetric theories, the trouble arises from the implicit running of thermal corrections to quadratic terms in the potential 3 leading to the strong renormalisation scale dependence of the high temperature effective potential.
Here v is a real background field, m 2 is the zero temperature mass parameter and Π T is the one-loop thermal contribution to the self-energy, i.e. the one-loop thermal mass. Denoting the interaction coupling by g 2 , the thermal contribution to the self-energy has the schematic form Π T ∼ g 2 T 2 . The implicit running of the Π T v 2 term is an effect at O(g 4 T 2 v 2 ) contributions. 3 For non-Z2-symmetric theories there is a similar troublesome thermal contribution to the tadpole. and is not cancelled by any one-loop, or resummed one-loop, contribution to the effective potential. It is cancelled by explicit logarithms which only appear at two-loop level.
Thus the cancellation of renormalisation scale dependence requires more work at high temperature than it does at zero temperature. This is a specific case of a more general feature of physics at high temperatures: due to the high occupation of infrared modes, the effective expansion parameter is the square root of the coupling, necessitating higher-order calculations. While the one-loop approximation may suffice at zero temperature, one must work harder at high temperatures to achieve the same degree of accuracy.
The realisation that thermal effects can correct the effective potential at leading O(g 2 ) order goes back to Ref.
[67], where the thermal correction in Eq. (1.1) was derived. The more general misalignment of the loop expansion and the coupling expansion at high temperature was further clarified in Ref. [68], in which the presence of a contribution at O(g 3 ) was demonstrated. However, the running of the tree-level potential starts at O(g 4 ), so it was not until calculations were pushed to this order [69,70] that issues of renormalisation scale dependence were directly tackled.
Ref. [70] pioneered the perturbative study of electroweak-type theories at high temperature. The approach was based upon organising perturbation theory as an expansion in powers of couplings, and not by the loop expansion. Strictly adopting such an expansion, the authors established that there is a unique way to handle the infamous ring or daisy diagrams, thereby fixing the O(g 3 ) term in the effective potential. Furthermore, the computation was extended to O(g 4 ). This required a two-loop computation, and the authors demonstrated that their results are -correctly -renormalisation group invariant to this order. In Refs. [71,72] this calculation was extended, allowing for a parametrically larger Higgs self-coupling.
A more detailed discussion of renormalisation scale dependence at high temperature was given in Ref. [73], formulated in the framework of (high temperature) dimensional reduction, a framework that was at least implicit in Ref. [70]. In this, an important distinction was made between the renormalisation scale dependence associated with the heavy non-zero Matsubara modes, and that associated with the light zero Matsubara modes. Dimensional reduction was further developed in Refs. [74][75][76] (for reviews see Refs. [77][78][79][80]). An early triumph of the approach was the determination of the phase diagram of the Standard Model electroweak sector [81][82][83]. Building on these earlier works, later studies using dimensional reduction to study hot electroweak theories have typically performed complete O(g 4 ) computations, for which the cancellation of renormalisation scale dependence serves as a useful crosscheck (c.f. for example [84][85][86][87][88][89][90][91][92][93][94][95][96][97]).
The majority of recent literature on cosmological first-order phase transitions utilises the daisy resummation scheme developed in Ref. [70], but falls short of carrying out a complete calculation at O(g 4 ), in particular missing logarithms of the renormalisation scale from two-loop thermal masses and vacuum diagrams (as well as logarithms related to field renormalisation). Hence, the results of these studies exhibit a strong renormalisation scale dependence, which cannot be removed by RG running of tree-level and one-loop terms, due to the incompleteness of the computation at O(g 4 ).
In this article, we aim to clarify the structure of the perturbative expansion at high temperature, as revealed through the lens of renormalisation scale dependence. We focus on the effective potential, as this carries much of the important thermodynamic information, though our general conclusions are not limited to the effective potential. The remainder of this article is organised as follows. To illustrate the crucial issues without unnecessary complications, we start by working with a simple example: the Z 2 -symmetric real scalar field, or φ 4 -theory. In Sec. 2 we review the RG improvement and scale independence of the effective potential at zero temperature. In Sec. 3 we show how this fails at high temperatures, leaving an uncancelled renormalisation scale dependence for the one-loop thermal effective potential. Further in Sec. 4 we show how this can be resolved by adding the leading two-loop contributions, which are O(g 4 ) in the coupling expansion. In Sec. 5 we turn to a realistic candidate theory of particle physics, the xSM, and demonstrate the numerical importance of our conclusions for predictions of equilibrium thermodynamics. In Sec. 6 we extend this discussion to the gravitational wave spectrum. Finally, we summarise and conclude in Sec. 7. For completeness, Appendix A collects some explicit technical details of the computation.

Scale independence at zero temperature
To demonstrate issues of scale dependence without unnecessary complications, we initially work with the simplest possible theory, that of a single real scalar field, φ, with a Z 2 symmetry, φ → −φ. The Lagrangian density is We regularise the theory using dimensional regularisation, and choose our counterterms according to the MS scheme, see Appendix A . The arbitrary renormalisation scale introduced in this way is denoted by µ, and the parameters m 2 (µ) and g 2 (µ) are the corresponding renormalised MS parameters. We denote the scalar self-coupling by g 2 , in analogy with gauge theories. In the following, we will not always show the explicit argument µ, but one should keep this renormalisation scale dependence in mind. We will assume throughout that the theory is perturbative, so that the loop-expansion parameter is small, g 2 /(4π) 1. At zero temperature, the phase structure of the theory can be determined from the minima of the perturbatively computed effective potential. To calculate this, one shifts the field by a constant, homogeneous background φ → v+φ, and then carries out a loop expansion for the fluctuating field φ, dropping linear terms [98]. At leading order, the result is just the tree-level potential evaluated on the background field, For positive m 2 > 0 this has a minimum at v = 0, whereas for negative m 2 < 0 this has two minima at v 2 = −m 2 /g 2 .
The perturbatively computed potential depends on the renormalisation scale, µ, through the Lagrangian MS parameters. For a generic background field, the leading order scale dependence takes the form 4 where we have used the one-loop beta functions, collected in Eqs. (A.7). At one-loop, the effective potential is corrected by fluctuations of the φ field, yielding the Coleman-Weinberg contribution [99] (see also Appendix A), Note the presence of the renormalisation scale µ explicitly within the potential. Here we use the freedom to shift the potential energy by a constant to fix V (0) = 0 for the full one-loop potential.
The explicit scale dependence of the one-loop piece of the potential is As one can see, this cancels the scale dependence of the parameters in the tree-level potential, shown in Eq. (2.5), The cancellation however only holds at leading order, receiving corrections from the one-loop running of parameters within the one-loop potential, and from two-loop corrections to the running of parameters in the tree-level potential. However these are suppressed relatively by g 2 /(4π), and hence can be neglected in a one-loop analysis.
In summary, at zero temperature the effective potential is independent of the renormalisation scale at one-loop order. In fact this holds order-by-order in , 5 or equivalently in g 2 for this theory. Thus renormalisation scale dependence of the zero temperature effective potential is always a higher order effect that can be neglected. The inclusion of the running of couplings within the potential is often called renormalisation group (RG) improvement, see for example Refs. [100][101][102]. Our analysis has been carried out in the MS renormalisation scheme, but the conclusions hold independently of this because the one-loop (and two-loop) beta functions, as well as the logarithm of the Coleman-Weinberg potential, are independent of the renormalisation scheme [103]. 4 Due to the absence of momentum-dependent divergences at one-loop order in this theory, there is no anomalous dimension; see Eq. (A.13). 5 The reduced Planck's constant, , appears as the loop counting parameter, though it is set to unity in natural units. Figure 2: At high temperature, the one-loop correction to the self-energy, or thermal mass, is of the same order as the tree-level mass. Both terms are of order O(g 2 T 2 ). We use the TikZ-Feynman package [104] to draw the diagrams.

∼
3 Scale dependence at high temperature At high temperature, the infrared modes of bosons become highly occupied. As a consequence their effective coupling becomes larger than at zero temperature. This causes a misalignment of the loop expansion and the expansion in powers of couplings, necessitating resummation of infinite classes of diagrams. For a recent review, see Chapters 3 and 6 of Ref [78].
In this section, we will investigate the renormalisation scale dependence of the resummed one-loop effective potential. As we will show, due to the enhancement of infrared modes, the cancellation of scale dependence which occurs at zero temperature no longer occurs at high temperature. Thus, at high-temperature, the resummed one-loop effective potential is strongly scale dependent.
At high-temperature thermal fluctuations become important, so that the tree-level potential is no longer the leading order approximation to the effective potential. As has long been known, near the critical temperature of a phase transition, the dominant thermal corrections are to the mass, and these compete with the tree-level mass [67]. Here and throughout we adopt the high-temperature approximation for thermal functions. For the Z 2 -symmetric scalar theory, the leading thermal contribution to the self-energy takes the form Near a phase transition it must be that m 2 and Π T are the same order of magnitude, so that This is illustrated in Fig 2. At such high temperatures, the loop expansion is no longer appropriate. Instead, the effective potential is computed as a power series in g 2 . For strong symmetry-breaking transitions, strong enough to provide the necessary departure from equilibrium for successful electroweak baryogenesis, there is the further parametric relation v ∼ T . For ease of counting powers, we will often use this parametric relation when discussing the magnitude of different terms in the effective potential. However, we will not assume this relation in our analysis. 6 The thermal effective potential at leading order in g 2 is instead, It is thus of order O(g 2 T 4 ) for v ∼ T . Note that v ∼ T also follows from demanding that the different terms of the effective potential are of the same order, Π T v 2 ∼ g 2 v 4 , as one would expect in the vicinity of a phase transition. In contrast to the effective potential, the renormalisation group equations are unaffected by finite temperature [107], being only dependent on the deep ultraviolet (UV). However, because the effective potential is modified at leading order, so too is its running, At subleading orders, half-integer powers of the coupling constant arise [68]. This feature is a further consequence of the high occupation of infrared modes, unique to perturbation theory at high temperature. The one-loop contribution of the zero Matsubara mode, appropriately daisy-resummed [70], gives the first correction to the leading thermal effective potential, suppressed relatively by g 2 . This amounts to, where again we add a (temperature-dependent) constant to fix V (0) = 0. We do this because we are interested in differences between the two phases, and not in the absolute value of the free energy or pressure. Eq. (3.6) is a correction to the effective potential of order O(g 3 T 4 ) for v ∼ T . Perhaps the most common approximation taken in the literature is a resummed one-loop approximation, in which thermal corrections to the effective potential are incorporated to one-loop order together with a resummation of daisy diagrams [30-34, 39, 45, 47, 50, 53, 54]. 7 For example, this approach has been adopted in numerical packages for analysing finitetemperature cosmological phase transitions [22][23][24]108]. In this approximation -and in addition with the high-T expansion that we assume -the effective potential takes the form

7)
refer to first-order phase transitions as these are of more interest to us. Our general conclusions regarding renormalisation scale dependence however will not depend on the details of the transition, and as we show later in this section, neither will they depend on the model. 7 Note that all these references use Arnold-Espinosa type daisy-resummation [70], where only the mass of the zero Matsubara mode is resummed, contrary to the approach by Parwani [69] in which all modes are resummed.
where, following Refs. [73,74], we have defined (γ is the Euler-Macheroni constant) Definitions of the various different terms in first line of Eq. (3.7) are standard, and are given in full in Appendix A, together with a derivation of the second line. The second line utilises a split between soft and hard terms, originating respectively from the zero and nonzero Matsubara modes. In this expression powers of g 2 are separated but zero temperature and high temperature pieces are mixed. The last term arises at one-loop, but is of order O(g 4 v 4 ), just as the Coleman-Weinberg potential at zero temperature. However, due to the enhancements of infrared physics, this is not the full O(g 4 ) piece of the thermal potential, but only a part thereof, as we will see in next section.
In analogy with the result at zero temperature, it might be guessed that the renormalisation group running of this one-loop approximation, Eq. (3.7), is subdominant. However, as we will show, this is not correct. Explicitly for this scalar theory, the leading order running is where we have dropped irrelevant constant terms. The scale dependence of the first two terms cancels exactly as in the zero temperature case since However, at the same (leading) order there is the following leftover uncancelled term, which is (dropping an irrelevant constant) This is the troublesome term foreshadowed in Eq. (1.1) in Sec. 1. The scale dependence does not cancel at leading order, in contrast to Eq. (2.9) at zero temperature. As a consequence, at high temperature the running of couplings is not subdominant in this approximation, being of the same order as the Coleman-Weinberg term itself. That is, at high temperatures the one-loop effective potential is scale dependent at order O(g 4 ), whereas at zero temperature the scale dependence of the one-loop potential was only of order O(g 6 ). In fact, this scale dependence is the lowest order that scale dependence could arise, being the same order as the running of the tree-level potential.
This equation, also holds to leading order in more complicated theories. In theories with non-Z 2 -symmetric scalars, there is also a term at this order arising from thermal contributions to the tadpole; see for example Refs. [80,96,97]. For a beyond the Standard Model theory with an extra Z 2 -symmetric scalar S the thermal self-energy of the Higgs scalar is schematically where g, g are SU(2) and U(1) gauge couplings, g Y is top Yukawa coupling (effect from other Yukawa couplings is numerically negligible), λ h the Higgs self-interaction and λ p the portal coupling between Higgs and scalar S. The numerical factor C depends on the representation of S under the gauge symmetries, for e.g. C = 2/3 in the real-singlet extended SM [80], C = 2 in the real-triplet extended SM [93] and C = 4/3 in the Two-Higgs Doublet Model (2HDM) [92]. 8 In the pure SM case, the leftover scale dependency is numerically dominated by the top quark since g 2 Y > g 2 , g 2 , λ h . The situation is even worse in those BSM theories where large portal couplings λ p /(4π) 1 are required to change the character of the electroweak phase transition from a crossover to a strong first order transition. The running of these large portal couplings prevents an accurate determination of the critical temperature (for example see Fig. 12 in Ref. [94]), as well as all other thermodynamic parameters.
For the case of radiatively-induced first-order phase transitions, our conclusions regarding the magnitude of the renormalisation scale dependence are in fact an underestimate. For the SM with light Higgs, this happens for λ h ∼ g 3 , (see for example Ref. [70], or for a recent discussion Ref. [109]), and in the vicinity of the transition the O(g 2 ) terms in the effective potential partially cancel, leaving an O(g 3 ) remainder. As a result the renormalisation scale dependence of Eq. (3.6) is of the same order as that of the tree-level potential, yielding an additional uncancelled scale dependence at leading order.

Scale independence at high temperature
The one-loop approximation to the thermal effective potential given in Eq. (3.7) is incomplete at O(g 4 ). This is the reason for the residual renormalisation scale dependence at this order. In the following, we will show explicitly how the scale dependence cancels in a complete calculation at O(g 4 ). This requires the computation of two-loop Feynman diagrams, including two-loop corrections to the thermal mass, as illustrated schematically in   Figure 3: At high temperature, the running of tree-level parameters is the same order as the running of the one-loop thermal mass, which in turn is the same order as explicit logarithms of renormalisation scale at two-loop order. All these terms are of order O(g 4 T 2 ).
Perhaps the most widely used way to derive the full O(g 4 ) thermal effective potential of a given theory is to utilise high-temperature dimensional reduction to a three-dimensional effective field theory (3d EFT); see Refs. [73][74][75]110] for the original literature and Refs. [56,77,78,80,96] for reviews. The technique of dimensional reduction is a systematic approach to the resummations required at high temperature, order-by-order in powers of couplings.
The O(g 4 ) result for the Z 2 -symmetric real scalar theory -that diagrammatically requires a two-loop determination -has been derived long ago [70,73,75,110]. Expanded to this order, the result reads where, following Refs. [73,74], we have introduced the notation where A is the Glaisher-Kinkelin constant and again we may use the freedom to add a constant to the potential to set V thermal (0) = 0. The full result for the two-loop effective potential within the 3d EFT can be found in Eq. (A.26) in the appendix. 9 This result differs from the one-loop thermal effective potential, Eq. (3.7), by the terms on the first line of Eq. (4.1), in square brackets. From the perspective of high-temperature dimensional reduction, these terms arise from two sources: from the two-loop correction to the thermal mass, as well as from two-loop vacuum diagrams within the 3d EFT. Despite arising from two-loop diagrams, the terms are clearly of the same size as the one-loop Coleman-Weinberg terms for v ∼ T , both in terms of powers of g 2 and in terms of powers of 1/(4π).
This result for the O(g 4 ) part of the effective potential depends explicitly on the renormalisation scale µ, through the terms L b (µ), as well as implicitly through the MS parameters, though this latter dependence is of higher order. The leading order renormalisation scale dependence is This exactly cancels the implicit renormalisation scale dependence of the leading order potential, Eq. (3.9), so that The cancellation however receives corrections which are of higher order. The leading corrections are O(g 5 ) due to the (implicit) renormalisation scale dependence of the O(g 3 ) term.
Cancellations analogous to that of Eq. (4.4) happen generically, also in more complicated theories, as indeed they must. This has been explicitly verified in the case of the SM in Refs. [70,[72][73][74]. Similar computations in theories beyond the SM include [92,94] (2HDM), [93,95] (real-triplet) and [80,97] (real-singlet); works that all utilise dimensional reduction to three-dimensional effective theories at high temperature. The cancellation of scale dependence is an important consistency check in the construction of these high temperature effective field theories.
In the computation of Eq. (4.1) both UV and IR logarithms arise, though in Eq. (4.1) the renormalisation scale dependence of the IR logarithms has already cancelled. The UV logarithms are related to the usual renormalisation group running of couplings. The latter, the IR logarithms, are instead related to the renormalisation group running within the 3d EFT [73], as is characteristic of EFTs [113]. Making use of this, the renormalisation scale of the IR logarithms can be replaced with a new scale, µ 3 . The scale µ 3 can then be set independently of µ, thereby reducing the occurrence of large logarithms and improving the convergence of perturbation theory. It also has the beneficial effect of separating out the effects of IR and UV scales, which we will make use of in the following.

The consequences for equilibrium thermodynamics
The residual renormalisation scale dependence at O(g 4 ) of the common one-loop approximation to the thermal effective potential, Eq. (3.7), implies that there are corresponding theoretical uncertainties for all physical quantities computed in this approximation. In this section, we focus on quantities which can be computed solely based on knowledge of the thermal effective potential. These are necessarily bulk equilibrium quantities, and we focus on the critical temperature T c and the phase transition strength α evaluated at T c . This latter quantity is proportional to the latent heat when evaluated at T c ; see Ref. [56] for a definition. We postpone to Sec. 6 discussion of the bubble nucleation rate, which is also necessary for determination of the gravitational wave spectrum, but which does not involve the effective potential in its computation.
As an explicit test of the magnitude of the residual renormalisation scale dependence, we investigate the simplest extension of the Standard Model with a first-order phase transition, namely, the Standard Model plus a real singlet scalar field (xSM). The computations require generalising the results of earlier sections to include the Standard Model field content. The details for this can be found in Refs. [80,97,114], where the high-temperature 3d EFT of this model was constructed and the effective potential computed to two-loop order. The calculation of equilibrium quantities is thus complete at O(g 4 ). 10 In computing T c and α c ≡ α(T c ) we follow the 3d approach outlined in Ref. [56]. That is, we perform dimensional reduction and then compute the relevant thermodynamic quantities in an -expansion within the 3d EFT, thereby ensuring order-by-order gauge invariance. To circumvent the difficulties of phase transitions which are radiatively induced at the soft (or ultrasoft) scale, we focus on two benchmark points with a two-step transition. The first step is a transition in the singlet direction, followed by a transition to the usual Higgs minimum. For the second step of the transition, there is a tree-level barrier between phases, leading to a stronger transition. We study only this second step. The coupling expansion is more slowly convergent for radiatively induced one-step transitions, being essentially an expansion in √ g [70,109], rather than in g as for the tree-level driven two-step transitions we focus on.
A parameter point in the Z 2 -symmetric xSM can be specified by the values of {M σ , λ m , λ σ } where M σ is the pole mass of the singlet scalar, λ m is the Higgs-singlet portal coupling and λ σ is the singlet self-coupling [80], the latter two are denoted as a 2 and b 4 respectively in Ref. [97]. We choose the following two illustrative benchmark points, BM1 and BM2  10 In fact, there are some contributions still missing at O(g 4 ), also left out in Ref. [74], which come from the temporal components of gauge fields in the second step of dimensional reduction and which we expect to be numerically small. One way to include these terms would be to perform the phase transition analysis at the soft scale, rather than integrating these fields out. Note that the loop expansion of the effective potential within the EFT is an expansion in powers of ∼ g, starting with the tree-level potential at O(g 2 ). In the above the order at which dimensional reduction has been carried out is matched by the order of the loopexpansion within the EFT.
In (i)-(iv) above, we have separated out the orders at which the (UV) dimensional reduction was carried out, from those at which the (IR) calculation within the 3d EFT was carried out. The former depends on the nonzero Matsubara modes, and the latter on the zero Matsubara modes. We have done this to highlight how the calculation factorises into these two parts, and so too does the renormalisation scale dependence. While the dimensional reduction step depends on the original renormalisation scale µ, the calculation within the 3d EFT depends on its own renormalisation scale µ 3 [73].
The one-loop approach (ii) is approximately equivalent to the usual one based on the one-loop daisy-resummed potential of Eq. (3.7) [56]. However, approach (ii) also includes the effects of one-loop wavefunction renormalisation, which enter through the dimensional reduction step, thus it includes all logarithms of scale which are present at zero temperature, including the effects of nonzero anomalous dimension; see Eq. (A.13). As a consequence, the magnitude of the renormalisation scale dependence of approach (ii) gives a conservative estimate for one-loop approaches in the literature, which typically neglect this effect. Further, by utilising dimensional reduction, we are able to maintain order-by-order gauge invariance and also to avoid double counting problems and an uncontrolled derivative expansion in the bubble nucleation calculation in Sec. 6.
In Fig. 4 we plot T c and α c as functions of the renormalisation scale µ in all four different approximations. In addition, for the O(g 4 ) approximation we show dependence on the 3d EFT renormalisation scale µ 3 . Note that in the 3d EFT running in terms of µ 3 starts only at two-loop order, and hence we show it only for the O(g 4 ) approximation. On the x-axis, we use the dimensionless ratio µ 3 /g 2 3 , where g 2 3 is the (dimensionfull) SU(2) gauge coupling in the 3d theory. Fig. 4 shows that only in the full O(g 4 ) approximation is the renormalisation scale under control. All other approximations show a striking dependence on renormalisation scale, amounting to a few tens of percent for the critical temperature, and a multiplicative factor of O(10) for α c . This is especially marked for BM2, for which the couplings are larger and the transition stronger. Comparing the one-loop approximation to the O(g 2 ) and O(g 3 ) approximations, one can see that the renormalisation scale dependence is somewhat smaller in the one-loop approximation but that it is nevertheless of the same order of magnitude, a numerical consequence of it being the same parametric order. In the full O(g 4 ) approximation, the renormalisation scale dependence is at least an order of magnitude smaller, demonstrating numerically that its parametric order has indeed been reduced. Dependence on µ 3 is stronger than dependence on µ, albeit not significantly. We conclude that for the equilibrium properties of the transition, the perturbative expansion is quantitatively under control only when all terms at O(g 4 ) are included, and this is not achieved at one-loop. This conclusion follows naturally in light of earlier sections, as O(g 4 ) is the minimal order calculation that admits RG-improvement.

The consequences for gravitational wave predictions
The consequences of an incomplete perturbative computation at O(g 4 ) are particularly severe for the gravitational wave spectrum produced by a first-order phase transition. This was demonstrated in Ref. [56], where it was found that missing O(g 4 ) terms, revealed through renormalisation scale dependence, give potentially the largest source of error for one-loop computations. 11 The amplitude of the gravitational wave spectrum depends very sensitively on the thermodynamic parameters, and these in turn have relatively large uncertainties. Although the theoretical arguments in Ref. [56] were general, the numerical study was restricted to a simplified version of the Standard Model Effective Field Theory. The arguments of this article at hand pinpoint where the renormalisation scale dependence arises. From this we expect the conclusions of Ref. [56] in this respect to hold rather generally.
In the previous section, we have presented four approximations with which we have computed the purely equilibrium quantities T c and α c . The gravitational wave spectrum depends additionally on the bubble nucleation rate, and the speed of the bubble wall growth, both of which are inherently real-time quantities, depending on the evolution of inhomogeneous field configurations. As a consequence, they are significantly more challenging to compute. In this article we have focused on equilibrium physics, and we refrain from overstepping this remit. Fortunately, at leading exponential order the calculation of the bubble nucleation rate reduces to a purely equilibrium calculation [121,122], O( 0 ) within the 3d EFT. For this we follow the approach outlined in Ref. [56], though approximating the nucleation prefactor simply by T 4 . At O( 1 ) real-time physics enters the bubble nucleation rate (through the dynamical prefactor), as it does for the bubble wall speed at leading order.
As a consequence, we are not able to carry out a full O(g 4 ) calculation of the gravitational wave spectrum. We can nevertheless compare different approximations to the dimensional reduction, though in each case the bubble nucleation calculation is carried out at O( 0 ) within the 3d EFT. As far as we are aware, there does not exist a complete calculation of the thermal bubble nucleation rate at O( 1 ) for any quantum field theory, yet it is necessary to extend to O( 2 ) in order to match the accuracy of the equilibrium calculations. For the bubble wall speed we simply take v w = 1 in the gravitational wave amplitude, noting that there remains considerable debate in the literature as to the leading effects which contribute to this quantity; see for example Refs. [123][124][125][126][127][128]. In the following, for simplicity, we will continue to use the same names for the approximations as in Sec. 5, though one should bear in mind these additional limitations with regard to bubble nucleation and the bubble wall speed.
The issue of renormalisation scale dependence is therefore more complicated for our calculations of the non-equilibrium quantities related to bubble nucleation. As before, the calculation factorises into a UV part (dimensional reduction), and an IR part (within the 3d EFT), the latter only carried out at O( 0 ). By carrying out the dimensional reduction at O(g 4 ), the µ-dependence cancels, as it did for T c and α c . Again this leads to significant improvements over the lower-order calculations. However, the µ 3 -dependence arising from the IR does not cancel, as this would require an O( 2 ) calculation of the bubble nucleation rate. This problematic µ 3 -dependence of the nucleation rate infects the GW signal, even if the equilibrium analysis is performed at O(g 4 ). 11 See also Ref. [115] where the magnitude of various uncertainties [116][117][118][119][120] due to macroscopic, rather than quantum field theoretic, physics were investigated. While certainly important, none of these are quite as severe as those we find here. Our approach to the macroscopic physics can be described as moderately diligent in the language of Ref. [115].  Fig. 4, dependence on µ 3 is significantly stronger due to the limited O( 0 ) accuracy of the bubble nucleation calculation. This is especially true for BM2, for which the scalar couplings are larger.
In Figs. 5 and 6 we present our results at BM1 and BM2 for the thermodynamic parameters that go into the computation of the gravitational wave spectrum; see Ref. [3] for definitions. In Fig. 5 we plot the percolation temperature (T * ) and the phase transition strength α * ≡ α(T * ), as functions of the RG scales. For these, we show all four of the aforementioned approximations, similarly to Fig. 4. In the O(g 3 ) approximation, the full range of RG-scale cannot be shown, since for some values the second step of the phase transition disappears altogether. Just as in Fig. 4, the O(g 2 ), O(g 3 ) and one-loop approximations show a strong RG-scale dependence, of the same order of magnitude for all three. The approximation using O(g 4 ) dimensional reduction shows a markedly weaker renormalisation scale dependence than the other three approximations, especially regarding µ. However, for T * and α * the µ 3 -dependence is much worse than it was for T c and α c , especially at BM2 which has larger scalar couplings. As discussed, this is due to the limited O( 0 ) accuracy of the nucleation calculation. In Fig. 6, we show the inverse duration of the transition (β/H * ) as a function of the RG scales. In addition, we replot β/H * against α * in Fig. 5, emulating the plots produced in PTPlot [3]. This revisualisation makes especially apparent the magnitude of the intrinsic uncertainty, for single benchmark points. For all three of the approximations with less than the full O(g 4 ) dimensional reduction, the intrinsic uncertainty is strikingly large. In particular for BM2, even in the best of our approximations, the variance in the (α * , β/H * )-plane is still very large.
Finally, using these thermodynamic parameters we plot the gravitational wave spectrum due to sound waves [3,129]. The gravitational wave signal of BM1 was already shown in Fig. 1   similar plot for BM2 is given in Fig. 7. For simplicity, in these figures we only show the O(g 4 ) and one-loop approximations; the other approximations have larger uncertainties, shown in Table 1. At BM1 the uncertainty band of the one-loop approximation spans over four orders of magnitude for the peak amplitude, making the prediction quite ambiguous. In the O(g 4 ) approximation, the uncertainty is significantly smaller, but is still one order of magnitude for the peak amplitude. We note that upgrading the one-loop to the O(g 4 ) approximation leads to an increase in the peak amplitude, and shifts the location of the peak to smaller frequencies -though these are expected to be model, and parameter point, specific trends. However, based on our argumentation of the previous sections, the reduction of RG-scale dependence can be expected to hold generically. For BM2, at which there are larger scalar couplings, the theoretical uncertainty in the gravitational wave signal is much larger, and as a prediction for experiments such as LISA it is very ambiguous. In the one-loop approximation the peak amplitude varies by six orders of magnitude, and even in the O(g 4 ) approximation it varies by four orders of magnitude, chiefly a consequence of the uncertainty in the bubble nucleation rate. From this we can conclude that, in order to pursue quantitatively accurate predictions, future work on the bubble nucleation rate is necessary.

Summary
The main conclusions of this article can be restated rather simply. As is well known, renormalisation scale invariance implies that, in the coupling expansion of the effective potential, lower order terms are linked by the implicit running of couplings to explicit logarithms of the renormalisation scale in higher order terms. At zero temperature, the coupling expansion of the effective potential takes the form This schematic formula presents a formal coupling expansion; in theories with multiple couplings there are expansions for each coupling, though in practice one typically establishes formal power counting rules for bookkeeping. At zero temperature the coupling expansion and the loop expansion coincide; A 2 is the tree-level term, A 4 arises purely at one-loop, and in general n-loop diagrams give the A 2(n+1) g 2(n+1) term. The running of couplings takes the form where B 4 is the one-loop term, B 6 is the two-loop term and so on. The one-loop running of A 2 g 2 is an O(g 4 ) effect and is cancelled exactly by explicit logarithms in A 4 g 4 , meaning that a one-loop calculation is sufficient to achieve renormalisation scale invariance at O(g 4 ), receiving corrections at O(g 6 ). Explicitly, the cancellation takes the form At high temperature, the enhancement of IR bosonic modes modifies the coupling expansion of the effective potential, V high-T eff = a 2 g 2 + a 3 g 3 + a 4 g 4 + a 5 g 5 + . . . . (7.4) Odd powers of g arise, and loop orders are mixed in the expansion coefficients. The coefficient a 2 receives contributions at both tree-level and one-loop, the coefficient a 3 comes from the (daisy-resummed) one-loop term, and the coefficient a 4 receives contributions at both (resummed) one-and two-loop. However, thermal effects do not modify the runnings of couplings, as this depends purely on the UV. So, one-loop running links the coefficients a i with a i+2 (and higher loop running with a i+4 and so forth), The leading order running starts at O(g 4 ), and hence this is also the minimal accuracy for any cancellation of the renormalisation scale. As a 4 receives contributions at (resummed) twoloop, no (resummed) one-loop calculation can achieve renormalisation group improvement at high temperatures. Further, to achieve a residual uncertainty of O(g 6 ), equivalent to a one-loop calculation at zero temperature, requires also computing a 5 , which receives IR contributions at three-loop order.

Outlook
At two benchmark points in the xSM, we have demonstrated the numerical importance of these theoretical considerations for the determination of properties of the phase transition.
We have found that calculations at less than O(g 4 ) accuracy show a very large renormalisation scale dependence, in agreement with previous studies in other models [56][57][58]94]. We also remarked that while we were able to perform a complete O(g 4 ) calculation for T c and α c , this is currently out of reach for the bubble nucleation rate and bubble wall speed, imparting a limiting source of uncertainty for the gravitational wave spectrum. Table 1 presents a summary of the magnitude of the scale dependence in the four different approximations we adopted. Though our numerical analysis was carried out merely at two benchmark points, the underlying arguments of Eqs. (7.1) and (7.4) are generic features of perturbative quantum field theory at zero and high temperature. Thus, we expect the qualitative conclusions to apply to the wide variety of models of cosmological phase transitions considered in the literature.
We thus infer that missing O(g 4 ) corrections to the thermal effective potential (and other thermodynamic quantities) are a crucially important, and likely a limiting theoretical uncertainty in predictions of the gravitational wave signal of cosmological phase transitions. Their inclusion is therefore a necessary hurdle to surmount, if one wants to make quantitatively reliable predictions. This applies also to electroweak baryogenesis, and other phenomena which depend on the dynamics of a cosmological first-order phase transition. We advocate that the magnitude of the theoretical uncertainty should be acknowledged and assessed in future analyses of cosmological phase transitions.
We remark that while our O(g 4 ) calculation shows significant improvements over the lower-order calculations, the theoretical uncertainty in the GW peak amplitude is still rather large: one order of magnitude at BM1 and four orders of magnitude at BM2. This motivates future work on the subject, with the aim of improving predictions further. In particular, we have identified the bubble nucleation rate as likely a limiting source of theoretical uncertainty, once O(g 4 ) corrections are included for the equilibrium quantities. This can be inferred from a comparison of Figs Table 1: Summary of theoretical uncertainties due to residual renormalisation scale dependence for the four different approximations. In each case ∆X/X ≡ (max(X) − min(X))/min(X), where the set X consists of the values obtained as the renormalisation scale varies by a factor of 4.
Here Ω refers to the gravitational wave peak amplitude and SNR to the LISA signal-to-noise ratio for a 3 year mission, computed using PTPlot [3].
corrections would be a first, and a significant advancement, whereas one must go to O( 2 ) in order to achieve parity with the equilibrium computations. 13 Nevertheless, higher order computations of equilibrium quantities would allow one to really test the convergence of perturbation theory, and the O(g 5 ) correction in particular may be sizeable, as was found for QCD [130,131]. The O(g 5 ) correction was also found to significantly improve agreement between lattice and perturbation theory for relatively weak phase transitions in the real scalar theory [96]. Computations of equilibrium thermodynamics at O(g 5 ) [89,91,96,130,131] and at even higher orders [132][133][134] have been performed for several theories. The underlying methodology for such higher-order calculations applies also to BSM electroweak phase transitions, though it nevertheless presents challenges.
It is important to emphasise that the discussion in this article is in the context of purely perturbative studies of thermal phase transitions. In fact, due to Linde's Infrared Problem [135] infinitely many (resummed) loop orders contribute to the effective potential at O(g 6 ), rendering perturbative approaches inherently incomplete. Perhaps the only tool which can overcome this problem is lattice Monte-Carlo simulations [76,81,[94][95][96][136][137][138]. Nevertheless, perturbation theory is a valuable guide, and one we can test our confidence in by comparison to the results of lattice simulations.
It may be argued, contrary to our conclusions, that in the context of BSM theories with unknown input parameters, precise calculations are unnecessary as the difference between different approximations may well be accommodated by a shift in the input parameters. To this perspective we present the following counterarguments: First, several orders of magnitude intrinsic uncertainty (as well as uncertainty regarding the order of the transition) is unsatisfactory on purely theoretical grounds, and naturally leads one to doubt the reliability of the calculations, and to aim to improve them where possible. Second, even in full parameter scans, the range of possible gravitational wave signals which can be produced by a given model depends relatively sensitively on e.g. the functional form of the thermal effective potential [37], and hence on the perturbative treatment. Third, by the time the second generation of gravitational wave detectors, such as LISA and Taiji, are due for launch, the LHC Runs 3, 4 and part of Run 5 will have been completed, as well as many other experiments searching for BSM physics. Thus, in the future scenario whereby particle physics experiments are able to point towards some specific type of BSM theory, or -being even more optimistic -even a narrow region of its parameter space, we want to be in a position to compute the thermodynamics relevant for GW predictions as accurately as possible, and for this it is important to acknowledge that accuracy less than O(g 4 ) is potentially unreliable. Finally, reversing this argument, if a stochastic GW background of primordial origin is observed, it is crucial to be able to make as accurate predictions as possible, for there to be any hope to reverse-engineer the underlying particle physics model -i.e. the LISA inverse problem [139].
where renormalised parameters are denoted without subscripts. Here, the wave function renormalisation Z = 1 and δZ = 0 as there are no divergent topologies with external momentum dependence at one-loop order. For the other counterterms, we define at one-loop order We ignore the vacuum counterterm required to cancel the field independent divergence of the effective potential. Bare parameters are independent of the renormalisation scale, and hence we can solve for the scale dependence or running of the renormalised parameters 14 (A. 6) In the limit → 0 we obtain the β-functions Shifting the field φ → φ + v by the homogeneous background v, produces the tree-level potential in Eq. (2.3), and the Lagrangian for the quantum field φ in the shifted theory reads The coefficient of the linear term vanishes at the tree-level minimum. The background fielddependent mass parameter M 2 (v) of the quantum field φ in the shifted theory can be used to derive one-loop quantum corrections to the effective potential, the Coleman-Weinberg potential [99]. The effective potential is the generator of all n-point 1-particle-irreducible correlation functions φ n with zero external momentum, To compute this, we employ a trick [141]: we shift v → v − ω and take a derivative with respect to ω at ω = v. This equals a single tadpole diagram in the theory with background v − ω (quantities of this theory are denoted by tilde) , (A.10) 14 Note that the scale dependence of bare parameters is required to vanish order-by-order in and g 2 , and hence β(g 2 ) starts at order g 2 .
where M 2 (v) and V φ 3 are those in Eq. (A.8) with v → v − ω replaced. Integration over momentum is performed in dimensional regularisation in D = 4 − 2 dimensions in the Modified Minimal Subtraction (MS) scheme [103]. Since V v 3 = d dω M 2 (v), an integration over omega retrieves (up to a constant) where Γ is Euler's gamma function and γ is the Euler-Mascheroni constant. Integrals of this type are computed by employing standard parametrisation tricksà la Feynman and Schwinger [103]. Divergent 1/ poles are cancelled by the tree-level counterterm contribution and the remaining finite part yields V CW (v) of Eq. (2.6). An alternative, and perhaps easier derivation of the β-functions follows from the Renormalisation Group equation -or Callan-Symanzik equation -that requires that the effective potential is scale independent dV eff /dµ = 0. By applying the chain rule, this leads to where we equate terms at O( ). The anomalous dimension γ is given by γ = 1 2 (Π ) −1 µ ∂ ∂µ Π , where Π ≡ d dk 2 Π 2 , i.e. the part quadratic in external momentum k of the self-energy function Π 2 . At one-loop order the self-energy in this theory is momentum independent (the leading momentum dependence arises only at two-loop from the sunset topology diagram) and hence γ vanishes. Then, requiring that the LHS and RHS are equal for any v allows one to solve for the β-functions, resulting in Eqs. (A.7).
Thermal corrections At high temperature, the one-loop contribution is given by (from now on we denote M 2 (v) → M 2 , keeping background field dependence implicit in our notation) where D-dimensional integration at zero temperature has been replaced by a sum-integral, utilising the Matsubara or imaginary time formalism of thermal quantum field theory [78,107]: where d = D − 1 and Euclidean four-momentum is defined as P ≡ (ω n , p) with the bosonic Matsubara frequency ω n = 2πnT , where n is integer. In addition, we have divided 15 the sum-integral into a temperature-independent Coleman-Weinberg piece and a thermal piece 16) where z ≡ M/T , the Bose-distribution is n B (E p , T ) = 1/(e Ep/T − 1) with E p = p 2 + M 2 and a b ≡ (4π) 2 Exp( 3 2 − 2γ). In Eq. (A. 16) we have used the high-T expansion (M/T 1) and dropped a constant term, independent of z. Note that V T is UV-finite, and does not contain logarithms of renormalisation scale.
Merely evaluating Eq. (A.16) is not sufficient to correctly capture the O(g 3 ) contribution to the thermal effective potential, as an infinite set of (daisy) diagrams contribute at this order. A minimal prescription for the necessary resummation [70] is to add the correction V T → V T + V daisy , where the daisy term -associated to the cubic term that is non-analytic in mass squared -reads where M 2 + Π T is the resummed background field-dependent mass parameter. This form is a result of the resummation of the mass of the zero Matsubara mode alone (by adding and subtracting the one-loop thermal mass correction to reorganise the perturbative expansion). Note that the second term in Eq. (A.17) just removes the corresponding cubic term with unresummed mass in Eq. (A.16). In total, the one-loop thermal effective potential then reads Note in particular, that in the high-T expansion, the ln[M 2 (v)] cancels between V CW + V T , leaving only temperature and renormalisation scale remaining inside logarithms. The thermal part V T can also be computed without a high-T expansion, by numerically evaluating the integral in Eq. (A.16). While the one-loop thermal part V T is explicitly free of logarithms of renormalisation scale, at leading order in the high-T expansion it matches the tree-level part in power counting, and hence the running of parameters inside V T is sizeable. We find it helpful to write the thermal effective potential in an alternative form by separating the soft (zero mode) and hard (non-zero mode) contributions 15 For a specific derivation, see [19,67], and for a more generic formula to separate vacuum and thermal contributions -by replacing a thermal sum by a complex contour integral including the Bose distributionsee Sec. 2.2 in Ref. [78] or Sec. 3.4 in Ref. [107]. and explicitly which completes the derivation of Eq. (3.7). Note that V CT cancels the divergent part of V hard . The particular form of the thermal effective potential in Eq. (A.23) helps to illuminate the connection to dimensionally reduced 3d EFT [74,75]. At one-loop order, this connection takes the form (A.25) The full O(g 4 ) effective potential for this theory, given in Eq. (4.1), can be constructed using high-temperature dimensional reduction [73,75,80,96,110]. For this, one needs the two-loop effective potential within the 3d EFT, where for convenience we have defined The residual µ-dependence is O(g 6 ). In addition to these UV logarithms, there arise IR logarithms within the full theory. These match against corresponding logarithms determining the running of parameters within the 3d EFT, allowing one to replace µ → µ 3 within all IR logarithms [73]; see Eq. (A.29). This is characteristic of EFTs more generally [113]. Finally, note that within Eq. (A.26) for the effective potential the µ 3 -dependence of the tree-level mass parameter cancels against that of the explicit two-loop term, leaving a residual µ 3 -dependence at O(g 5 ).
The advantages of the EFT approach are even clearer at the next order [142]. For a complete determination of the O(g 5 ) effective potential in this theory, one needs to compute three-loop diagrams. However, only ordinary d-dimensional integrals within the EFT are needed at this order, and no new sum-integrals, i.e. Eq. (A.26) needs extending but not Eqs. (A.28) to (A.30).