Gravitational waves and tadpole resummation: Efficient and easy convergence of finite temperature QFT

We demonstrate analytically and numerically that"optimized partial dressing"(OPD) thermal mass resummation, which uses gap equation solutions inserted into the tadpole, efficiently tames finite-temperature perturbation theory calculations of the effective thermal potential, without necessitating use of the high-temperature approximation. An analytical estimate of the scale dependence for OPD resummation, standard Parwani resummation (Daisy resummation), and dimensional reduction shows that OPD has similar scale dependence to dimensional reduction, greatly improving Parwani resummation. We also elucidate how to construct and solve the gap equation for realistic numerical calculations, and demonstrate OPD's improved accuracy for a toy scalar model. OPD's improved accuracy is most physically significant when the high-temperature approximation breaks down, rendering dimensional reduction unusable and Parwani resummation highly inaccurate, with the latter underestimating the maximal gravitational wave amplitude for the model by 2 orders of magnitude compared to OPD. Our work highlights the need to bring theoretical uncertainties under control even when analyzing broad features of a model. Given the simplicity of the OPD compared to two-loop dimensional reduction, as well as the ease with which this scheme handles departures from the high-temperature expansion, we argue this scheme has great potential in analyzing the parameter space of realistic beyond the Standard Model models.

Unfortunately, 4D perturbation theory converges slowly at finite temperature (for a recent discussion of the convergence see Ref [62]).Naively, the distribution function diverges for long wavelength modes [63].This issue is delayed if the theory is resummed such that the strongly coupled, long wavelength behavior is screened by hard thermal loops [64][65][66].However, even after such a "Daisy resummation" the perturbation theory converges slowly and as a result a one-loop calculation can predict a peak gravitational wave amplitude that varies by mul-tiple orders of magnitude [14,67] for reasonable choices of the renormalization scale.
A known solution is to integrate out heavy Matsubara modes, which results in a 3D effective field theory.If one defines such a theory, both in its matching and the expansion in the dimensionally reduced theory, at next-to-leading order (NLO) in the appropriate coupling constant, 1 the uncertainty in the gravitational wave spectrum can in some cases be reduced to the percent level [14].Further, recent work has demonstrated that such perturbative calculations in the dimensionally reduced theory reproduce 3D lattice results very well, indicating that the infrared problems at higher loops may not be numerically important [68,69]. 2 The difficulty in dimensional reduction is its tractability.For even the Standard Model effective theory, defining the effective potential at noticeably better accuracy than (convenient) standard 4D methods requires the calculation of O( 102 ) diagrams at finite temperature.While a recent package adds automation to this process [70], one still has to monitor whether it is appropriate to have multiple 1 There is something of an inconsistency in the literature as to what prescription corresponds to what order.We will use the convention in this paper that NLO is when a calculation is performed accurately in dimensional reduction to O(g 4 ), such that performing the resummation at NLO and calculating the effective potential at NLO within the effective theory corresponds to correctly defining the theory up to O(g 4 ). 2 Two caveats deserve to be mentioned.First the timescale of the phase transition had strong uncertainties in the NLO dimensionally reduced theory.Second, the power counting used is slightly different to what has been used before, expanding in the ratio of the quartic coupling to the gauge coupling, x ∼ λ/g 2 , rather than the gauge coupling g 2 .
dynamical fields and whether to use the soft or ultrasoft potential for different regions of the parameter space.Further, it is difficult to go beyond the high-temperature expansion in dimensional reduction as there is no longer a hierarchy between the lightest Matsubara mode and the soft scale to justify an effective field theory.This is an important limitation, since strong phase transitions require sizable couplings to the scalar undergoing the transition, which in turn leads to large field-space dependencies for particle masses and hence breakdowns in the high-temperature expansion.Therefore, it is strongly motivated to find more convenient 4D calculational methods that can achieve similar levels of accuracy.This would be of significant utility in examining the large theory space of beyond the Standard Model (BSM) scenarios with strong phase transitions, electroweak or hidden.
One candidate for such a resummation method, called "partial dressing", was first developed three decades ago [71] (see also [72][73][74][75]) in the context of simple a ϕ4 toy model.It represented a simple analytical way of resumming the most important higher-order corrections to the thermal propagator without double counting.More recently, this method was revisited and adapted for numerical application to the full SM with an additional scalar [76].Referred to as "optimized partial dressing" (OPD), this method was also shown to be easily applied outside of the high-temperature approximation, promising a more accurate treatment of strong phase transitions in this important regime.However, significant further work is required before OPD could become a gold standard for finite-temperature calculations.A systematic and rigorous study of the scheme's convergence and validity is outstanding.There are open questions on the detailed analytical construction of the gap equation as a function of scalar vacuum expectation value (VEV), and the method of its numerical solution. 3Finally, gauge bosons were not yet consistently included in the system of gap equations, and OPD's relationship to renomalization group (RG) improvement of the effective potential is unclear.
The OPD scheme consists of two steps.The dominant contributions from many higher-order diagrams are conveniently included in a gap equation for the thermal mass.The convenience arises from the fact that the diagrams do not need to be evaluated analytically, one merely needs to solve the gap equation, analytically in some approximations but generally numerically.It is also straight forward to handle cases where the hightemperature expansion breaks down as one can merely use the full one-loop thermal functions within the gap equation.The second step involves inserting the full thermal mass into the tadpole, rather than the full one-loop effective potential as it was rigorously demonstrated that this is the way to prevent double counting of higher order diagrams [71].Finally, missing diagrams are easily identified and can be added by hand.A point to note is that while resummed two-loop potentials have been evaluated without resorting to high-temperature expansion in [77], the current work differs from the former since it includes subleading diagrams via the correct OPD resummation procedure, in addition to dealing with full thermal functions.Apart from convenience and handling cases where the high-temperature expansion breaks down, it is interesting to note that a (nonstandard) dimensional reduction calculation with gap resummation finds a critical Higgs mass at which electroweak symmetry begins to become first order [78] as well as a critical end point in QCD at finite density [79], unlike standard dimensional reduction or 4D perturbative calculations.This suggests solving the gap equation may even be the missing ingredient to accurately characterize the SM phase diagram within a perturbative treatment.
In this work, we carefully study OPD in the context of two simple test models-a single-field ϕ 4 theory, and a two-field ϕ 4 scalar field theory with a discrete Z 2 × Z 2 symmetry, only one of which is broken.We first analytically show that an RG-improved effective OPD thermal potential is parametrically as accurate as the RG-improved potential calculated using two-loop dimensional reduction, which formally motivates fully developing the OPD calculation in this and subsequent works.We then clarify construction of the gap equation and the method of its numerical solution, demonstrating that the corresponding effective potential improves on conventional Daisy resummation schemes [64][65][66] in accuracy and precision.Specifically, we compute the thermal potential near the critical temperature, compute the nucleation temperature and strength of the gravitational wave signal, and evaluate the dependence of these predictions on the renormalization scale to assess the degree of convergence.Unlike dimensional reduction, OPD is applicable in physical regimes where the high-temperature approximation does not apply, but when applied to scenarios where this approximation is valid, OPD yields results that are very close to the two-loop dimensional reduction calculation.This is in line with expectations derived from our analytical results. 4 Our simple test calculation can be seen as either obtaining results for a particularly simple hidden sector that undergoes a strong phase transition, or as a toy model for SM extensions with extra scalars that achieve a strong electroweak phase transition.At any rate, obtaining careful results for this simple scenario allowed us to improve and more completely understand several analytical and numerical aspects of the OPD method, setting the stage for future work to further develop OPD to include RG-improvement, the contributions of fermions and gauge bosons, momentum dependent contributions to the gap equations and nonrenormalizable operators in the Lagrangian.
The structure of this paper is as follows.In Secs.II and III, we review the different resummation schemes in perturbative finite temperature field theory and computation of the gravitational wave signal.In Sec.IV, we study analytically the scale dependence of physical predictions in the ϕ 4 theory and the two scalar field theory.The numerical implementation and analysis, including detailed construction of the gap equation away from the origin, is discussed in Sec.V. We finally conclude our paper with a discussion in Sec.VI.

II. PERTURBATION THEORY AT FINITE TEMPERATURE AND RESUMMATION METHODS
In this section we briefly review finite-temperature perturbation theory and the different resummation schemes we compare in this paper.
The form of gravitational wave spectra generated by a cosmological first-order phase transition ends up being quite sensitive to the precise description of the effective potential.For example, in the case of the Standard Model effective field theory or simple SM extensions with scalars, in regions where the theory has a strong phase transition the uncertainty in the peak gravitational wave amplitude can be multiple orders of magnitude [14,67].
The large uncertainty can be understood in two steps.First, finite-temperature two-loop (and sometimes higher-order) contributions can be similar in size to zerotemperature one-loop pieces.Second, any modest uncertainty in the temperature at which a phase transition occurs is amplified substantially in the prediction of gravitational wave observables.To see the second point, consider that the peak amplitude of a gravitational wave spectrum from a cosmic phase transition (see Sec. III) scales as where α ∼ ∆V /T 4 is the change in the trace anomaly and β/H * is the inverse timescale of the transition [80][81][82][83].The difference in pressure between the two phases, ∆V , tends to decrease with temperature and, from dimensional analysis, may typically scale quartically with the inverse temperature.The inverse timescale typically scales exponentially, so if the uncertainty in T is not too large one can take this as scaling linearly with the inverse temperature.Altogether, this accounts for a remarkable scaling relationship for the peak amplitude such that an O(25%) uncertainty in the percolation temperature can lead to two orders of magnitude uncertainty in gravitational wave observables.The parameters responsible for producing a gravitational wave amplitude can independently constrained by measurements at future colliders.However, if the gravitational wave observables are to provide meaningful constraints on the parameters, the theoretical uncertainties need to be brought under control.In fact, as we will later see, even broad questions about a model such as "what is the largest possible gravitational wave amplitude consistent with this model across its parameter space?" can vary by orders of magnitude comparing crude and more sophisticated calculation techniques.
We follow standard procedure in attempting to measure the importance of neglected higher-order terms by measuring the renormalization scale dependence of various observables [14].In doing so one needs to make a somewhat arbitrary choice of what range of values for the renormalization scale one should use in the loop calculation.Even for the generous range of scale parameter variation considered in [67], it is possible for the next-toleading order predictions to lie outside the uncertainty bands of the leading order prediction.However, the scale variation still gives a good parametric estimate of the size of the missing NLO terms.The conventional resummation methods that result in these large theoretical uncertainties were developed by Parwani [64] as well as Arnold and Espinosa [65,66].In these schemes, the long distance behavior is screened by including a thermal mass term, Π, such that the resummed potential has the form The thermal mass and the finite-temperature one-loop potential are readily calculated using standard methods. 5t is useful to recall the form of Π in the high-temperature expansion for a ϕ 4 theory with quartic coupling g: Both of these traditional resummation schemes only resum the leading O(T 2 ) term in the above thermal mass expansion, but they differ in some details.Arnold-Espinosa adds a ring term to the effective potential, which can be found by performing the substitution in Eq. ( 3) in the high-temperature-expanded thermal potential only.Parwani, on the other hand, substitutes M 2 + Π into the full effective potential, including the zero-temperature Coleman-Weinberg potential.For reasons that are beyond the scope of this quick discussion, the field tended to use the Arnold-Espinosa scheme, but the Parwani scheme actually has better scale dependence, since the thermal mass contribution in the zerotemperature Coleman-Weinberg piece induces a partial cancellation with the finite-temperature potential piece.We therefore use the Parwani scheme [referred to as truncated full dressing (TFD) in [76]] to minimize the scale dependence of traditional Daisy resummation, leading to the most conservative assessment of the benefits of OPD or DR.The Parwani resummation scheme correctly reorganizes the theory such that all pieces up to third order in the SU (2) L gauge coupling or larger are included in the potential.However, g 2 is reasonably large in the Standard Model, and terms at least in the next order O(g 4  2 ) are needed to bring uncertainties in the gravitational wave amplitude under control.Organizing perturbation theory with Parwani resummation to include all terms of O(g 4  2 ) is highly nontrivial, as naive methods lead to double counting.
To improve the accuracy of finite temperature perturbation theory, dimensional reduction at next to next-toleading order appears to provide a recipe, reducing theoretical uncertainties to the percent level [14].Dimensional reduction relies on the observation that in imaginary time, the quantum field theory is identical to a three-dimensional theory with a compactified time di-mension whose size is determined by the temperature [85].One can then integrate out the heavy Kaluza-Klein like modes, known as Matsubara modes, leaving an effective field theory in three dimensions.If a scale hierarchy persists between the remaining states, typically a soft and ultrasoft scale of order gT and g 2 T /π respectively, the soft states can be integrated out leaving behind a simpler effective field theory again.
Dimensional reduction naturally incorporates resummation through the matching relations between the four dimensional theory and the effective dimensionally reduced theory.Calculating the relevant self-energy diagrams in these matching relations can be done at multiloop levels of accuracy in addition to defining the effective potential at multiple loops within the effective theory without double counting.The whole process makes it straightforward albeit work intensive to organize perturbation theory into powers of an effective coupling constant (or a ratio of constants as argued for in ref. [69]).At next-to-leading order [or O(g 4 )] the theoretical uncertainties in the gravitational wave amplitude for the Standard Model effective field theory are effectively under control.The drawback of dimensional reduction is both the enormous practical difficulty of the schemeover a hundred diagrams are required to see noticeable improvement compared with conventional methods-and the fact that it conventionally assumes the validity of the high temperature expansion.
Another method of resummation-partial dressing or gap resummation calculates many higher-order diagrams numerically by solving a gap equation, which can be pictorially represented as follows: where the double lines represent the resummed mass and the single line the tree-level mass.It was demonstrated three decades ago [71] that in order to avoid double counting of higher-order diagrams, one needs to insert the resummed mass into the tadpole and integrate, where the first term, V 1-loop , is the zero plus finitetemperature one-loop potential and the second term is the finite-temperature potential due to two loop sunsets.This term is neglected in the solution to the gap equation.The method can also be applied when the high-temperature regime breaks down, in particular by keeping the high-temperature approximation in the gap equation but using full thermal integrals in the potential [76] (referred to as OPD), and is much easier to use than dimensional reduction as the gap equation only requires one to define the one-loop effective potential at finite temperature.It is therefore highly attractive if it can be demonstrated to provide substantial improvement over conventional methods, i.e.Parwani resummation.
In this work, we seek to ascertain systematically for the first time whether gap resummation, or OPD, performs better than Parwani resummation and is comparable to dimensional reduction at next-to-leading order for the purposes of obtaining gravitational wave predictions.We also clarify how the gap equation should be constructed and numerically solved away from the origin.

III. REVIEW OF GRAVITATIONAL WAVE SIGNAL CALCULATION
Any uncertainties in the prediction of the thermal potential, arising from the slow convergence of perturbation theory at finite temperature, become amplified when calculating the gravitational wave observables.In this section we will review the calculation of the gravitational wave peak amplitude due to the acoustic contribution calculated using the sound shell model [86][87][88][89].We will later consider the two scalar field model and assume that the Standard Model particle content is coupled strongly enough to maintain kinetic equilibrium with our two new scalar fields but weakly enough that we can ignore their effect on the gravitational wave phenomenology.In the sound shell model, the spectrum is completely determined by the temperature of the transition, the fluid velocity, the mean bubble separation, the life time of the sound waves and the fraction of energy released that becomes converted to sound waves.To a good approximation, these quantities can be related to four macroscopic parameters -the transition temperature, the bubble wall velocity, the trace anomaly normalized by the critical density and the inverse lifetime of the transition A more diligent calculation irons out errors that are orthogonal to our analysis [90] so we ignore them in this work.It should be noted, however that in the sound shell model, the four thermal parameters result in two observables, the peak amplitude and frequency.More careful simulations yield a more optimistic picture where all four thermal parameters can be extracted from the precise shape of the power spectrum [91], but conservative analyses claim fits to three parameters [92].In any case, the sound shell model predicts a spectral shape, of the form with a peak frequency Hz, (9) and a peak amplitude [82,[86][87][88][89]] where v w is the bubble wall velocity and g * counts the number of relativistic degrees of freedom in the plasma at the time of transition.For the efficiency we assume a relativistic bubble wall velocity which means we can, to a good approximation, relate the trace anomaly to the kinetic energy fraction and the suppression factor from the finite lifetime of the source is where τ sw = R * / Ūf with the fluid velocity and the mean bubble separation having the form U 2 f ∼ 3 4 κα and R * = (8π) 1/3 v w /β, respectively.Finally, to derive the percolation temperature T p one has to solve the equation, IV. ANALYTIC COMPARISON OF RESUMMATION METHODS FOR ϕ 4 THEORIES Considering single-ϕ 4 and two-ϕ 4 theories with Z 2 symmetry and working in the MS renormalization scheme, we now follow a similar procedure as [67] in the high-temperature limit to analytically show that the effective potential computed in an RG-improved one-loop OPD calculation has parametrically similar scale dependence to the RG-improved two-loop dimensional reduction calculation, greatly improving on standard Parwani resummation.This motivates full development of the OPD calculation, including the progress we make in constructing and consistently solving the gap equation and comparing scale dependence across different calculation schemes in the numerical analysis of Sec.V.

A. Single-field ϕ 4 theory
Let us consider a simple real ϕ 4 theory with a discrete Z 2 symmetry.This theory admits phase transition that is second-order and therefore cannot be considered for gravitational wave phenomenology.However, this simplicity also makes the differences between the resummation schemes in obtaining predictions for the thermal potential transparent, and will serve as a useful pedagogical warmup.The tree-level potential has the form, The one-loop correction at zero temperature in addition to the tree-level potential gives the Coleman-Weinberg potential, where µ is the renormalization scale in the MS scheme and M 2 = m 2 + g 2 ϕ 2 /2 is the field-dependent mass.No physical quantity will depend on the choice of the renormalization scale.Therefore, the above dependence is merely an artefact of the truncation of perturbation theory at finite order.As such, a convenient measure for the convergence of perturbation theory is the residual scale dependence -the implicit scale dependence arising from running couplings in the tree-level potential should cancel the explicit scale dependence in the Coleman-Weinberg potential leaving corrections from the implicit scale dependence in the latter term which is formally higher order.Let us demonstrate this explicitly.Our renormalization group equations at one loop have the form The scale dependence of our potential to O(g 4 /π 2 ), dropping two-loop-size terms, then has the form which sums to a field independent cosmological constant that we can ignore.At finite temperature, the one-loop correction in the high-temperature expansion has the form, where a b = 16π 2 exp 3 2 − 2γ E , γ E being the Euler-Mascheroni constant.It is trivial to see that nothing in our one-loop theory cancels the implicit scale dependence of the finite-temperature piece.Moreover, the size of the uncanceled corrections are as large as the scale dependence of the tree-level pieces in powers of g/π and could be relatively larger when T /ϕ is large, that is the infrared limit.On top of the explicit infrared divergence that appears at higher-loop, finite-temperature perturbation theory converges slowly in part because of a mismatch of the order of the loop expansion and the size of a term in powers of g n /π m as well as the infrared enhancement of uncanceled pieces in terms of T /ϕ.

Parwani resummation
Historically, resumming daisy diagrams was to cure the infrared divergences inherent in finite-temperature field theory [64].However, since this means including higherloop diagrams in the effective potential, we will also see an improvement in the scale dependence.Parwani resummation works through replacing M 2 → M 2 + Π where Π = 1 24 g 2 T 2 is the thermal mass to lowest order in the high-temperature expansion.A partial cancellation of the scale dependence occurs The opposite signs in the above terms is the origin of the cancellation.To achieve a full cancellation, we require the missing two-loop term which is the sunset diagram.At high temperature it is In the above we have performed a high-temperature expansion.The full expression that we use in our numerical calculation we put into Appendix B. A straightforward calculation shows that the scale dependence of the T 2 ϕ 2 term cancels.The leading order uncanceled piece is

Dimensional reduction
At finite temperature in the imaginary time formalism, a thermal field theory is equivalent to a Kaluza-Klein theory with a compactified imaginary time direction of size 1/T .Integrating out this tower of Matsubara modes leaves one with an effective theory in three dimensions.Doing so automatically includes resummation by construction via the matching relations.Although famously rigorous and consistent, it can be quite formidable technically in a realistic model.In the case of ϕ 4 theory, the resulting effective potential is simple enough to be written in a closed form, even at NLO, where Here, µ 3 is the scale dependence in the effective theory and For a fair comparison, we take the larger of the dependencies on µ 3 and µ as reflective of the residual scale dependence at NLO.To calculate the dependence on the former, it is useful to write The dependence on µ turns out to be subdominant and has the form so we focus on the µ 3 dependence and find This is of order g 4 /π 3 , much smaller than the O(g 2 ) residual dependence of the Parwani calculation.

Gap resummation
Higher-order diagrams can be included in a simple gap equation where the thermal mass is defined as the second derivative of the resummed potential The resulting thermal mass cannot be substituted into the full potential without double counting diagrams.Instead, one includes the resummation of the tadpole , and then integrates over the field to acquire the potential at the end.In the above, everything was written in the high-temperature limit to make the problem analytically tractable.However, the power of gap resummation is that the gap equation can be solved numerically without a high-temperature expansion.The scale dependence of the resummed mass and the tadpole can be written to O(g 4 ), After substituting the derivative of the gap equation with respect to the scale into the scale dependence of the tadpole, the remaining scale dependence cancels up to O(g 4 ) and the residual piece is O(g 6 ), where ζ = 0 corresponds to the direct result of the above treatment, and ζ = 1 is obtained by including sunsets in the gap equation.Of course, this is the scale dependence of the tadpole, not the potential.After integration, the residual term is the same order as in DR with a slightly different prefactor, suggesting that the use of gap equations is competitive with NLO DR, where we have, in the final term, expanded an arctanh function that is actually well-behaved in the infrared limit.

B. Two-field ϕ 4 theory
Let us now consider a two-scalar field model.In principle, this is the minimal model that could produce an ob-servable gravitational wave signature, as the portal couplings can produce a modest thermal barrier.In practice, the peak amplitude tends to be very small.Nevertheless, the predictions of this model can be treated as realistic phenomenological predictions, perhaps existing in some dark sector, under the proviso that any couplings keeping the system in kinetic equilibrium with the visible sector can be sufficiently small that their effect on the potential is negligible.The potential for our model is We will consider the case where the second scalar does not acquire a VEV throughout the transition, as its function is to provide the thermal barrier.In this case the field dependent masses have the simple form Note that we include ϕ 2 when performing derivatives with respect to the potential and only set ϕ 2 to zero at the end.To keep equations compact, we do not show their contribution here.Finally, the renormalization group equations have the form 1. Parwani resummation In the Parwani scheme, the thermal masses have the form Let us now put together the scale dependence of each piece of the potential in the Parwani regime.First the tree-level potential, (51) and the zero-temperature piece of the Coleman Weinberg has the form (52) The above piece cancels the field dependent part of the scale dependence of the tree-level potential.
The quadratic-temperature dependent piece has three parts.First from the high-temperature expansion of the one-loop potential, which at lowest order in the hightemperature expansion is This partially cancels the piece arising from the thermal masses in the Coleman Weinberg potential Adding the two terms together leads to a residual piece with Mi 2 = M 2 i + Π i , and we also show the next-toleading O(T ) part of the high-temperature expansion for reasons that will become clear.The final contribution is from the leading power sunset term, which in the high-temperature expansion has the form The sunset actually cancels the O(T 2 ) term in Eq. ( 55) completely, so we focus on the O(T ) contribution The residual scale dependence in V Parwani is the sum of the O(T ) terms in Eqs. ( 55) and (57).If there is strong first-order phase transition, we expect g 12 ≫ g 1 , g 2 .We therefore show the leading order field-dependent term in powers of g 12 : The leading order uncanceled term is therefore of the same order as in the single-field case.

Dimensional reduction
We use DRalgo [70] to derive the potential for this model at O(g 4 ).That is we include two-loop calculations in the dimensionally reduced theory and a NLO (two-loop) resummation (NNLO in the nomenclature of DRalgo).For our analytic comparison, it is easiest to work with the soft rather than the ultrasoft potential, though the appropriate potential to use depends upon where one is in the parameter space.The full soft potential is Here the dimensionally reduced couplings, masses and fields are given by respectively.Again, µ 3 is the renormalization scale in the three-dimensional theory and A ∼ 1.2 is the Glaisher number.It is straight forward to show that the residual scale dependence of (λ 1,3d , λ 2,3d , λ 12,3d ) is of O(g 6 i ).The masses cancel at the same order with the exceptions The residual scale dependence then for NLO DR is

V. NUMERICAL IMPLEMENTATION OF GAP RESUMMATION AND RESULTS
We now review how to set up the numerical OPD calculation.In particular, we specify that an iterative approach should be used without any mass derivatives in the gap equation, clarifying some ambiguities from the original numerical treatment [76].An important point to note is that in our numerical implementation we will not rely on high-temperature approximation and will be solving the gap equation away from the origin, neither of which was done in the previous section.We then present numerical results for the thermal potential and its gravitational wave observables for a representative benchmark point in the two-field ϕ 4 theory, comparing OPD and Parwani resummation to demonstrate the improved accuracy and precision of the OPD calculation.To facilitate comparison to dimensional reduction, which only works when the high-temperature approximation is valid, we also compare effective potentials for a scenario that does not yield a strong phase transition, showing that the OPD calculation yields a result very close to the dimensional reduction calculation, both of which differ from the Parwani result.This supports the expectation from our analytical derivation that one-loop OPD is competitive with two-loop dimensional reduction, and motivates development of RG improvement for OPD to fully realize the OPD accuracy promised by the analytical results of the previous section.(As explained in the introduction, these numerical studies all utilize effective potentials without RG improvement, as formulating a consistent RG-improvement scheme for OPD that is not restricted to the high-temperature approximation is the subject of upcoming work.)

A. Construction and solution of gap equation away from origin
Our procedure for OPD resummation is as follows.Note we restrict ourselves here to the scenario where only one scalar acquires a VEV during the phase transition: i instead of Π i to denote corrections to each scalar field's mass beyond the tree-level m 2 i , since it in general includes both zero-and finitetemperature corrections.For a given temperature, the mass corrections δm 2 i are obtained by numerically solving a set of coupled algebraic gap equations on a grid of field values along the excursion of the symmetry breaking field, in this case ϕ 1 , • The continuous functions δm 2 j (ϕ 1 , T ) are obtained by interpolating the solutions to the gap equation on a grid of ϕ 1 VEVs, which are then substituted in the first derivative of the zero + finite-temperature one-loop potential, plus the finite-temperature twoloop sunset term.
One of the crucial aspects of OPD is its numerical efficiency while going beyond the high-temperature approximation, achieved by using full thermal integrals in the potential but high-temperature approximations in the gap equation. 6This was justified in [76] by arguing that the gap equation only matters for phase transitions in regions of field-and parameter-space where components of the plasma become light and hence the high-temperature expansion is valid, while further out in field space where the high-temperature expansion fails the thermal mass is accurately treated as small in the full potential, meaning the gap equation has a much smaller effect and its error is parametrically suppressed.As we outline in the next subsection, we have systematically verified that this assumption is in fact correct, justifying the use of the high-temperature expansion in the gap equation and the enormous simplification it brings.The algebraic equation ( 76) can be solved iteratively: for a fixed ϕ 1 and T where the solution starts at δm 2 j (ϕ 1 , T ) n=1 = 0 and converges to a fixed value after a handful of iterations.
The above is almost identical to the numerical OPD procedure originally laid out in [76], with one impor-6 Note that the full thermal integrals can be efficiently approximated to high precision by a piecewise defined function joining the high-and low-T approximations, which for the bosonic case is given by J piecewise n 2 K 2 (y n) for y 2 less [more] than 0.22, where K 2 is the modified Bessel function of the second kind.We limit n ≤ 3 but more precision is easily obtained by including more terms.tant difference.The authors of [76] found that the purely algebraic, iterative method of solution is problematic, since it yielded multiple oscillating solutions to the gap equations away from the origin in field space when applied to the SM with simple scalar extensions.This was solved in an ad hoc manner, by using an alternative formulation of gap equations which involved keeping mass derivatives in the gap equation and using them to constrain δm j (ϕ 1 + ∆ϕ 1 , T ) based on the previous solution δm j (ϕ 1 , T ) on the ϕ 1 -grid, turning the algebraic gap equations into a set of differential equations.While this yielded unique and apparently reasonable solutions most of the time, the procedure was very vulnerable to numerical errors due to the singular nature of the resulting gap equation near field values where ∂ 2 V /∂ϕ 2  1 flips sign, i.e.I: Physical and MS parameters of the two-field model for our benchmark numerical calculation at the renormalization scale µ = 600 GeV (check Appendix A for the relation between the two).The fixed value of the renormalization scale is only used to plot Figs. 1 and 2.
when passing through the potential barrier at the critical temperature.The inclusion of derivative terms was also not strictly justified by the original proof of the validity of partial dressing in [71], even though their effect was expected to be subleading.
In our careful analysis of OPD as applied to the much simpler two-field ϕ 4 theory, we found that this differential version of the gap equation yielded numerical solutions of the effective potential that contained unacceptable artifacts for the sizeable couplings that yield a first-order phase transition.On the other hand, the simpler iterative approach always yielded unique and reasonable solutions to the system of gap equations.We therefore use the solution method of Eq. ( 78) in our analysis.We hypothesize that the nonconvergence of the gap equation solution in [76] was caused by applying OPD to the full SM with extra scalars, without consistently including gauge bosons in the system of gap equations (rather just including their O(T 2 ) contributions in the scalar gap equations).
Our purely local and iterative implementation of the gap equation is therefore fully consistent with the original formulation of partial dressing [71].This is ultimately a fortunate development for the practical application of OPD, as the purely iterative numerical implmentation is much simpler and faster.We will show how to consistently expand OPD to include gauge bosons in an upcoming publication.
B. Numerical results for benchmark two-field ϕ 4  theory We will now consider the two-field ϕ 4 theory with benchmark parameters shown in Table I, numerically computing the effective potential and gravitational wave signal in Parwani and OPD resummation to compare the two schemes.We find similar behaviour for other parameter points with strong phase transitions, so these results are representative.
This numerical study is done by choosing the MS parameters at zero temperature for an arbitrary renormalization scale µ (chosen to lie near the ϕ 2 mass) which reproduces experimental observables at one-loop and using these parameters to compute the thermal potential and gravitational signal.This procedure is then repeated for different choices of µ to estimate the scale dependence and hence theoretical uncertainty of our results [67].The experimental observables for our model are the vev ⟨ϕ 1 ⟩, pole mass M 1 , M 2 of the scalar fields at the symmetry broken vacuum and the quartic couplings g 2 12 and g 2 2 .The details of the one-loop matching of the parameters to the experimental observables are given in Appendix A.
Figure 1 shows the thermal potential for our benchmark point at T = T c when the true and false minima are degenerate, as well as the corresponding solutions to the gap equation.As one can see the iterative method does give a smooth, unique solution exhibiting correct physical behavior, whereby the mass corrections are maximum at the origin and decreases with ⟨ϕ 1 ⟩ since the fields acquire more mass reducing their participation in the thermal plasma.This is markedly different from the constant thermal masses assumed in Parwani resummation, which makes use of the lowest-order high-temperature expansion far away from the origin where it is no longer justi- fied.
On the other hand, we also verified, for this and other choices of parameters, that the use of the high-T expansion in the gap equation for OPD was valid.Compared to solving the gap equations with full thermal functions, we only found meaningful differences in the region of field space where M 2 /T 2 is large, M 2 being the resummed mass.For example, for this particular benchmark this difference shows up when M 2 2 /T 2 ≳ 3, significantly larger than the still sizeable but more moderate values relevant for our phase transition.Even then, the difference in δm 2 i (full thermal potential) obtained with the high-temperature approximation in the gap equation vs full thermal functions in the gap equations are at most ∼ 10% (1%).This confirms the original argument for the high-temperature approximation in the gap equation made in [76].
OPD by construction is more accurate than Parwani, as it takes into account both proper counting of diagrams and higher-order corrections.These higher-order corrections do make a significant change in the profile of the thermal potential.OPD seems to predict a larger barrier for the thermal potential when compared to Parwani as seen in Fig. 1.This tends to be a general feature of OPD throughout the parameter space and results in the maximal gravitational wave amplitude being orders of magnitude larger than what is predicted by Parwani as can be seen in Fig. 2.
In addition to better accuracy, the numerical results also demonstrates significantly improved precision, as signaled by the reduction of scale dependence when using OPD.Figures 3 and 4 show the variation of the thermal parameters and peak gravitational wave amplitude as a function of the renormalization scale.An estimate of the variation of the critical and percolation temperature , ∆T = Tmax−Tmin Tmax on varying the renormalization scale gives ∆T OPD ∼ 2% compared to ∆T Parwani ∼ 15%.Similarly, almost a factor of 1.7 reduction is observed in the variation of the strength of the phase transition, ∆α with ∆α Parwani ∼ 49% and ∆α OPD ∼ 29%.Of course, it is to be noted that not all the thermal observables show such an improvement.In particular it can be seen in Fig. 3 that ϕ c has a larger variation in OPD (∆ϕ c,OPD ∼ 12.6%) compared to Parwani (∆ϕ c,Parwani ∼ 3.4%).This on the other hand has minimal effect on variation of the gravitational wave observables since both α and β/H * have stronger dependence on T p rather than ϕ c .In particular, this also means that OPD's prediction of the gravitational wave peak frequency has much better precision.While the results here show the superiority of OPD compared to Parwani, one should note that in a realistic model one would expect even better scale dependence.This is due to the fact that to achieve a first-order phase transition in this model, a very large g 12 coupling is required to compensate for the small number of degrees of freedom that become massive during the transition.The portal coupling is so large that the uncertainty from matching is significant, despite this being a zerotemperature uncertainty.This is an unfortunate artifact of the toy model we used to develop this analysis.The µ dependent lines for the thermal parameters are almost in parallel, demonstrating that as we are near the nonperturbative regime where zero-temperature uncertainties are nearly out of control, even if the finite-temperature uncertainties are greatly improved by the OPD resummation scheme.Therefore, we expect the gravitational wave prediction to be significantly improved in OPD for more realistic models; furthermore, future work to extend OPD to consistently include RG improvement of the effective potential, and possibly including momentum-dependent terms in the thermal mass resummation, would allow OPD to give more accurate gravitational wave predictions for this particular toy model as well.
Finally, as a last check we should also compare the results to that of dimensional reduction.Unfortunately due to the large couplings involved the high-temperature expansion for M 2 badly breaks down near the critical temperature.This makes the use of standard forms of dimensional reduction inappropriate.Particularly problematic is the two-loop sunset diagram that scales as ϕ 2 T 2 times a large logarithm in the high-temperature expansion (see the last line of equation 59).This term actually dominates when not Boltzmann suppressed and is responsible for the double barrier behaviour we see in Fig. 5.Note the ultrasoft potential does not do better as it is actually unbounded (and contains strange kinks).For this model, the high temperature expansion -and therefore dimensional reduction -breaks down for any couplings large enough to catalyze a strong first-order phase transition, so it makes little sense to compare the predicted thermal parameters quantitatively.In this comparison, our calculation gives the best result by default, and while this is physically significant in terms of OPD expanding our capacity to make reliable predictions for phase transitions that are not in the high-temperature regime, it is instructive to consider a case where a quantitative comparison to DR at O(g 4  12,phys ) can be made.We therefore compare the potentials directly for differing values of the portal coupling g 12 , without requiring a strong phase transition.We expect that, for small to medium sized g 12 , OPD and DR should agree with each other and disagree with Parwani; but at some threshold DR and Parwani, which both depend upon the high-temperature expansion, should diverge from OPD, which does not.We find this to be the case, as illustrated by three benchmarks in Fig. 6.In fact, when the high-temperature expansion breaks down, DR gives a similar prediction to the Parwani calculation.This reinforces the two main points of our analysis: first, as demonstrated analytically in Sec.IV, that OPD improves the four-dimensional perturbative calculation to the point where it is of similar accuracy as DR, and second, that OPD improves upon DR by not relying on the high-temperature approximation, which is not respected in many physically relevant phase transitions.

VI. DISCUSSION AND CONCLUSION
Understanding the era of cosmological electroweak symmetry is a central question of the next generation of theoretical and experimental efforts.Conventional  methods of modelling cosmological electroweak symmetry breaking suffer very difficult theoretical uncertainties.Accurate calculations involving next-to-leading order dimensional reduction are very difficult and, to this day, there exists no global treatment of a BSM scenario involving extra dynamical scalar fields.Further, off the shelf dimensional reduction is in the high-temperature regime, whereas strong phase transitions catalyzed by thermally induced barriers necessitates large couplings where the high-temperature regime is expected to be invalid.
In the context of all of these issues, we are motivated to investigate and further develop the OPD scheme utilizing a recursive solution to the gap equation for scalar masses (though it is worth keeping in mind that the solution of the gap equation can also be used as a replacement for calculating the self energy used in matching relations in the dimensional reduction paradigm as well).OPD includes sizeable contributions neglected by the Parwani scheme, and is therefore required for more accurate calculations.In the context of theoretical uncertainties, our results seem to suggest regions for cautious optimism with regard to OPD.The analytic calculation in Sec.IV indicates a similar precision, i.e. scale dependence at fourth coupling order, to dimensional reduction, improving on the second-order dependence of Parwani resummation, while the ease of going beyond the hightemperature approximation in OPD promises additional advantages for phase transitions with sizeable couplings compared to DR. Numerical results in Sec.V for the thermal parameters show an improvement for OPD compared to Parwani, giving further grounds for optimism.All of the thermal parameters except the critical vev display significant reduction in theoretical uncertainty when using OPD.This anomalous behavior for the critical vev does not significantly influence the final gravitational result due to larger dependence on the other thermal parameters, and at any rate we suspect that this is an artifact of the very large couplings necessary for a strong FIG.6: OPD (black) and Parwani (dashed) potentials for a different benchmark from Table .I, where the high-T approximation is at least close to valid, in order to quantitatively compare to DR (purple).⟨ϕ 1 ⟩ = 400 GeV, M 1,pole = 80 GeV, M 2,pole = 500 GeV, g 2 2,phys = 3.60, µ = 500 GeV with g 2 12,phys = 0.30, 1.78, 2.26 (top, middle, bottom).Crucially, OPD and DR agree for the first two portal couplings, demonstrating the improvement of OPD compared to Parwani.At higher g 2 12 (bottom), the high-T approximation used by DR and Parwani breaks down, and only OPD can remain reliable.
phase transition in our toy model.We also expect this behaviour to improve once a consistent RG-improvement scheme for the OPD effective potential is implemented.
In summary, both analytical and numerical results clearly indicate the superiority of OPD to the Parwani scheme, which miscounts subleading corrections that are numerically important at large coupling and nonzero scalar vevs, both generally found in phase transitions.Compared to DR, OPD should have parametrically similar theoretical uncertainty once RG-improvement is implemented, with the further advantage of working beyond the high-temperature regime, which is again ubiquitous for strong first-order phase transition.
Our study represents a first step in a systematic and rigorous development of the OPD scheme for precision high-temperature calculations, and there is a clear itinerary of future directions that must be pursued to apply OPD to realistic extensions of the SM.The highest priority next steps include a careful examination of how OPD could be combined with renormalization group improvement of the effective potential to further reduce the theoretical uncertainties; understanding the importance of full momentum dependent self-energy in gap equations; consistent inclusion of gauge bosons in the system of gap equations, and addition of sunset equivalent diagrams in the gauge sector; as well as appropriate modification to the gap equation for multiple symmetry broken fields to handle arbitrary field excursions during the phase transition.Finally, being able to include nonrenormalizable operators would be very useful, given the importance of SMEFT to constrain BSM physics at colliders and elsewhere.We are currently working to address these issues and plan to present them in future publications.

FIG. 1 :
FIG. 1: Upper panel: The effective potential with one-loop matching for the benchmark in Table I and µ = 600 GeV with the solid line denoting OPD and the dashed line referring to Parwani resummation.Lower panel: Numerical solution of the gap equation for the mass correction δm 2 i (ϕ 1 , T c ) at the critical temperature for the same benchmark.Here T c = 287 GeV for Parwani and T c = 265.6GeV for OPD.The orange and blue curves are for ϕ 1 and ϕ 2 respectively.

- 20 FIG. 2 :
FIG.2:The maximal gravitational wave amplitude for all parameters fixed except g 12 .Solid line denotes the prediction of OPD and dashed line the predictions of Parwani.

FIG. 3 :
FIG.3:The scale dependence (µ) of the critical temperature and VEV for the benchmark in TableIwith dashed line corresponding to Parwani and solid line corresponding to OPD with both augmented by two-loop sunsets.

- 21 FIG. 4 :
FIG. 4: Scale dependence (µ) of the thermal parameters and the peak gravitational wave amplitude as predicted by the sound shell model.Solid/dashed lines correspond to OPD/Parwani, respectively.

FIG. 5 :
FIG.5: Soft potential in NLO dimensional reduction using the benchmark in TableI(solid purple) at T = 191.7 GeV as well as the problematic HT calculation of the two-loop sunset which dominates (dashed purple).