Higher orders for cosmological phase transitions: a global study in a Yukawa model

We perform a state-of-the-art global study of the cosmological thermal histories of a simple Yukawa model, and find higher perturbative orders to be important for determining both the presence and strength of strong first-order phase transitions. Using high-temperature effective field theory, we calculate the free energy density of the model up to O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(y5T4), where y is the Yukawa coupling and T is the temperature. The locations of phase transitions are found using the results of lattice Monte-Carlo simulations, and the strength of first-order transitions are evaluated within perturbation theory, to 3-loop order. This is the first global study of any model at this order. Compared to a vanilla 1-loop analysis, accurate to O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(y2T4), reaching such accuracy enables on average a five-fold reduction in the relative uncertainty in the predicted critical temperature Tc, and an additional ∼ 50% strong first-order transitions with latent heat L/Tc4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ L/{T}_c^4 $$\end{document}> 0.1 to be identified in our scan.

1 Introduction Extensions of the Standard Model (SM) of particle physics can lead to additional phase transitions in the early universe, with consequences for open problems such as baryogensis and dark matter.This in turn can produce a gravitational wave background detectable at current and future detectors such as LISA [1], Taiji [2], TianQin [3], DECIGO [4], NANOGrav [5] and the International Pulsar Timing Array [6].For a sizeable gravitational wave background to be produced requires a first-order phase transition, unlike the smooth crossovers corresponding to electroweak symmetry breaking [7] and colour confinement [8] in the SM.
We study perhaps the simplest model giving rise to a first-order phase transition, by coupling a real scalar field (ϕ) to a Dirac fermion (ψ) via a Yukawa interaction, where, in the mostly minus metric and using the Feynman slash notation, / ∂ = γ µ ∂ µ , L Yukawa = −yϕ ψψ.(1.4)In principle, this Yukawa theory may be coupled to the SM through the Higgs portal.The resulting model has been widely studied as providing a simple dark matter candidate [9][10][11][12][13][14][15][16][17][18][19][20], and is sometimes known as singlet fermionic dark matter.Phase transitions in this model have been studied in refs.[21][22][23].However, in what follows, we will neglect Higgs portal couplings, and focus on the minimal model of ϕ and ψ, as our primary interest is the reliability and convergence of perturbative approaches to first-order phase transitions.Our model can still have more direct physical relevance in the context of feebly interacting massive particles [24,25].
At the high temperatures necessary for a phase transition, with πT large compared to m and m ψ , the usual loop expansion breaks down due to the hierarchy of scales introduced between the temperature and relevant masses.Whereas the mass of the zero Matsubara mode has thermal contributions only at the soft scale O(yT ), the nonzero Matsubara modes receive thermal contributions to their masses at the hard scale O(πT ).This scale hierarchy causes some Feynman diagrams to be larger than suggested by loop counting.Effective field theory, applied in the framework of high-temperature dimensional reduction [26][27][28][29][30], offers a systematic means to account for this through resummations, and the only framework which has been pushed to higher perturbative orders.For an introduction, see ref. [28].
Perturbative calculations of the thermodynamics of cosmological phase transitions have been observed to suffer from huge theoretical uncertainties in a variety of models, often amounting to a factor of between a hundred and a million for the predicted gravitational wave peak amplitude.Numerically, often the largest theoretical uncertainty appears as a strong dependence of physical observables on the renormalisation scale [31][32][33][34][35][36].Relatedly, computations using different approximations have been found to yield drastically different predictions [37,38].However, this uncertainty is significantly reduced with higher order dimensional reduction, demonstrating the importance of such calculations.
High-temperature dimensional reduction has been applied at 2-loop order to a wide range of models, including the Standard Model [29] and a number of its extensions (see for example refs.[39,40]), and at 3-loop and partial 4-loop order to QCD [30,41] and the Standard Model in the electroweak symmetric phase [42].The generic rules of dimensional reduction were put down in ref. [29], and the first software package to carry this out, DRalgo, was published recently [43].The dimensional reduction of the Yukawa model considered here has not previously been carried out.Our results in fact provided a correction for DRalgo, to be included in an upcoming release.
Motivated by making reliable predictions for gravitational wave experiments, recently there has been growing interest in quantifying and reducing theoretical uncertainties for cosmological phase transitions [31,34,[44][45][46][47][48][49].A number of recent studies have tested the validity and accuracy of perturbation theory, utilising dimensional reduction to push to higher orders [31,33,[50][51][52][53][54].While these studies showed that lower-order calculations can suffer from large theoretical uncertainties, this was based on studying a small number of benchmark parameter points.One is led to wonder how generic these conclusions are.For example, in a broad scan of the parameter space of a model, do the global features of the results change between successive perturbative orders; in short, does the blob move? 1 In this work, we aim to address this and related questions for our Yukawa model.

Dimensional reduction
Equilibrium thermodynamics of quantum field theories can be studied in the imaginary time formalism [55].The fields of the partition function are defined in three infinite spatial directions and one compact direction with length 1/T .In the compact direction, bosons, including the scalar ϕ in our model, satisfy periodic boundary conditions and can be expanded as where τ parameterise the compact direction and x the spatial directions.On the other hand, fermions such as ψ satisfy anti-periodic boundary conditions and can be expanded as These Fourier modes are referred to as Matsubara modes [56].The free part of the Euclidean equilibrium action can be written as where γ i and γ 0 are Euclidean gamma matrices [57].
Thus one can view the equilibrium thermodynamics of the four dimensional (4d) theory as the vacuum dynamics of infinitely many 3d fields [58], the fields φ n with effective squared masses m 2 n = m 2 +(2πT n) 2 and the fields ψ n with effective squared masses m 2 ψ +(2πT (n+1)) 2 .At temperature T the effective coupling constant of the n = 0 light bosonic mode increases as ∝ T /m 0 , reflecting its high occupancy n B ∼ T /m 0 .The naive perturbative expansion for this mode therefore breaks down at high temperatures.A resummed perturbative expansion is nevertheless possible, though its convergence is slower than at zero temperature.
At high temperatures some form of resummation is necessary, yet the specifics depend on what assumptions are made regarding the relative magnitudes of different parameters [59].In this paper, the following power counting prescriptions were adopted for the parameters of the Yukawa theory, fixing the relative sizes of the couplings such that: (i) the loop expansion parameters at zero temperature are of similar size, and (ii) the thermal contributions to the effective potential are of similar size to the tree-level terms, as expected for a phase transition (additional details given in appendix A): Using the hierarchy of scales between masses and the temperature, we can construct an effective field theory (EFT) for only the light modes.This construction is referred to as dimensional reduction due to the EFT being defined in three dimensions, and is the most general 3d theory obeying the same internal and spatial symmetries as the 4d theory, as well as having the same number of light bosonic degrees of freedom: (2.5) In principle, higher dimensional operators could be added to this effective Lagrangian, though these turn up at one higher order in y than we work, using the power counting prescription of equation (2.4) [28,52].At a first-order phase transition, the different terms in equation (2.5) balance to yield two minima with a maximum between.For this to be possible, the different terms in the potential must be the same order of magnitude, such as m 2 3 ϕ 2 3 ∼ λ 3 ϕ 4 3 .Extending the power counting to the 3d couplings, m 2  3 ∼ m 2 ∼ y 2 T 2 and λ 3 ∼ λT ∼ y 2 T , this relation leads to so that relatively strong phase transitions are possible in this model.Note that the powers of temperature arise through first scaling the 3d scalar field to canonically normalise it, and then absorbing all explicit factors of T into the parameters of the EFT [28,29].At high temperatures and at leading order (LO), the free energy density of the Yukawa theory is f = −π 2 T 4 /20, independently of the background scalar field.For first-order phase transitions, it is the difference between phases ∆f which determines the dynamics.For homogeneous phases, the free energy density is equal to the thermal effective potential, and its expansion takes the form ∆f where the superscript + indicates that these terms require resummation to compute.The coefficients c n are all of O(1) in our power counting, and are temperature dependent.In this work, we have evaluated the free energy density up to the O(y 5 T 4 ) term, i.e. calculating up to c 5 .This is then compared to the result at leading order, with the free energy density at O(y 2 T 4 ), i.e. just c 2 .For brevity, we will refer to results using the full calculation as O(y 5 ) and leading order results as O(y 2 ).
Loops of the hard energy scale πT (i.e. the dimensional reduction) yield an expansion in integer powers of y 2 .The half-integer powers come from loops within the soft scale EFT, where the energy scale is yT (or √ λT ).The expansion parameter for the soft scale is α 3 ∼ y; see equation (3.20).
To derive the effective couplings of the EFT, we matched the off-shell correlation functions at soft external momenta: zero Matsubara modes with p ∼ yT .We evaluated the connected, one-particle irreducible correlation functions, denoted by Γ (k) in the 4d theory and Γ (k) 3 in the 3d effective theory.The calculations yield generic soft-scale observables accurate up to the O(y 5 T 4 ) term for the free energy density.These are given in appendix B. By matching these physical observables, we arrive at expressions for the 3d effective parameters: (2.8) 2  (2.9) (2.10) (2.11) where couplings either do not have running at O(y 4 ), or are barred couplings defined to make transparent the renormalisation scale invariance at O(y 4 ): where β χ = β b χ + β f χ denotes the beta functions, given in appendix D, L b (Λ) and L f (Λ) are logarithms of the renormalisation scale, defined following ref.[29] and given in appendix B, γ E is the Euler-Mascheroni constant, and c is a constant arising from 2-loop sum-integrals, also given in appendix B. Similarly to the barred couplings, we defined: where γ ϕ = − 1 φ 0 dφ 0 /d(log Λ) = −2y 2 /(4π) 2 .All couplings and fields here are MS renormalised.These dimensional reduction relations extend the leading-order results of ref. [60].
These expressions for the 3d effective parameters now enable perturbative analysis of the EFT to yield results applicable to the 4d Yukawa model.The scale Λ can be replaced with two separate renormalisation scales, µ and µ 3 which may be chosen independently [29].The Λ scales explicitly shown in equations (2.9) and (2.10) reproduce the renormalisation group running of the 3d EFT, allowing us to exchange for some new scale µ 3 that can be chosen independently.Given the superrenormalisability of the 3d EFT, this running can be upgraded to be exact by replacing gλT3/2 → g 3 λ 3 and λ 2 T 2 → λ 2 3 in front of the logarithms of Λ in equations (2.9) and (2.10).The remaining Λ scale dependence carried implicitly in the 4d MS parameters is then denoted as µ for clarity.
3 Phase transitions and latent heat

Conditions for first-order phase transition
The pattern of phase transitions in this 3d EFT has been studied in ref. [61].The existence and order of any phase transitions can be determined as follows.One can find a basis in which g 3 (T ) = 0 for all T by shifting ϕ 3 → ϕ 3 − g 3 /λ 3 .The bare potential then takes the following form: where we have defined The mass counterterm is given in equation B.17, and the possible tadpole counterterm cancels.This reduces the theory to the Z 2 -symmetric ϕ 4  3 theory coupled to a finite Z 2 -breaking external field σ3 (T ).A phase transition, if it exists, occurs at where the symmetry is restored.This defines the critical temperature, T c .As a consequence of equation (3.3), the properties of the phase transition, such as its order and strength, may depend only on m2 3 and λ 3 . 2Further, by dimensional analysis, there can only be non-trivial dependence on the dimensionless ratio m2 3 /λ 2 3 .Using mean field theory, one can see that for m2 3 /λ 2 3 ≪ −1 there is a first-order phase transition, and for m2 3 /λ 2 3 ≫ 1 there is a crossover.In between, the line of first-order phase transitions ends in a second-order phase transition at m2 3 = m2 3, * .The value of m2 3, * has been determined using lattice Monte-Carlo simulations to be [61,62] m2 3, * = [0.0015249(48)+ with the number in parenthesis being the statistical uncertainty in the last digits, and the whole expression is evaluated at the critical temperature.
Due to the Z 2 symmetry at the critical temperature, ϕ 3 → −ϕ 3 , the condition σ3 = 0 is exact, and is unmodified by loop corrections within the 3d EFT which can yield odd powers of y.Therefore, the perturbative expansion of T c is one in powers of y 2 , all odd powers of y cancelling due to the symmetry, The power counting of other couplings and parameters are related to y using the prescription of equation (2.4).The LO result T (0) c can be computed with an unresummed 1-loop calculation, while the next-to-leading order correction a 2 requires 2-loop dimensional reduction to be carried out with the free energy evaluated to at least O(y 4 T 4 ).The coefficients a n can be obtained analytically by solving equation (3.3) in a strict expansion in powers of y.
The order of the transition is determined by the sign of m2 3 − m2 3, * at the critical temperature.m2 3 is determined by purely hard-scale physics and m2 3, * is determined by purely soft-scale physics.However, the soft physics is known from lattice simulations, so the only remaining uncertainties in m2 3 − m2 3, * come from the hard scale.

Phase transitions for Yukawa theory
For a given parameter point in our Yukawa model, the critical temperature (if it exists) and the order of any phase transition can be determined by combining equations (3.3) and (3.4) with the expressions for the 3d effective parameters, equations (2.8) to (2.12).Using the dimensional reduction relations at leading order, equation (3.3) for the critical temperature has a simple analytic solution, about which subleading corrections can be obtained by a strict perturbative expansion.Our higher order results include the subleading term suppressed relatively by y 2 ; see equation (3.5).In what follows, we give explicit analytic results for the critical temperature to this order.
To simplify the results, we introduce 4d tilded parameters, defined in analogy with those in 3d by transforming ϕ → ϕ−g/λ to eliminate the cubic scalar term.In terms of the (untilded) parameters of the original basis, the (tilded) parameters of the new basis read and y and λ are unchanged.
The coefficient of the symmetry-breaking term in the 3d EFT then reads Solving σ3 (T c ) = 0 in a strict expansion, using the free energy to O(y 5 T 4 ), we find the critical temperature, At leading order, this can be expressed in terms of the original 4d parameters as For y = 0, the scalar and fermion decouple, and √ T σ3 becomes independent of T for any choice µ ∝ T , so σ3 cannot change sign for the scalar-only theory, at least up to this order.As a consequence the only possible phase transition for the scalar-only theory is of second order, where σ3 = 0 for all temperatures and the transition takes place when m2 3 − m2 3, * goes through zero [62].Similarly, σ3 cannot change sign in the Z 2 limit of this model, where σ = mψ = 0, but y ≠ 0.
If T 2 c is negative, then there is no phase transition.If it is positive, equations (3.2) and (3.4) can be used to determine m2 3 − m2 3, * and hence the order of the phase transition.The 3d effective mass parameter has a perturbative expansion in integer powers of y 2 , which, when evaluated at the critical temperature, reads where we have chosen both renormalisation scales proportional to the temperature, with µ = AπT , and µ 3 = BT .Renormalisation scale dependence can be used to estimate the size of missing higher order terms, and so to give an intrinsic measure of the theoretical uncertainty in a prediction.In a strict expansion of a physical quantity, any renormalisation scale dependence cancels exactly, order by order.However, by solving for the running of the parameters to a higher order than the underlying calculation, a formally higher order renormalisation scale dependence can be introduced.
In practice, we have solved the 1-loop renormalisation group equations perturbatively, working to one higher order than the thermodynamic calculation; see appendix D for the beta functions.That is, denoting a 4d MS parameter by κ, we have solved for the running up to and including the O(κy 2 ) term for our LO thermodynamics, and up to the O(κy 4 ) term for our higher order thermodynamics.As before, y here is used as a shorthand to include terms of equivalent sizes according to the power counting scheme in equation (2.4).Writing a given MS parameter κ as a function of t = log (µ ′ /µ), the higher order perturbative solution to the running is where κ a runs over all the parameters in the model.Note that while 2-loop contributions to the beta functions would contribute at O(κy 4 ), we do not include them, as we are merely using the running to estimate uncertainties.These running couplings were then included in expressions for the thermodynamic quantities, such as eqs.(3.8) and (3.10).By then varying the renormalisation scale over a range, we can obtain a measure of the intrinsic uncertainty of the predictions.We chose the renormalisation scales as µ = AπT c , where A varies over 7 equidistant values from 1/2 to 2. This range follows standard practice in the field [33,63].The arithmetic mean of the results is taken as the predicted physical quantity, and the minimum and maximum values yield the range of uncertainty.For simplicity, we fix B = 1 throughout.

Scan of parameter space
A global study of the parameter space was carried out.Without loss of generality, units were chosen such that m 2 = 1, and σ was set to zero (this can be fixed by a shift ϕ → ϕ + constant).The distributions of the other parameters y 2 , g 2 , λ and m ψ were chosen according to where U(a, b) is the uniform distribution on the interval [a, b), and the signs of y and g were chosen randomly and with equal probability.Perturbativity at zero temperature requires that all three couplings are small compared with ∼ 8π (see e.g.[64]), hence the maximum of the range, and a log distribution was chosen to explore a wide range of possible ratios of couplings.
A linear distribution was chosen for the fermion mass parameter because the boundaries between different behaviours occur at O(1) values of m ψ , and no special behaviours are expected (or observed) in the m ψ → 0 limit.Further, large values of m ψ are more likely to invalidate our high-temperature approximation.
If treated as physical parameters, to be consistent with the accuracy of the rest of our calculations, the values generated should be matched to MS parameters at zero temperature at 1-loop order [29].However, as we are primarily interested in the perturbative convergence of thermodynamic predictions, rather than interpreting the predictions in a cosmological context, we neglect this zero-temperature matching between MS and physical parameters, instead treating this as a toy model.We therefore carry out a different scan, directly over MS parameter space at the input renormalisation scale µ 0 = 1, carrying over their values directly to the phase transition evaluation.
Fig. 1 shows the phase transition evaluation of the Yukawa theory, using the full O(y 5 ) calculation as per equation (2.7), for a scan of 20,000 parameter points using the distributions summarised above.The axes are chosen as combinations of parameters which bound the first-order phase transitions identified.We find that first-order phase transitions are limited to parameter space where g 2 λ ≳ 2m 2 .At leading order, this condition follows from the requirement m2 3 − m2 3, * < 0 for first-order phase transitions, which can be written as Figure 1: Predicted nature of scalar singlet Yukawa theories studied with free energy evaluated to O(y 5 T 4 ) for our scan of parameters.First-order phase transitions were found to be restricted to parameter space where g 2 λm 2 ≳ 2 and y 4 λ ≲ 0.5.There are two additional directions in parameter space not shown here, λ and m ψ /m. with the second term positive definite.
Parameter points are also shown which have no leading order solution for the critical temperature (in orange), or are predicted to be crossovers (in black).Following this are parameter points, shown in purple, where the high temperature approximation does not hold, defined as T c < 2 ⋅ max(m, m ψ ).This precise condition is motivated by the convergence of the high-temperature expansion for the 1-loop thermal functions [57].
It can also be seen from fig. 1 that for values of y 4 λ ≳ 0.5, our power counting assumptions can begin to break down.For the points in green, the correction to the critical temperature from formally subleading terms in the O(y 5 T 4 ) free energy density and from running is larger than the leading order value.We have denoted as non-perturbative such points which satisfy for some renormalisation scale µ in the range considered, where T c (µ) is the critical temperature computed with the O(y 5 T 4 ) free energy density and O(y 6 ) running, and T (0) c is the critical temperature evaluated at O(y 2 ) with no running.At leading order, an analogous non-perturbative condition can be defined, with T c (µ) now the critical temperature computed to O(y 2 ) with running at O(y 4 ).
For the parameter points in cyan, λ is negative after running, which can be interpreted a breakdown of the model.Finally, grey points have negative m 2 after running, which indicates large corrections from perturbative renormalisation group running.
The theories identified as first-order phase transitions were required to be consistently evaluated as such across the range of renormalisation scales µ = AπT c , where A varies from 1/2 to 2. This is further discussed in section 4.3.
As there are only four independent parameters, it has been feasible to attain a reasonably dense scan of the relevant parameter space with a simple random scan, and to illustrate features of the results with projections onto a plane in parameter space.While the distribution of our results inevitably depends on the details of our scan choice, all our calculations were carried out on the same set of parameter points, allowing us to reveal differences in predictions based on the O(y 2 ) and O(y 5 ) approximations.In theories with higher dimensional parameter spaces, we would expect adaptive search algorithms to become necessary for efficient study of the phase transitions of a model [65].

Latent heat evaluation
After determining the type of phase transition at a parameter point, if any, attention is next turned to the strength of the phase transition in terms of the latent heat.Strong phase transitions are of particular interest for detectable gravitational wave signals.
The latent heat, L, of a first-order phase transition can be determined by the following thermodynamic relation evaluated at the critical temperature: where f is the free energy density of the full 4d theory and ∆f ≡ f + − f − , where f + and f − are the free energy densities in the higher and lower temperature phases respectively.This can be expressed in terms of the effective parameters of the 3d EFT as: where σ3 is given in equation (3.1) and ⟨ φ3 ⟩ is the linear field condensate.The definition of the former is structurally akin to a beta function.It depends only on the UV physics of the hard scale, and can be read off from equation (3.7).The field condensate depends on the IR physics of the soft-scale EFT.At leading order, it is related to: where is the field value of the tree level minima about the Z 2 -symmetric origin.
The jump in the linear condensate has been evaluated to three loops within the EFT, and is given by [61]: where µ 3 is the 3d EFT MS renormalisation scale, and we have introduced the loop-expansion parameter α 3 : Note that λ 3 ∼ y 2 T , while v 2 0 ∼ T so that α 3 ∼ y.L/T 4 c is a dimensionless measure of the latent heat, and therefore of the strength of the first-order phase transition.Similarly to the critical temperature, we evaluated this using free energy evaluated to both O(y 2 T 4 ) and O(y 5 T 4 ).Parameters were run to higher orders, and varied over a range of renormalisation scales, to measure the uncertainty in the predictions.
Fig. 2 plots L/T 4 c for consistent first-order phase transitions, with a maximum value of ∼ 0.5 found in our scan at O(y 5 ).The latent heat is plotted against y 2 m/g, with a strong positive correlation between the two.Since theories with y 2 m/g ≳ 0.5 begin to be non-perturbative in our analysis (Fig. 1), a follow up investigation could identify additional, stronger, first-order phase transitions at larger values, for which a different power counting may be necessary, as y 2 m/g ∼ y is no longer small., is plotted for all the first-order phase transitions in the parameter space studied within our scan, evaluated at O(y 5 ).It is found to be positively correlated with the ratio y 2 m/g.

Comparing results to leading order
The calculations in this work have been carried out to three powers of y beyond leading order, as set out in equations (2.4) and (2.7).This section will investigate the size and properties of the sub-leading corrections in the context of the properties of first-order phase transitions.
A number of works have found large discrepancies between calculations of phase transition parameters at low perturbative orders.For example, in the Z 2 -symmetric limit of the realscalar-extended SM, refs.[37,38] have found huge discrepancies between computations of the effective potential at order O(y 2 T 4 ) and O(y 3 T 4 ).Ref. [33] found significantly reduced uncertainties in this model at order O(y 4 T 4 ).Around the critical temperature, properties of the phase transitions are highly sensitive to additional sub-leading terms, contributing to large corrections to the properties of expected gravitational waves.

Additional strong phase transitions
Across different renormalisation scales, there may not be only quantitative differences in the values of properties such as latent heat, but also qualitative differences in the predicted nature of the phase transition.This can mean for example that a particularly theory has a firstorder phase transition at one renormalisation scale, but becomes a crossover when evaluated at another scale.
This has been treated in ref. [66] through a characterisation from "ultra-conservative" to "liberal", based on the consistency of results with respect to scale.Here we adopt an approach that would be considered "conservative" in this sense, where a theory associated with a set ) and not at O(y 2 ), using power counting as per eqs.(2.4) and (2.7).For consistent first-order phase transitions at both O(y 2 ) and O(y 5 ), the changes due to loop corrections for the latent heat and critical temperatures are plotted as red arrows.There are no parameter points that are consistent first-order phase transitions at leading order but not at O(y 5 ).
of parameters is considered to be a consistent first-order phase transition if it is true across the range of renormalisation scales µ = AπT c , where A varies over 7 equidistant values from 1/2 to 2. The theories shown in blue in fig. 1 are all consistent first-order phase transitions in this sense, whereas points in other categorisations were evaluated to have the corresponding prediction (e.g.crossover) at one or more of the renormalisation scales considered.
It can be seen from figure 3 that a large number of parameter points are only found to be consistent first-order phase transitions with the higher order O(y 5 ) calculation, and not when evaluated at leading order, O(y 2 ).The higher loop order calculation significantly changes the landscape of strong first-order phase transitions.The additional first-order phase transitions at O(y 5 ) in our scan correspond to a 57% increase in FOPTs with L/T 4 c > 0.1 (from 229 to 359), or a 270% increase in FOPTs with L/T 4 c > 0.2 (from 19 to 71). Figure 4 plots the same points, now also indicating the nature of the theories which are consistent first-order phase transitions at O(y 5 ), but not at O(y 2 ), for at least one renormalisation scale.Over half of these are due to non-perturbative corrections when evaluated at leading order, with running at O(y 4 ), as defined in section 3.3.
Most of the remainder are predicted to be crossovers ( m2 3 − m2 3, * > 0) for at least some part The leading order calculation failed to identify the majority of the strongest firstorder phase transitions found at O(y 5 ).At one or more renormalisation scales, many of these were found to be not perturbative at leading order -where contributions from the next order are larger than the leading order value of T c .Other strong phase transitions were categorised as crossovers at leading order. of the renormalisation scale range at leading order.Two example parameter points are shown in figure 5.In both cases, the O(y 2 ) evaluation indicates either a crossover ( m2 3 − m2 3, * > 0) or a first-order phase transition ( m2 3 − m2 3, * < 0), depending on the renormalisation scale used.Therefore, neither point is considered to be a consistent first-order phase transition at leading order.Using the O(y 5 ) calculation, the parameter point in figure 5a is confirmed to be a first-order phase transition, whereas the example in figure 5b is found to be consistently a crossover across the renormalisation scales considered.If their natures were determined using only the O(y 2 ) calculation at a single renormalisation scale µ = πT , the parameter point in figure 5a would be determined as a crossover, and the point in figure 5b a first-order phase transition, both opposites of the conclusions of the O(y 5 ) calculation.positive, for phase transitions at lower critical temperature, and negative for those at higher critical temperature.The critical temperature corrections are generally negative.This pattern of flow means that a given set of theories occupying some parameter space, is transported as a whole in the latent heat and critical temperature plane, as a result of including higher order contributions.In other words, the blob moves, at least in the context of our scan.

Movement of phase transitions
As well as in the space of the physics predictions of the theory, it is also possible to study the transport of predictions in parameter space when higher order contributions are added.Figure 6 plots theories predicted to be first-order phase transitions at O(y 2 ) and O(y 5 ),  ), consistent first-order phase transitions expand in our scan's theory parameter space to include points with larger Yukawa couplings relative to λ.
with power counting as before per eqs.(2.4) and (2.7), against combinations of parameters y 4 /λ and g 2 /(λm 2 ).These are the same combinations as fig. 1, where they are found to be correlated with the type of phase transition predicted.The figure shows that, for our scan, within the higher order calculation the set of consistent FOPTs extends in the direction of larger y 4 /λ, shown in blue.Therefore, in parameter space, the blob grows, but does not move, in the context of our scan.
At larger still values y 4 /λ ∼ π 2 , 1-loop contributions from the Yukawa interaction compete with the tree-level scalar self-coupling in the zero-temperature effective potential.In this case a new set of power counting relations, and consequently resummations, are required for the perturbative description [36,67,68].

Reduced errors
Figs. 7 and 8 show relative errors for theories that are consistent first-order phase transitions, when evaluated at both O(y 2 ) and O(y 5 ) in the counting of eq.(2.7).The uncertainties were determined by varying the renormalisation scale between µ = πT /2 and µ = 2πT .Fig. 7 shows that nearly all of the theory space benefits from reduced uncertainty at higher order.It is instructive to focus on the strongest first-order phase transitions, as only these are expected to yield gravitational wave signatures which are observable by next-generation detectors [69,70].This can be seen in fig.8.
The errors are in general larger for the strongest transitions.For first-order phase transitions in our scan that are both strong (with L/T 4 > 0.1) and consistent in both calculations, There is generally a reduction in error in predictions of the latent heat at O(y 5 ) compared to O(y 2 ), in our scan, for theories that are consistent first-order phase transitions at both orders.the higher order calculation reduces the average relative range in critical temperature from 18% to 4%, in latent heat from 19% to 6% and in m2 3 − m2 3, * from 19% to 7%.If we considered even stronger phase transitions, with L/T 4 > 0.2, the higher order calculation reduces the average relative range in critical temperature from 28% to 7%, in latent heat from 40% to 16% and in m2 3 − m2 3, * from 13% to 5%. Figure 9 provides an example of a strong first order phase transition.The uncertainties in the critical temperature and the strength of the phase transition were significantly reduced at O(y 5 ) compared to the leading order O(y 2 ) evaluation.Note that below the scale choice A ∼ 0.4, these quantities are not evaluated at O(y 2 ) as the parameter point is found to receive non-perturbative corrections to the critical temperature.Since this is outside of the range 1/2 ≤ A ≤ 2, this parameter point is still considered a consistent first-order phase transition at O(y 2 ).
As well as a reduction in uncertainty, it can be seen from figure 8 that a small number of theories have larger estimated errors when evaluated for the higher order calculation.This reveals an underestimate of the uncertainty associated with the perturbative calculation when evaluated at leading order, and is another improvement offered by the calculation computing free energy density up to O(y 5 T 4 ).Interestingly, the lower half of the O(y 2 ) plot resembles  c are plotted with their relative errors for theories that are consistent FOPTs when evaluated with f at both O(y 2 T 4 ) and O(y 5 T 4 ).There is a large reduction in relative error at higher order, especially for strong phase transitions, for our scan.There are also a number of points where the LO result underestimated the error.Figure 9: Uncertainties in T c and L/T 4 c are significantly reduced when evaluating using the free energy density to O(y 5 T 4 ) here shown for a benchmark point with a strong transition.
the O(y 5 ) plot, when the former is scaled up by a factor of 5 ∼ 10, i.e. an order of magnitude.
A "conservative" approach was adopted here, where a first-order phase transitions is required to be consistently identified across a range of renormalisation scales.A similar story emerges if a more "liberal" approach is taken, such as evaluating the nature of the theory only at a single renormalisation scale, for example setting A = 1.Another option is to include a theory as a first-order phase transition if it is evaluated as such for any subset of the renormalisation scale range considered.Using each of these alternative approaches, applied to our scan, the calculation at O(y 5 ) evaluates additional strong phase transitions not identified at leading order using the same approach, combined with reduced relative error in the critical temperature and latent heat evaluated.

Discussion
Let us outline an answer to the central question of our study of cosmological phase transitions: how important are higher perturbative orders for a global parameter scan of a model?Previous studies have shown that higher orders are vital for making reliable predictions for cosmological phase transitions at individual benchmark parameter points in a model [31,33,[50][51][52]54].Yet it was unclear if these conclusions would extend to a global parameter scan, or whether higher perturbative orders would merely reshuffle parameter points.
We have tackled these questions using a Yukawa model as our guide.The simplicity of this model has allowed us to push to O(y 5 T 4 ) accuracy for the free-energy density.To do so, we have carried out high-temperature dimensional reduction, constructing a 3d effective field theory of only the lightest bosonic modes, and have included relevant infrared contributions up to three loops.
In short, higher orders reduce errors across the board, but are most important for describing the strongest phase transitions, where errors are largest.At leading order, the renormalisation scale dependence of the strongest phase transitions we have identified yield error estimates over 100%, or change the character of the transition from first order through second order to crossover.Only at higher loops can the order of these strongest transitions, and their thermodynamic properties, be computed reliably.This is expected to significantly impact predictions of the resulting gravitational wave signals, as these depend sensitively on the thermodynamics of the transition.Yet, additional nonequilibrium physics enters the computation of the gravitational wave spectrum, through the bubble nucleation rate and bubble wall speed.The calculation of higher orders for these quantities is in its infancy: for the nucleation rate the O(y 3 ) term has recently moved into reach [60,71,72], while for the bubble wall speed even a complete leading order calculation would extend the current state-of-the-art, which relies on a leading-log expansion [73].
Future work could extend the model to include couplings to the Standard Model through the Higgs portal.Carrying through to explicitly evaluating the consequences for gravitational waves would then yield results directly applicable to the widely studied singlet fermionic dark matter candidate, extending the conclusions of previous works [21][22][23] to higher orders.While the question remains as to how to compute higher orders for the nonequilibrium quantities which enter the gravitational wave spectrum, already some relevant information can be gleaned by upgrading the equilibrium calculation.
A natural question to ask next is: to what extent do we expect our conclusions to apply to other models?Crucially, strong first-order phase transitions are associated with large couplings and large errors also in scalar extensions of the Standard Model [33,51,52,54,74], as well as in the Standard Model Effective Field Theory [31].Therefore, it seems plausible that the broad brushstrokes of our conclusions may apply to many other models, at least when the structure of the power counting relations is similar, that is, when the transitioning field is of the soft scale m eff ∼ yT .However, there may be qualitative differences for symmetry-breaking transitions, where the transitioning field instead lives at the supersoft scale m eff ∼ y 3/2 T /π [59,60,75,76], or for very strong phase transitions where additional hierarchies of scale occur [36,54].

A Power counting prescription
At high temperatures, thermal fluctuations (loops) qualitatively change the structure of the effective potential.At such temperatures, the vanilla loop expansion breaks down but one can still make progress by counting powers of couplings.To do so in a theory with multiple couplings, it is useful to fix power counting relations between the couplings.
To start with we assume that all three interaction terms are approximately equally perturbative.This assumption anyway underlies the usual loop expansion at zero temperature, and amounts to parametrically equating the loop expansion parameters for each interaction The latter two terms are squared because each additional loop requires two cubic interactions, but only one quartic interaction, and follows from standard graph theoretic results [77].The inverse powers of m in the last term are necessary to cancel the dimensions of g, and arise from loop integrals of ϕ.Note that we have not made explicit the factors of 1/(4π) arising from loop integrations.
In addition, we will assume that the leading effects of high-temperature fluctuations are of the same size as the tree-level scalar potential.For the scalar tadpole and mass, this amounts to This assumption is natural in the vicinity of a phase transition, as thermal fluctuations must balance against tree-level terms.
Finally, we must make an assumption about the fermion mass, m ψ .For the fermion to have any effect on a phase transition, its mass cannot be large compared with the temperature.Otherwise its effects would be exponentially (Boltzmann) suppressed.For simplicity, we will power count the fermion mass the same as the scalar mass, m ψ ∼ m .
(A.3) Our analysis will nevertheless include m ψ ≪ m as a special case.In summary, we adopt the following power counting prescription: A generic equilibrium observable θ has an expansion of the form θ/θ LO = 1 + c 1 y + c 2 y 2 + c 3 y 3 + . . ., and in this work we calculate up to and including the c 3 y 3 term.For the effective potential and latent heat, which are O(y 2 T 4 ) at leading order, this means evaluation to O(y 5 T 4 ), or equivalently O(λ 5/2 T 4 ), with errors of O(y 6 T 4 ).For the critical temperature, O(my −1 ) at leading order, the evaluation is to O(my) as there is in general no contribution at O(my 2 ), with errors of O(my 3 ).

B Evaluating and matching correlation functions
We evaluated the connected, one-particle irreducible correlation functions, denoted by Γ (k)  in the 4d theory and Γ (k) 3 in the 3d effective theory.These observables were then be used to derive matching relations between the theories to find expressions for the effective parameters.
In constructing the EFT [78,79], we are aiming to capture the averaged effect of the hard scale modes (∼ πT ) on the soft scale modes (∼ yT ).As such, in the computation of correlation functions we are free to cut off the soft modes in any way, as long as we do so on both sides of the matching relations.This observation underlies the strict perturbation theory approach of Braaten and Nieto [28], which we follow.Essentially one expands all loop integrals in powers of the soft scale, and then uses dimensional regularisation to cut off any infrared divergences.In this context, matching correlation functions is equivalent to directly integrating out the hard modes [75,79].
The one, two, three and four-point correlation functions are given by the sum of diagrams in figs.10, 11, 12 and 13 respectively.These were generated using FeynArts [80].Only diagrams and terms up to order O(y 5 ) in the power counting prescription eq.(2.4) were included in the calculations.Tadpole and mass terms are treated as interactions, according to strict perturbation theory, due to the hierarchy of energy scales between their coefficients and πT which characterises the free Lagrangian.The results of the loop sum integrals can be found in the literature and are listed in appendix C.
The evaluation of the correlation functions are given in equations (B.4) to (B.11).The Euler-Mascheroni constant is denoted by γ E , and A is the Glaisher-Kinkelin constant.Zerotemperature counterterms, given in appendix D, are used in the evaluation to cancel the temperature independent divergences.In order to simplify the formulae, we have also introduced notation following [29]: Figure 10: The Feynman diagrams contributing to the one-point correlation function, Γ (1) (0) up to O(y 4 T 3 ), with power counting set out in eq.(2.4).Γ (3)  (p, q, −p, −q) ≈ g + δg − (p, q, r, −p, −q, −r) up to O(y 4 ).
Γ (4)  (p, q, r, −p, −q, −r) ≈ λ + δλ − In evaluating these equations we expanded assuming p ∼ yT .Performing the same calculation in the effective theory, and using the vanishing of scaleless integrals in dimensional regularisation, we find trivially that, Γ Following ref. [29], we match the leading momentum dependence of the correlation functions through a field normalisation, as well as equate the 1 to 4 point correlation functions at zero momentum between the full and effective theories.Using these relationships we arrive at the expressions (2.8) to (2.12) for the 3d EFT parameters.The following counterterms are required in the 3d theory to cancel against the temperature dependent poles from the full theory, where thermal four momenta are denoted by uppercase letters, P = (p 0 , p), with spatial momenta p and Matsubara frequencies p 0 = 2πT n for bosonic sum integrals (denoted by P ) and p 0 = (2n + 1)πT for fermionic sum integrals (denoted by {P }).Massless 1-loop bosonic integrals can be evaluated using the following master sumintegral, (D.10) The beta functions for the dimensionless couplings are corroborated by refs.[11,82].Barred couplings explicitly demonstrating renormalisation scale invariance for the 3d EFT parameters in eqs.(2.8) to (2.12) are defined as: where for each parameter, the full beta function is equal to the sum of the bosonic and fermionic parts β = β b + β f , with: (D.17) For m ψ and y there is some arbitrariness in this split, which has been chosen in order to minimise the presence of log(2) terms in the dimensional reduction relations, equations (2.8) to (2.12).

Figure 2 :
Figure 2:The strength of the phase transition, represented here by the dimensionless L/T 4 c , is plotted for all the first-order phase transitions in the parameter space studied within our scan, evaluated at O(y 5 ).It is found to be positively correlated with the ratio y 2 m/g.

1 •Figure 3 :
Figure3: The blue points show theories which are only consistent first-order phase transitions (i.e.valid across renormalisation scales) at O(y 5 ) and not at O(y 2 ), using power counting as per eqs.(2.4) and (2.7).For consistent first-order phase transitions at both

Figure 4 :
Figure 4:The leading order calculation failed to identify the majority of the strongest firstorder phase transitions found at O(y 5 ).At one or more renormalisation scales, many of these were found to be not perturbative at leading order -where contributions from the next order are larger than the leading order value of T c .Other strong phase transitions were categorised as crossovers at leading order.

Figure 3 Figure 5 :
Figure3also shows the change in latent heat and critical temperature for all consistent first-order phase transitions identified at both leading order and O(y 5 ).The red arrows visualise distinctive flows of the higher order contributions in the latent heat and critical temperature plane for our parameter space scan.The latent heat corrections are generally

Figure 6 :
Figure 6: From O(y 2 ) to O(y 5), consistent first-order phase transitions expand in our scan's theory parameter space to include points with larger Yukawa couplings relative to λ.

Figure 7 :
Figure7: There is generally a reduction in error in predictions of the latent heat at O(y 5 ) compared to O(y 2 ), in our scan, for theories that are consistent first-order phase transitions at both orders.

Figure 8 :
Figure 8: T c /m and L/T 4c are plotted with their relative errors for theories that are consistent FOPTs when evaluated with f at both O(y 2 T 4 ) and O(y 5 T 4 ).There is a large reduction in relative error at higher order, especially for strong phase transitions, for our scan.There are also a number of points where the LO result underestimated the error.