Real scalar phase transitions: a nonperturbative analysis

We study the thermal phase transitions of a generic real scalar field, without a $Z_2$-symmetry, referred to variously as an inert, sterile or singlet scalar, or $\phi^3+\phi^4$ theory. Such a scalar field arises in a wide range of models, including as the inflaton, or as a portal to the dark sector. At high temperatures, we perform dimensional reduction, matching to an effective theory in three dimensions, which we then study both perturbatively and on the lattice. For strong first-order transitions, with large tree-level cubic couplings, our lattice Monte-Carlo simulations agree with perturbation theory within error. However, as the size of the cubic coupling decreases, relative to the quartic coupling, perturbation theory becomes less and less reliable, breaking down completely in the approach to the $Z_2$-symmetric limit, in which the transition is of second order. Notwithstanding, the renormalisation group is shown to significantly extend the validity of perturbation theory. Throughout, our calculations are made as explicit as possible so that this article may serve as a guide for similar calculations in other theories.


Introduction
Since the Hot Big Bang, the universe may have passed through a number of different phases. In the Standard Model (SM) of particle physics, electroweak symmetry breaking and colour confinement took place at temperatures of approximately ∼160 GeV [1] and ∼155 MeV [2][3][4] respectively, though both these transitions are smooth crossovers. Extensions of the Standard Model may lead to a wide variety of phases and phase transitions in the early universe. Such phase transitions may have an importance for baryogenesis [5][6][7][8], and may lead to a detectable signal of gravitational waves [9], allowing the possibility to probe particle physics in a completely new way. The gravitational waves produced by first-order phase transitions even offer the possibility of studying dark sectors that are uncoupled to the Standard Model [10][11][12]. Scalar fields, whether fundamental or effective, often lie at the heart of phase transitions, acting as order parameters which take different expectation values in different phases. In this article, we study the simplest scalar theory which gives rise to a first-order phase transition. Its Lagrangian is given by, where φ is the scalar field, µ runs over 0,1,2,3 and we have used the mostly minus metric signature. It has been referred to variously as an inert, sterile or singlet scalar, and in the most part we will refer to it as a singlet scalar. The singlet scalar may couple to the SM, or other fields, in which case the full Lagrangian takes the form where the portal sector contains SM-singlet interactions. The singlet extended SM is referred to as the xSM [13]. Indeed couplings to gravity are always present, couplings to the Higgs field are generically expected on effective field theory grounds, and likewise for couplings to sterile neutrinos [14] in minimal extensions of the Standard Model such as the νMSM [15][16][17]. On the other hand, being a single real scalar, there can be no (renormalisable) couplings to gauge fields, or to the charged chiral fermions of the Standard Model. In this article, we will consider phase transitions in the singlet direction, focusing on the case where the effects of all other fields can be accounted for by modifying the effective couplings of the infrared modes of the singlet.
We are motivated to study this model, in part, because it arises in a wide range of cosmological and particle-physics model building: providing a possible dark matter candidate [18][19][20][21], acting as the inflaton [17,[22][23][24][25] and providing for electroweak baryogenesis [26]. As a consequence, it has been the focus of many collider searches [13,[27][28][29][30]. An additional recent attraction to this work has been the possibility of observing a gravitational wave background from such a first-order phase transition [31][32][33][34][35][36][37][38][39][40][41] at future gravitational wave detectors, such as LISA [42], DECIGO [43], BBO [44] and Taiji [45]. Finally the simplicity of the real scalar theory is itself a motivation for its study, as it offers the possibility of carrying out relatively high order calculations explicitly, in turn allowing us to test the convergence and reliability of perturbation theory.
This last point, to test the convergence and reliability of perturbation theory, is another important motivation for this work. The gravitational wave spectrum produced by a first-order phase transition depends sensitively on the thermodynamics of the transition, which in turn is difficult to calculate reliably. A recent work [46] investigating this in a minimal extension of the Standard Model found that typical (one-loop) calculations suffer from huge multiplicative uncertainties in the gravitational wave peak amplitude, of order O(10 2 − 10 3 ) depending on the strength of the transition. In addition, Ref. [46] identified uncertainties of unknown, but potentially large, numerical importance. In light of all this, progress in the theoretical methodology is necessary to make robust predictions for future gravitational wave experiments.
Fundamentally, studies of high temperature physics are hampered by the strong coupling of light bosonic modes, which arises as a collective effect of their high occupancies. If at zero temperature λ is the perturbative expansion parameter, 1 at high temperatures, T , the effective expansion parameter for light modes, with mass m ≪ T , is modified as, (1.4) 1 Here λ stands for a generic expansion parameter, and not merely the coupling in this specific theory.
-2 -For such light modes, the effective coupling constant is much larger than the corresponding zero temperature coupling constant, and for sufficiently light modes the perturbative expansion breaks down altogether [47], as it does for non-Abelian gauge bosons in the symmetric phase.
For scalar fields at sufficiently high temperatures, thermal corrections to the effective mass grow quadratically and always dominate over the zero temperature mass. At such temperatures the lightest scalar modes have an effective mass m 2 ∼ λT 2 , implying that the effective coupling constant of such modes is reduced from λ → √ λ. By utilising effective theory techniques, such expansions in √ λ have been carried out to relatively high orders in several theories (see for example Refs. [48][49][50][51]).
However, in the vicinity of a phase transition, the situation is somewhat more difficult. Near the critical temperature, there is an approximate cancellation between the tree-level and thermal contributions to the mass of the field undergoing the transition, so that m 2 λT 2 and perhaps even m 2 ≪ λT 2 . This means that the effective coupling constant can be larger even that √ λ, and the perturbative expansion for these light modes can break down altogether. In this case, the only reliable method of calculation is lattice Monte-Carlo simulation.
Importantly, the potentially nonperturbative physics of the light bosonic modes is universal. Due to the hierarchy of scales, one can construct an effective field theory (EFT) for just the light modes [52][53][54]. This EFT is defined in 3d, and hence the construction is called (high temperature) dimensional reduction. The EFT depends on physics at shorter scales only through its effective parameters. Thus by studying the EFT with generic parameters, one arrives at results applicable to a wide range of 4d particle physics models. In this way, one needs only perform lattice Monte-Carlo simulations once, and the results can be recycled again and again; see for example Refs. [40,[55][56][57] where simulations of the SU(2)-Higgs 3d EFT were recycled.
In this article, we make use of dimensional reduction to perform a comprehensive study of phase transitions in which the infrared dynamics is governed by the singlet scalar. We perform dimensional reduction explicitly to next-to-leading order (NLO) starting from the 4d singlet scalar theory, and additionally we demonstrate the leading contributions due to couplings to other fields. In utilising dimensional reduction for carrying out the high temperature resummations, we differ from the widely adopted approach of (one-loop) daisy resummation [58,59] (see Ref. [60] for a recent review). We do so because dimensional reduction offers several advantages over daisy resummation (see Ref. [46] for a recent discussion), in addition to isolating the universal dynamics of the light bosonic modes.
Within the 3d EFT, on the one hand we perform a state-of-the-art perturbative calculation, computing the equilibrium properties of the low-energy effective theory to three-loop order, and utilising renormalisation group improvement. On the other hand, we put the theory on the lattice and calculate its properties at the critical temperature with relatively high statistics and utilising exact lattice-continuum relations. Our lattice Monte-Carlo simulations are carried out for parameter choices ranging over two orders of magnitude, allowing us to determine precisely where perturbation theory is reliable, and where it breaks down.
In outline, for a complete, nonperturbative calculation of the thermodynamics of a given particle-physics model, the following steps are performed: (I) Vacuum renormalisation: matching relations between physical observables and Lagrangian (MS) parameters [61].
(II) Dimensional reduction: matching relations between 4d theory parameters and effective 3d theory parameters [53,54,62] (III) Perturbative study in 3d effective theory: computation of thermodynamic quantities such as the jump of the order parameter and the latent heat [52,53,63] -3 -(IV) Framework for lattice Monte-Carlo simulations in 3d: lattice continuum relations [63][64][65] and possible O(a) improvement [66][67][68][69] (V) Lattice Monte-Carlo simulations of thermodynamic quantities in 3d: computations on fixed lattices, statistical analysis and taking the continuum limit [70] We perform all five steps in this article, as was also done in Ref. [71] for an MSSM-like model. One of our aims in this is simultaneously to give a bird's eye view of the technical steps required, and to flesh out all the relevant details, so as to offer a concrete guide to performing such calculations in other models. For relevant reviews, see Refs. [50,72,73].
Step (I) does not involve any thermal physics, and hence is rather separate to the others. As the calculation is standard, we present it in Appendix C. The following section, Sec. 2, consequently begins with Step (II), dimensional reduction. In Sec. 3 we analyse the phase diagram of the real singlet model at the broadest level: using only symmetries and results from the literature regarding the second-order phase transition in the Z 2 symmetric model. In Sec. 4 we carry out Step (III), a perturbative analysis of the thermodynamics of the phase transition. In Sec. 5 we perform Steps (IV) and (V), carrying out Monte-Carlo simulations for six different parameter points. Finally in Sec. 6 we discuss the implications of our results, in particular, we compare lattice and perturbation theory to better understand the limits of the validity of perturbation theory.
For completeness, we should note that in this article we do not study bubble nucleation or bubble growth, but content ourselves with studying purely equilibrium properties of the phase transition. Within dimensional reduction, a framework for studying bubble nucleation on the lattice has been developed in Refs. [68,74]. As argued in Ref. [46], following Refs. [75,76], dimensional reduction also offers a natural framework for semiclassical calculations of bubble nucleation. We plan a follow-up paper to this in which we study the bubble nucleation rate in this model.

Dimensional reduction
The equilibrium thermodynamics of quantum field theories can be studied in the imaginary time formalism (for a review, see Refs. [72,77]). In this case the fields live in R 3 × S 1 , i.e. three infinite spatial directions and one compact 'Euclidean time' direction, with length 1/T . In the compact direction bosons, such as the scalar φ, satisfy periodic boundary conditions and hence may be expanded in Fourier modes as, 2 (2.1) In this context the Fourier modes are referred to as Matsubara modes [78]. Here τ ∈ (−1/2T, 1/2T ] parameterises the compact direction, x ∈ R 3 parameterises the spatial directions and i runs over the spatial indices, 1, 2, 3. In terms of these modes, the quadratic part of the action for this scalar field takes the form, The rest of the scalar potential also provides interactions between Matsubara modes. In sum, one can view the equilibrium thermodynamics of a theory in 3+1 dimensions as the vacuum dynamics of a theory in 3 Euclidean dimensions containing infinitely many fields, ϕ n , with squared masses m 2 + (2πT n) 2 .
At high temperatures, when 2πT ≫ m, there is a hierarchy between the masses of the n = 0 and the n = 0 modes. This hierarchy causes some Feynman diagrams to become parametrically larger than their loop counting would suggest, hence requiring resummation. Dimensional reduction is a means to carry out such resummations in a systematically improvable way. It consists of integrating out the heavy (or ultraviolet (UV)) n = 0 modes, to derive an effective theory for the light (or infrared (IR)) n = 0 mode.
The practical steps of dimensional reduction were worked out independently in Refs. [52,54] and Ref. [53]. Here we follow the approach of Refs. [53,62] more closely, performing the matching in strict perturbation theory 3 for the generic real scalar theory. In outline the steps are: 1. Write down the most general 3d theory obeying the same internal and spatial symmetries and containing the same number of bosonic field degrees-of-freedom as the original 4d theory.
2. Calculate the static correlation functions for the operators in the Lagrangian in both the original theory and the effective theory. In both cases use strict perturbation theory.
3. Determine the coefficients of the 3d effective theory by matching the results of the two theories for momenta p ∼ √ λT .
In the following we will outline these steps, first for the pure, real singlet scalar model. Afterwards we will consider the effect of interactions with Standard Model, and other, fields.
Note that the approach we have adopted makes use of the high-temperature approximation in computing thermal loop integrals. While the validity of this follows naturally from the hierarchy of scales assumed in dimensional reduction, it is possible to avoid this approximation if necessary [79,80].

Minimal model
We start by considering the scenario where the real singlet scalar field is uncoupled to any other field, i.e. the Lagrangian is purely L singlet given in Eq. (1.1). This will allow us to set out the method, as well as to show many explicit details, given the simplicity of the model.
The accuracy of the matching procedure can be assessed by counting powers of coupling constants. The couplings of the 3d effective theory will consist of a sum of vacuum contributions and thermal contributions from the n = 0 modes which have been integrated out. We will restrict ourselves to sufficiently high temperatures such that the vacuum contributions are not parametrically larger than the thermal contributions. This amounts to the following power counting prescription, These conditions will generally be satisfied in the vicinity of the critical temperature, at which point thermal and vacuum contributions are balanced. There remains a final scaling relation between g, λ and T . For the theory to be perturbative at zero temperature it must be that g ≪ m, and hence we have that g ≪ √ λT . As a default assumption for the dimensional reduction, we will take g ∼ λT , however we will leave this final scaling relation more freedom than the others, because the order of the phase transition will depend on it.
We perform the matching at NLO, in which all effective couplings are calculated up to O(λ 2 ) multiplied by appropriate powers of T to make up the dimensions. This amounts to one-loop matching for the cubic and quartic couplings, and two-loop matching for the tadpole coupling and mass [54,81], though depending on the scaling relation between g and λ some Z 2 -breaking terms may be dropped as subdominant. Reaching this order is crucial for cancelling the leading renormalisation scale dependence.
For the pure real scalar theory, Step 1 is rather straightforward. The 3d effective theory is simply, We have used the subscript 3 to denote quantities pertaining to the 3d effective theory. Note this theory is a Euclidean field theory, hence the positive sign between derivative and potential terms. For Step 2, the philosophy is the following. The coupling constants of the low energy effective theory account for the effect of the UV modes of the full theory, which have been integrated out. Thus, in matching the two theories, we are free to treat the IR contributions in any way we choose, as long as we do so in the same way in both theories. That is because such IR contributions do not contribute to the coupling constants of the effective theory, and by treating the IR contributions in the same way in both theories, their contributions will cancel exactly. Given this freedom, the simplest way to treat the IR contributions is simply to cut them off in dimensional regularisation. Further, we are free to expand around any constant background field, including φ = 0, as the choice of background will only affect the IR. These choices significantly simplify practical calculations.
Strict perturbation theory in the full 4d theory is defined by the following split, where we have explicitly shown the counterterms, but not the powers of the scale Λ, which make up the dimensions in dimensional regularisation. Note that we adopt the MS renormalisation scheme, in which case there is no need for a field renormalisation counterterm, as it receives no divergent contributions. The omission of the terms proportional to σ and m 2 from the free Lagrangian is justified, within strict perturbation theory, by the smallness of these coefficients in comparison with the scale 2πT , which characterises the free Lagrangian. Strict perturbation theory in the low energy effective theory is defined by the following split, Due to the superrenormalisable nature of this theory in 3d, there are only a finite number of divergent diagrams. The counterterms needed are those shown here, which appear only at two-loop order. There are further divergent diagrams at four-loop order, but these are independent of the field configuration, contributing only to the cosmological constant, which we omit.
To derive the matching relations, we must choose a suitable set of observables to match between the full theory and the effective theory. Following Ref. [54], we choose these to be the connected, one-particle irreducible (1PI) correlation functions, denoted by Γ (k) in the full theory and Γ (k) 3 in the effective theory. These are equal to minus the sum of all connected, 1PI Feynman diagrams with k legs. We expand around zero background field, not the minimum of the potential, as the difference is anyway projected out by the IR cut off. The correlation functions are evaluated with soft external momenta: zero Matsubara modes and small spatial momenta, p ∼ √ λT . The one-point correlation function is given, up to two-loop order, by the sum of diagrams shown in Fig. 1, generated using FeynArts [82]. The results of the relevant loop integrals can be found Figure 1: Feynman diagrams contributing to the one-point function up to two-loop order. They are shown here in the same order that they appear in Eqs. (2.9) and (2.13).
in the literature and are listed in Appendix A. Evaluating the diagrams in the strict perturbative expansion, we find, where our notation for momenta and loop integration follows Refs. [53,62], and is given in Appendix A. We use the symbol ≈ to denote an equality which holds only in strict perturbation theory up to some loop order. In Eq. (2.10) γ E is the Euler-Mascheroni constant and A is the Glaisher-Kinkelin constant. In order to simplify the formulae, and following Ref. [54], we have introduced the notation, In going from Eq. (2.9) to Eq. (2.10) we have used the counterterms given in Appendix C, which cancel the temperature-independent divergences. Being temperature independent, these are the same counterterms that one would find in a calculation at T = 0.
Performing the same calculation in the effective theory, we find, (2.13) Again, the notation follows Refs. [53,62] and is defined in Appendix A. In going from Eq. (2.12) to (2.13) we have used that scaleless integrals vanish identically in dimensional regularisation. Note that the pure zero Matsubara parts of Eq. (2.9) give exactly the same scaleless integrals as in Eq. (2.13). Thus, one can see that the cancellation of the purely IR physics would also occur for other regularisation schemes. Possible contributions from IR-UV cross terms are essentially projected out by the strict perturbative expansion, which alternatively can be shown to cancel upon performing resummation [54].
Calculations of the remaining correlation functions, those with two, three and four external legs, are very similar and are given in Appendix B. The only real change is that for the two-point -7 -function, the leading O(p 2 ) momentum dependence must also be calculated. With the results for the correlation functions in hand, we can turn to Step 3 of dimensional reduction, matching.
For Step 3, we must match to find the field, φ 3 , as well as the four parameters of the effective theory, σ 3 , m 2 3 , g 3 and λ 3 . Naively, one might think that equating the four correlation functions would not give enough conditions, but in fact an extra condition arises from the leading momentum dependence of the two-point function. This extra condition essentially matches the kinetic operators, while the other four conditions match the potential terms.
To match the momentum dependence of the correlation functions, we must allow for a field normalisation factor between the zero Matsubara mode, ϕ 0 , and the low energy field operator, φ 3 [54]. Demanding that the p 2 part of the two-point functions agree, once this normalisation is taken into account, implies 14) where Π is the self-energy of the zero Matsubara modes in the full theory, as computed in strict perturbation theory, and on the second line we have introduced the dash to denote the derivative with respect to p 2 . The power of T arises so that φ 3 is canonically normalised in three dimensions.
Once the normalisation of the effective field operator is taken care of, matching the Lagrangian parameters simply amounts to equating the correlation functions at zero momentum and taking into account overall powers of T to make up the dimensions. Doing this leads to the following matching equations, From Eqs. (2.14)-(2.19), we can read off the effective couplings. Using the explicit expressions for the correlation functions in Appendix B, and expanding up to O(λ 2 ) and we find the effective parameters. Before writing the matching relations, we first make a couple of judicious modifications following Refs. [53,54]. We run the MS parameters of the 4d theory from the matching scale Λ, introduced by dimensional regularisation, to some new renormalisation scale µ. The beta functions are collected in Appendix C. This running may be essential to minimise large logarithms in perturbation theory, but can also be used to investigate any scale dependence of our result. After this, the tadpole and mass parameters retain some Λ-dependence, which will eventually be cancelled by two-loop diagrams in the 3d EFT. To extend this cancellation, which occurs at O(λ 2 ), to all-loop orders, we rewrite the coefficients of the Λ-dependent terms and the 3d counterterms in terms of the 3d effective parameters. This possibility is a consequence of the superrenormalisability of the 3d theory.
After these modifications, and assuming g ∼ λT , we arrive at our final result for the matching -8 -relations to the 3d EFT up to O(λ 2 ): where, in order to simplify the formulae, and following Ref. [54], we have introduced the constant To make clear the renormalisation scale dependence of the result, we have defined barred couplings κ, with κ ∈ {σ, m 2 , g, λ}, which are renormalisation scale invariant at one-loop order, where β κ ≡ dκ/d log µ denotes the one-loop beta functions, given in Appendix C. The parameters on the right hand side of Eqs. (2.20) to (2.26) have been run to the renormalisation scale µ, e.g. λ = λ(µ), whereas the 3d effective parameters are defined at Λ. There also remain temperature dependent 1/ǫ poles from the computation in the full theory, which must cancel against identical poles in the EFT. Matching these, we obtain The matching relations, Eqs. (2.21) to (2.28), pass several nontrivial checks. First, the counterterms, which have been derived from loop computations in the full 4d theory at finite temperature, are the just the right counterterms required to cancel UV divergences within the 3d EFT; see Sec. 4.2. Second, all dependence on the renormalisation scale Λ cancels up to the order that we have calculated. In particular, the matching relations for the field, the 3-point and the 4-point couplings are independent of the renormalisation scale Λ, For the tadpole and mass, σ 3 and m 2 3 , instead we find, Eqs. (2.29) to (2.32) are in fact the exact beta functions of the 3d EFT; see Sec. 4. As a consequence, the Λ-dependence arising from loop calculations within the 3d EFT cancels the Λ-dependence of -9 -the matching relations. This allows us to exchange Λ for some new scale µ 3 , chosen for example to minimise large logarithms within perturbative calculations in the 3d theory. Thus the original scale Λ (described as the matching scale in Ref. [53]) all but disappears, to be replaced by two renormalisation scales, µ and µ 3 which may be chosen independently.
We have also performed some checks against the literature. We find agreement with Refs. [52,53,59,83] (excluding what is presumably a sign error in Eq. (60) of Ref. [52]), which provide checks in the Z 2 -symmetric limit. Further, we find agreement with Ref. [81], which provides a check of all our correlation functions at one-loop order. Finally, we are grateful to the authors of Ref. [84], with whom we have cross-checked the full matching relations.

Higgs interactions
A singlet scalar may couple to the Standard Model Higgs, H, via the following portal couplings, In fact, generically such terms should be included as they are renormalisable and do not explicitly break any symmetries of the theory. Additionally, if the field φ is interpreted as the inflaton, there must be nonzero couplings to Standard Model particles in order to reheat the universe after inflation; see for example Ref. [85]. If the temperature T is above, or around, the electroweak symmetry breaking scale, then the zero Matsubara modes of the Higgs will also enter the 3d EFT. So too will the gauge bosons of the SM, though these do not couple directly to the singlet scalar. Further, such Higgs-portal couplings will give corrections to Eqs. (2.21) to (2.28) arising from the nonzero Matsubara modes of the Higgs, which at leading order amount to In calculating this correction, we have used the one-loop thermal correlation functions from Ref. [81]. However, a complete calculation of the dimensional reduction of the SM plus singlet is beyond the scope of this article; for which see the forthcoming Ref. [84]. If the singlet scalar appears around the electroweak scale, the dynamics of the coupled system will be a complicated interplay of the fields, perhaps involving a two-step transition [86][87][88][89]. This possibility has been studied in Refs. [31,[33][34][35][37][38][39][40]. While the full two-step transition goes beyond the scope of this article, our analysis is in principle applicable to a first step in the singlet scalar direction.
In Eqs. (2.34) and (2.35) we have assumed the temperature is above, or around, the electroweak symmetry breaking scale. Below this temperature, the Higgs takes a nonzero vev, modifying the effective couplings of the singlet at tree-level. As the temperature lowers further, perturbative excitations of the Higgs become exponentially suppressed, so that temperature-dependent Higgs corrections to the singlet EFT can be neglected, leaving only the tree-level effects of the Higgs vev.

Other possible interactions
Due to the lack of gauge charges, a real scalar field cannot couple to gauge fields or to the charged chiral fermions of the Standard Model. However, a variety of other interactions are possible.
Yukawa interactions with Dirac or Majorana fermions are possible, of the form, where A is a flavour index, y A are the Yukawa couplings and ψ A are the fermion fields. Such a Yukawa term may arise in models with sterile neutrinos [14], constructed to explain the observed small neutrino masses and to provide a dark matter candidate [17,21,[90][91][92][93][94][95][96][97][98][99][100]. In the νMSM, for example, the scalar field φ plays the role of the inflaton [17]. Such terms may also arise in simplified dark matter models involving a fermionic dark matter candidate and real scalar portal [101]. In the following, for simplicity, we consider the ψ A to be Dirac fermions. Fermions do not have zero Matsubara modes, due to the different boundary conditions implied by Fermi-Dirac statistics. As a consequence, they cannot enter the 3d effective theory. However, the Yukawa couplings will give corrections to the matching relations of the low-energy effective theory, which at LO amount simply to Here we have assumed that the masses of the fermions are small compared with the thermal mass scale πT , which is naturally the case if the y A are perturbative and the fermion masses are generated through this Yukawa mechanism. If, on the other hand, the fermions are much heavier than the thermal mass scale, their contributions to m 2 3 will be Boltzmann suppressed and can be neglected. In the intermediate case where the fermion masses are of the same order as the thermal mass scale, the full temperature dependence of the fermionic thermal functions must be retained.
Finally, we should mention that φ must inevitably couple to gravity, modifying the singlet Lagrangian as, where g is the metric determinant (not to be confused with the scalar cubic coupling), R is the Ricci scalar and ξ 1 and ξ 2 are coupling constants. Treating the metric as a classical, slowly-varying background field, we can see that the Ricci scalar corrects the tadpole and mass terms of the scalar φ already at tree-level, 4 and hence can potentially drive the scalar field through a phase transition. This has been studied in the context of a rapid quench, causing a spinodal, or tachyonic, instability for the scalar field in the Z 2 -symmetric version of this model [104,105]. On the other hand, a sufficiently slow time variation of the Ricci scalar can be incorporated in the effective parameters of the 3d EFT. Finally, we note that the full non-Z 2 -symmetric scalar field theory, with both Yukawa interactions and nonminimal couplings to gravity has recently been studied in Ref. [25], in which the scalar φ acts as the inflaton field. Under the assumption that the Z 2 -symmetry breaking terms are small, this model was shown to approximately reproduce the spectral index and tensor-to-scalar ratio of the Z 2 -symmetric Starobinsky inflation model.
In the following, we will focus on the thermodynamics of the 3d EFT containing only the real singlet scalar. We will remain agnostic about the full UV theory, and as a consequence will treat the 3d effective parameters, {σ 3 , m 2 3 , g 3 , λ 3 }, as unknown functions of temperature T .

Phase diagram
As the temperature changes, the couplings of the low-energy effective theory change accordingly. Cosmological history therefore traces out a path through the space of effective couplings, a path parameterised by the temperature. Observables, such as the free-energy density, or the expectation value of the field operator, will also change with temperature, following the variation of the effective couplings. In this way, the hard, nonzero Matsubara modes, which dictate the temperature dependence of the effective couplings, can drive the soft, zero Matsubara mode through a phase transition.
In this article, we are particularly interested in the case where this thermal evolution leads to a first-order phase transition. For there to be a first-order phase transition, there must be a coexistence of phases in some temperature range. Homogeneous phases of the real scalar theory may be distinguished by the value of the field condensate, where V 3 denotes the volume of space, and Z 3 denotes the partition function of the 3d EFT. The field condensate acts just as the density does in a liquid ↔ gas transition. At the critical temperature two different phases are equally likely to occur, their free energies being equal. In the case where the 3d effective field theory consists only of the real scalar field, there is a freedom to shift the field by a (temperature dependent) constant. Because of this freedom, one can find a basis in which g 3 (T ) = 0 for all T , which can be achieved by shifting after which the bare potential takes the form, where we have used that λ 3 δσ 3 = g 3 δm 2 3 , and have dropped the constant additive term. This shift reduces the general theory to a Z 2 symmetric Ising-like theory in the presence of a finite external fieldσ 3 (T ). Thus, if there is a phase transition, the critical temperature, T c , occurs at at which point the whole partition function possesses a Z 2 symmetry. Thus, at this point, either there are two identical phases related by the Z 2 symmetry, in which case there is a first-order phase transition, or there is only a single phase, in which case there is either a higher-order phase transition, or a crossover. Due to the symmetry, this is an exact equation determining the critical temperature: it is not corrected at any order in the loop expansion of the 3d EFT. By using the beta functions of the 3d effective parameters, Eqs. (2.29) to (2.32), one can see that this equality holds independently of the renormalisation scale. The nature of the transition, at this critical point, can be found from the value of the coefficient of φ 2 in Eq. (3.3) (minus the counterterm) which we have denoted by r. At tree-level the transition goes from first-order, through second order, and then crossover as the value of r goes respectively from negative, through zero, to positive. This is illustrated in Fig. 2. Beyond tree-level, the value of r for which the transition is second order, r * , shifts away from zero. In Ref. [106], a lattice Monte-Carlo study of the Z 2 -symmetric theory, it was found that, Combining this result with our previous arguments for the generic non-Z 2 -symmetric theory, we Figure 2: From left to right, the evolution of the tree-level potential in the 3d effective theory for first-order (r < r * ), second-order (r = r * ) and crossover (r > r * ) transitions. Here we adopt the basis in which the cubic term is zero, i.e. φ 3 → −g 3 /λ 3 + φ 3 , and choose µ 3 = 0.07860(18)λ 3 , such that r * (µ 3 ) = 0. In all cases the system starts at φ 3 < 0 at high temperatures and transitions to φ 3 > 0 as the system cools.
find that where r and r * are to be evaluated at the critical temperature, i.e. where Eq. (3.4) holds. The phase diagram of the theory is shown in Fig. 3. By using the beta functions of the 3d EFT, Eqs. (2.29) to (2.32), in the definition of r, one can see that these conditions giving the order of the phase transition hold independently of the renormalisation scale. Eq. (3.6) is exact, up to the statistical uncertainty in the determination of r * . As one considers weaker and weaker first-order phase transitions, i.e. as r tends towards r * from below, both the jump in the field expectation value and the screening mass tend towards zero. For asymptotically weak first-order transitions, the approach to the second order point is determined by universality, with the universality class being that of the 3d Ising model. 5 For example, the difference in the field condensate between the two phases, ∆ φ 3 , and the screening mass, m s , are given by with critical exponents β = 0.3258 (14) and ν = 0.6304(13) [108]. Note that, due to the Z 2 symmetry at the critical temperature, condensates of even powers of the field, are equal in both phases and hence, at the critical temperature, where, for clarity, we have shown the shift of Eq. (3.2) explicitly. For first-order phase transitions, condensates of odd powers of (φ + g 3 /λ 3 ) all have nonzero discontinuities across the transition.

Phase transition in perturbation theory
Perturbation theory is applicable rather generically to the nonzero Matsubara modes, as long as the theory is perturbative at T = 0. However the infrared physics of the zero Matsubara mode can become nonperturbative at high temperatures. In this section, we will investigate the applicability of perturbation theory to the infrared EFT, and apply it to the computation of various equilibrium thermodynamic properties. In particular, in Sec. 4.1 we will consider the general form of the loop expansion within the 3d EFT. Then in Sec. 4.2 we will compute the effective potential to threeloop order. Using this result, we will compute the free-energy and the discontinuity of the order parameter in Sec. 4.2, followed by the latent heat in Sec. 4.5.

The loop expansion near T c
For a consistent perturbative expansion, one should expand around a background field which is a minimum of the tree-level action. For homogeneous phases, this can be achieved by shifting the origin of the field to set σ 3 to zero. The couplings in the 3d effective theory both have positive mass dimension. As such, each successive loop order comes with either a factor of λ 3 or g 2 3 , as well as compensating powers of the tree-level mass parameter, making the dimensionless combinations We are interested in the case m 3 ∼ m 3,c . Using Eq. (3.4) we find that both loop-expansion parameters are of the same order, and the combination of parameters which dictates the convergence of the 3d loop expansion near the critical temperature is Shifting the field back to a generic point with nonzero σ 3 , the loop-expansion parameter becomes complicated by the introduction of cubic roots, but as this is simply a change of basis, it does not change the underlying physics. For a generic basis, but at the critical temperature, i.e. atσ 3 = 0, we are also able to find a simple form for the loop-expansion parameter, where v 0 is the field-value of the tree-level minima about the Z 2 -symmetric origin. From the above, it would appear that the loop expansion breaks down in the Z 2 -symmetric limit, g 3 , σ 3 → 0. This is related to the fact that the phase transition is of second order in the Z 2 -symmetric limit, and hence the scalar mass and the jump in the order parameter both go to zero. As these quantities enter the loop-expansion parameter inversely, the loop expansion breaks down in this limit. Thus, one must either find some other expansion parameter or resort to nonperturbative methods.
At the critical temperature in this model, the two phases are identical and hence the expansion parameter is the same in both phases. This is unlike the case of the Standard Model Higgs field, for which the presence of perturbatively massless gauge bosons causes the loop expansion to fail around the symmetric minimum, though it may converge well around the broken minimum [52,63,109].

The effective potential
In calculating the effective potential of the 3d EFT, we start with the bare Lagrangian, given in Eqs. (2.3) and (2.4). Unlike in Sec. 2, where we were only interested in UV physics at the scale ∼ 2πT , in studying the IR physics of the phase transition we must make a different split between free and interacting terms. In particular, the tadpole and mass terms now enter the free Lagrangian.
In the study of first-order phase transitions, the convex effective potential defined by a Legendre transform [110] is not particularly relevant. Instead it is appropriate to define the effective potential as the result of integrating over all non-zero momentum modes, following Ref. [111] (see also Ref. [72]). Only the 1-particle-irreducible diagrams contribute [111], of which there is one at one-loop, two at two-loop and six at three-loop order. The one-and two-loop diagrams are shown in Fig. 4, and the three-loop diagrams can be found in Ref. [112]. The necessary one-and two-loop integrals are collected in Appendix A.
For the effective potential up to two-loop order, we find - 15 -where, for the second and third derivatives of the tree-level potential, we have defined The required loop integrals are collected in Appendix A. At two-loop order, there is one logarithmically divergent diagram: the sunset. The corresponding counterterms are, Due to the superrenormalisability of the theory, the δσ 3 and δm 2 3 counterterms are in fact exact, to all orders in . However, the constant δV 3 counterterm will receive further contributions at four-loop order.
Note that these counterterms match the temperature-dependent divergences that we found in dimensional reduction, Eqs. (2.27) and (2.28). Further, by demanding that the bare parameters are independent of the scale Λ, we recover the beta functions of the 3d parameters that we found in dimensional reduction, Eqs. (2.29) to (2.32).
After inserting the results for the loop integrals from Appendix A, Eq. (4.5) becomes where µ 3 is the MS renormalisation scale of the 3d EFT. The Z 2 -symmetric limit of Eq. (4.8) agrees with that in Refs. [52,53]. The three-loop contribution can be constructed using the results of Ref. [112], and reads, where ξ = 0.03074157526289594 . . . is the result of performing a 1d integral numerically.

Condensates in the -expansion
As introduced in Sec. 3, the different phases of the system may be distinguished by spatially-averaged expectation values, or condensates, of the field operator or powers of the field operator. Condensates of operators which appear in the action can be conveniently computed by taking derivatives of the vacuum energy density, ε vac = − log Z 3 /V 3 , where Z 3 denotes the partition function of the 3d EFT and V 3 denotes the volume of space [63]. For example, for the field condensate we have Expectation values of higher powers of the field can be found by taking derivatives with respect to m 2 3 , g 3 and λ 3 . Taking derivatives with respect to the bare couplings yields UV divergent results. On the other hand, taking derivatives with respect to the renormalised couplings introduces counterterm corrections which result in finite results for the condensates.
In the following we consider expectation values of the form, where the subscripts b and s denote the broken and symmetric phases respectively. This language is used in analogy with symmetry-breaking transitions, although there is no symmetry which distinguishes between the phases here (except in the Z 2 -symmetric limit). The vacuum energy is the sum of all connected diagrams, whereas the effective potential is the sum of the connected, 1PI diagrams. In the vicinity of a minimum of the tree-level potential, the 1particle-reducible (1PR) diagrams, which are missing from the effective potential, can be generated from the 1PI diagrams by performing a strict expansion [113] (see also Refs. [63,64,109,114]). In this, one expands the effective potential and the vev in powers of , and solves for the vev order-by-order in . In gauge theories, this approach gives rise to gaugeinvariant results order-by-order in [113,115]. However, for radiatively-induced transitions, whereby tree-level, O( 0 ), and loop-level, O( ), contributions compete on a par, the expansion breaks down. For instance, in the SU(2)-Higgs model at two-loop order in the -expansion, the critical temperature is IR divergent and the latent heat receives an unphysical imaginary part [109].
In calculating the condensates, it is convenient to shift the origin of the field according to Eq. (3.2), so that the third derivative of the tree-level potential vanishes, About this origin, and at the critical temperature the tree-level minima are located at v (0) = ±v 0 , where v 0 is given by Eq.
where all the functions on the right-hand sides of the equations are evaluated at the tree-level minimum. Evaluating these, we find The explicit log µ 3 term at two-loop order cancels the running of v 0 at order 2 , and further the absence of a log µ 3 term at three-loop order is a result of the absence of running of the one-loop correction to ∆ φ 3 , which only depends on λ 3 . Numerically Eq. (4.16) reads where we have introducedμ 2 3 = µ 2 3 /(3λ 3 v 2 0 ), and have indicated the size of the four-loop corrections. The loop-expansion parameter, α 3 , is given in Eq. (4.3). As can be seen, the expansion coefficients are all O(1), suggesting that the magnitude of α 3 should give a reliable estimate of how well the series converges. As α 3 is inversely proportional to the jump in the order parameter, the series converges more quickly for stronger transitions.

Renormalisation group improvement
The effective potential, or rather its φ 3 -derivative to avoid the cosmological constant, satisfies the renormalisation group (or Callan-Symanzik) equation, where the nonzero beta functions are taken from Eqs. (2.31) and (2.32). Due to the superrenormalisability of the theory, this equation is exact. Expanding it in , leads to an infinite set of relations linking the running of couplings at O( n ) to explicit logarithms at O( n+2 ). Using this, one can deduce the explicit µ 3 dependence present at four-and five-loop order from the running of the couplings at two-and three-loop order, However, Eq. (4.18) is an exact equation, and hence it seems reasonable to solve it exactly, rather than order-by-order in . Doing so resums the most ultraviolet sensitive higher order contributions, giving the renormalisation group improved (RGI) effective potential. By incorporating this nonperturbative information, one would hope to improve the accuracy or convergence of perturbative results. The construction of the RGI effective potential was presented in this context in Ref. [52], though we will adopt an alternative approach, as follows.
Correlation functions, such as the linear condensate ∆ φ 3 , satisfy an identical renormalisation group equation. This can be solved order-by-order in , to find the scale dependent parts of the 4and 5-loop contributions, where µ 3,0 is some initial scale, at which v 0 and α 3 are given. This defines our RGI perturbative calculation. Note that this improvement is naturally incorporated upon performing dimensional reduction, if the exact running of the parameters is used to cancel the dependence on the matching scale Λ, replacing it with a new renormalisation scale µ 3 ; see the end of Sec. 2.1. Interestingly, while both g 3 and λ 3 are independent of µ 3 , the expansion parameter α 3 is not. While this is clear in the form given in Eq. (4.3), to see it in the form given in Eq. (4.2) one should note that changing µ 3 shifts the minimum of the potential, and the basis transformation required to shift back to remove the tadpole induces a µ 3 -dependence in the cubic coupling.

The latent heat
The latent heat of the transition, L, can be determined from the following thermodynamic relation, where f is the free energy density. Noting that ∆f = T ∆ε vac , and using the chain rule, we can expand this relation in terms of the field condensates, 24) all evaluated at the critical temperature, at which point ∆f = 0. This expression for the latent heat can be simplified significantly by shifting the origin of the field following Eq. (4.4). In this basis, the coefficient of the cubic term vanishes by construction and the even power condensates vanish at the critical point due to the Z 2 symmetry; see Eq. (3.3). Thus, the latent heat simplifies to,

Phase transition on the lattice
Perturbation theory can only take us so far. As discussed in Sec. 4.1, at the critical temperature the loop-expansion parameter for the 3d EFT scales as α 3 ∝ λ 3/2 3 /|g 3 |. So, as one approaches the Z 2 -symmetric limit, the loop-expansion parameter grows without bound, signalling the complete breakdown of perturbation theory.
To reliably study the phase transition for both small and large α 3 , we resort to lattice Monte-Carlo simulations; for relevant overviews see Refs. [70,[116][117][118], and for recent studies in other models see Refs. [89,119,120]. The first step, which we carry out in Sec. 5.1, is to find explicit relations between the bare parameters of the lattice Lagrangian, and those of the continuum theory in the MS renormalisation scheme. Once this is done, measurements from Monte-Carlo simulations can be directly interpreted in terms of MS observables.
In Sec. 5.1 we derive the lattice-continuum relations, applicable in the limit a → 0. Following this we discuss the Monte-Carlo simulations in Sec. 5.2, with details of the algorithms given in Appendix D. Results for the latent heat, extrapolated to the continuum limit, are presented in Sec. 5.3, with additional plots of the continuum extrapolations given in Appendix E.

Lattice-continuum relations
The simplest lattice Lagrangian which converges to Eq. (2.3) in the continuum is where x labels lattice sites, i denotes the link from the site x to its nearest neighbour in one Cartesian direction, and the sum over i denotes a sum over Cartesian directions. 6 We have added the subscript L to the tadpole and mass counterterms because, due to the different regularisation of divergent integrals, these quantities differ from their MS counterparts. Conversely, due to the absence of momentum-dependent divergences, or divergences of diagrams with three or four external legs, the field renormalisation and three-and four-point couplings may be chosen to be equal to their MS counterparts. Due to the superrenormalisability of the theory, all divergences of the theory turn up at finite loop order and hence can be calculated analytically, in lattice perturbation theory. Thus, it is possible to derive the exact relationship between the parameters in the lattice and continuum theories, in the limit a → 0. To do so, we follow Refs. [64,65] in computing the effective potential in lattice perturbation theory, and equating the result to the effective potential in the continuum, in the MS scheme.
The computation of the effective potential in lattice perturbation theory, mirrors almost exactly that in Sec. (4.2). The diagrams and combinatorics are the same; only the values of the loop integrals are different. A notable consequence of this difference is that with lattice regularisation, unlike with dimensional regularisation, there are linear divergences at one-loop, which demand compensating one-loop counterterms.
The one-loop order contribution to the effective potential is where the lattice loop-integral J is defined in Eq. (171) of Ref. [63]. In the infinite volume limit, this lattice loop integral may be carried out as a Fourier integral over the first Brillouin zone, the result of which can also be found in Ref. [63].
Expanding for small a, we find the following one-loop counterterms, were Σ is a numerical constant, the result of a dimensionless integral. It is given analytically in Eq. (170) of Ref. [63], and its numerical value is Σ = 3.17591153562522.... Progressing to two-loop order involves nothing qualitatively new. The lattice-regularised, twoloop effective potential is equal to its MS counterpart in the a → 0 limit if we choose the following counterterms (including their finite parts), Here ζ is another numerical constant, the value of which is approximately ζ ≈ 0.0884801 (called C 3,U in Ref. [68], see also Ref. [65]). There are also divergent additive contributions to the (cosmological) constant counterterm at three-and four-loop order, but there are no divergent, field-dependent contributions at higher-loop orders. Thus, the two-loop results for the tadpole and mass lattice counterterms are exact in the a → 0 limit. Finally, we would like to comment on the lattice-continuum relations for the four special condensates, defined as in Eq. (4.10). Now that the lattice-continuum relations have been found for the bare lagrangians, those for the condensates follow simply by differentiation with respect to the renormalised MS parameters [63]. Due to the absence of σ 3 -dependence in the lattice-continuum relations above, the result for the one-point condensate is simple, as a → 0.
-20 -  Figure 5: Probability distributions, p, of lattice Monte-Carlo measurements ofφ 3 at the critical temperature. Fig. (a) shows the dependence on the lattice volume at fixed lattice spacing, a. Fig. (b) shows the dependence on the lattice spacing, a, at fixed lattice volume. The coexistence of clearly separated phases shows the first-order nature of the transition.

Monte-Carlo simulations
With the lattice-continuum relations in hand, one can carry out Markov-chain Monte-Carlo simulations. For this we wrote a simulation code in C. Starting from an arbitrary initial lattice field configuration, we perform successive updates on the field such as to sample the Boltzmann-weighted probability distribution of the theory. The list of configurations forms the Markov chain. After a large number of updates of the whole lattice (typically a few million), the Markov chain converges to a relatively accurate representation of the Boltzmann probability distribution. From this measurements of physical quantities, such as the field condensates, can be made. In order to ensure the field configuration does not get stuck in one phase, we utilised multicanonical methods [116,121]. Details of the algorithms we adopted are discussed in Appendix D. Binning measurements ofφ 3 to form a histogram gives an approximation to the probability distribution for this observable; see Fig. 5. At the critical temperature this probability distribution shows a two-peaked structure, characteristic of a first-order phase transition. To extract physical results relevant for continuum physics requires taking first the infinite volume limit, and then the zero lattice spacing limit. In Fig. 5a we show the effect on the probability distribution of increasing the lattice volume, while keeping the lattice spacing fixed, and in Fig. 5b we show the effect of decreasing the lattice spacing while keeping the lattice volume fixed.
At the critical temperature there is only one dimensionless parameter which characterises the 3d EFT, and this can be chosen to be g 3 /λ 3/2 3 (alternative choices are (r − r * )/λ 2 3 and α 3 ). This is because, of the four initial Lagrangian parameters, one parameter can be fixed by shifting the field origin, σ 3 = 0 say, a second parameter can be fixed by the condition of being at the critical temperature, m 2 3 say, and a third parameter can be fixed by a choice of units, λ 3 = 1 say. This is analogous to what happens in the SU(2)-Higgs theory at the critical temperature, where the dimensionless parameter which controls the character of the phase transition is x ≡ λ Higgs,3 /g 2 SU(2),3 [70]. In our case the most closely corresponding dimensionless combination would be x singlet ≡ λ 3 /|g 3 | 2/3 .
We chose six parameter points to simulate on the lattice, starting at g 3 /λ 3/2 3 = −1.2 and decreasing by powers of 2 until g 3 /λ 3/2 3 = −0.0375. In all cases, the MS renormalisation scale was -21 -taken to be µ 3,L = 1.066496λ 3 . 7 To extrapolate to the continuum limit for each such parameter point, we simulated between five and eight different lattice spacings, and for each lattice spacing we simulated between five and eight different lattice volumes. Thus, all in all, our data set consists of more than 200 different simulations. The continuum-extrapolated results are independent of the renormalisation scale.
The different lattice spacings, a and volumes, L 3 , should all be chosen to satisfy a ≪ ξ ≪ L , (5.9) where ξ is the correlation length of the system, or the inverse screening mass. For parameter values where the EFT is perturbative, the screening mass will be close to the tree-level mass, in which case ξ ≈ 1/m 3 . Including also the one-loop corrections, for σ 3 = 0, we find In principle the two-loop corrections can be constructed from the results of Ref. [112], though we have not done so. In the opposite limit, near the second-order phase transition where perturbation theory does not work, Eq. (3.8) should provide a better estimate. In choosing appropriate lattice spacings and volumes, we have used a combination of these two estimates. However, these approximations are not perfect and, especially for small values of −g 3 /λ 3/2 3 , additional lattice spacings and volumes were necessary to attain reasonable continuum limits. A more robust alternative, which we did not attempt, would be to directly measure ξ on the lattice, by the exponential decay of correlation functions with distance; see for example Ref. [70].
Our methods for analysing the simulation data are fairly standard, and generally follow Refs. [70,118]. Error bars for simulation data points show statistical errors, calculated using jackknife resampling on blocked measurements, with each block much larger than one autocorrelation time.

Latent heat on the lattice
The latent heat is proportional to the change in the linear field condensate, ∆ φ 3 , with the proportionality constant being dependent on the details of the full 4d theory, but not on the dynamics of the 3d EFT; see Sec. 4.5. Thus, for a theory which is weakly coupled at T = 0, this proportionality constant can be calculated perturbatively. The possibly nonperturbative dynamics of the zero Matsubara mode only enters the latent heat through ∆ φ 3 . In a finite volume, the linear condensate can be defined as where P (φ 3 ) denotes the probability density of being in a stateφ 3 andφ min 3 = −g 3 /λ 3 denotes the position of the minimum of the probability density between the two phases. In a small enough volume, states between the two phases are not uncommon, however as the lattice volume grows such states become exponentially rarer; see Fig. 5a.
In Fig. 6 we show our results for g 3 /λ 3/2 3 = −1.2 and −0.6. Plots for the other four parameter points are collected in Appendix E. The continuum-extrapolated results are collected in Table 1.
Figs. 6a and 6b show the infinite volume extrapolations at fixed lattice spacing. As the theory is gapped, at large volumes the infinite volume limit is approached exponentially fast, with corrections ∼ exp(−L/ξ). The lines shown in these figures show least-squares fits to the data of the form (d) Figure 6: Taking the infinite volume, followed by the zero lattice spacing limits of lattice data for ∆ φ 3 at two parameters points. For the fits to the lattice spacing dependence here, we show both a linear fit excluding the data point with largest a, and a quadratic fit to all the data. In Fig. 6c the linear and quadratic fits give reduced χ 2 ≈ 1.5 and 0.9 respectively. In Fig. 6d these are instead 5.6 and 1.5. We take the difference between these fits as a measure of the systematic error introduced in the extrapolation.  Table 1: Continuum-extrapolated lattice results for the discontinuities in the linear and cubic condensates. The UV divergences of the cubic condensate has been removed, following Ref. [63]. Errors quoted are statistical followed by systematic.
where the c i are fit parameters. In all cases the reduced χ 2 values for the fits are in the range ∼ [0. 2,5]. As a measure of the systematic uncertainty associated with the infinite volume extrapolations, we repeat the succeeding analysis using only the largest volume simulations at each value of a, rather than the extrapolation of the fit. The results of the infinite volume extrapolations are the data points in Figs. 6c and 6d, with the error bars resulting from the least-squares fit. The lines in these figures in turn show the extrapolations to zero lattice spacing. As a consequence of the lattice-continuum relations of Sec. 5.1, for sufficiently small am 3 , the lattice results should approach the continuum limit linearly, where we have denoted the MS parameters of the theory by κ = {σ 3 , m 2 3 , g 3 , λ 3 , µ 3 }. For each parameter point we attempt to extrapolate to a = 0 linearly in a, but consider also higher powers in a Taylor expansion in a if the linear fit is poor. In such cases we truncate at the lowest power which gives a reduced χ 2 of order one.
Some, but not all, of the lattice spacing fits reveal the presence of significant nonlinear terms in am 3 . By comparing the six different parameter points, we suggest an explanation in terms of the coefficients C i (κ) in Eq. (5.12). For the strongest transition, g 3 /λ 3/2 3 = −1.2, the fit in Fig. 6c shows C 1 > 0. Progressing to weaker transitions, shown in Figs. 6d, 13c, 13d, 14c and 14d, one can see that C 1 changes sign around g 3 /λ 3/2 3 ≈ −0.6, remaining negative all the way to g 3 /λ Fig. 6c suggests the curvature is negative, C 2 < 0, for g 3 /λ 3/2 3 = −1.2. Likewise progressing to weaker and weaker transitions, one can see that C 2 changes sign somewhere between g 3 /λ 3/2 3 = −0.6 and g 3 /λ 3/2 3 = −0.075. These considerations can explain, for example, the comparatively poor linear fit in Fig. 6d as due to the smallness of C 1 at this parameter point, and the comparatively good linear fits in Figs. 13c and 13d as due to the smallness of C 2 . For the parameter points where a nonlinear fit is shown, we use the difference between this and the linear fit as a measure of the systematic uncertainty in the extrapolation. The relatively large systematic errors introduced by the extrapolation suggest that for future studies, it would be worthwhile to use the O(a) improved lattice-continuum relations [66][67][68][69].

Discussion
In this work, we have studied the phase transitions of a generic real scalar field. We have been agnostic about couplings to other fields, such as to the Higgs, to sterile neutrinos and to gravity. However, we have focused on the case where the infrared dynamics of the phase transition is dominated by the real scalar; the contributions of all other particles to the transition being simply to modify the effective couplings of the infrared EFT. With this restriction, we have characterised the order and strength of the phase transition nonperturbatively.
Our results can be applied to a variety of 4d particle physics and cosmological models. As long as the phase transition is dominated by a real scalar field, no new simulations need be performed to nonperturbatively determine the phase diagram and latent heat of the 4d model in question. This is a key advantage of the EFT approach. One needs only to compute the matching relations to the 3d EFT, discussed with examples in Secs. 2.2 and 2.3. In a similar way the lattice results of the SU(2)-Higgs 3d EFT [70,74,122], have been repurposed to determine nonperturbatively parts of the phase diagram (and other quantities) for the xSM [40], the two higgs doublet model (2HDM) [55] and the triplet scalar extended SM (ΣSM) [57].
We leave for future work the computation of the bubble nucleation rate within the 3d EFT, from which one can determine the nucleation temperature, the duration of the phase transition and the change in the trace-anomaly through the transition. These in turn are crucial ingredients in determining the gravitational wave spectrum resulting from a first-order phase transition.
-24 -   ). In Fig. 7a, no uncertainty bands are plotted for the tree and 1-loop approximations as there is no explicit µ 3dependence at these orders. While both unimproved and RGI perturbation theory agree well with the lattice results at small couplings, the unimproved approach breaks down badly above α 3 ≈ 1, whereas the RGI approach continues to work remarkably well.
-25 - Theoretically, an important result of this work is a quantitative answer to the question: how well does perturbation theory fare in this model, in comparison to nonperturbative lattice results? We can answer this question as a function of the loop-expansion parameter, which is also the only dimensionless parameter on which the strength of the phase transition depends nontrivially. Our answer is shown in Figs. 7 and 8, which also summarise our results. Fig. 7 shows the lattice and perturbative results across the full range of parameter space that we have studied. The horizontal axis is the loop expansion parameter of the 3d EFT, α 3 , as defined in Eq. (4.2), and approximately equal to α 3 ≈ λ 3/2 T /(4π|g|). 8 In Fig. 7b the perturbative results plotted use the renormalisation group improvement of Sec. 4.4, whereas those in Fig. 7a use unimproved perturbation theory. Results at three different renormalisation scales are plotted in order to estimate the uncertainty due to missing higher-loop orders.
For small α 3 , there is very good agreement between lattice and perturbation theory, at least when the one-loop term is included. For the smallest two values of α 3 , the 2-and 3-loop perturbative results lie entirely within the lattice error bands. By α 3 ≈ 0.25, though the agreement is still good, the scale dependence has grown significantly, and there are small discrepancies with the lattice result, even at 3-loop order. The unimproved and RGI perturbative calculations differ little at these smaller values of α 3 .
At larger α 3 , the unimproved perturbative results deviate increasingly from the lattice, and at around α 3 ≈ 1 this perturbative calculation breaks down altogether, showing a strong dependence on the renormalisation scale. At these large values of α 3 , higher-loop orders do not improve the situation. By contrast, the RGI perturbative results remain remarkably under control all the way up to the largest expansion parameter we study, α 3 ≈ 2. The scale dependence of the RGI results is much smaller, 9 the agreement with the lattice results is much better, and each additional loop order improves the agreement. Fig. 8 shows this more clearly: it shows that the 3-loop RGI calculation has an accuracy of a few percent even at α 3 ≈ 2.
The apparently miraculous convergence of RGI perturbation theory at expansion parameters as large as α 3 ≈ 2 deserves explanation. One perhaps natural explanation is that this is due to the superrenormalisability of the 3d EFT, and the consequent exactness of the renormalisation group equations. Nevertheless, this point deserves further study in other superrenormalisable 3d theories.
Overall, the discrepancies found here between perturbation theory and the lattice are smaller than those that have been found in studies of non-Abelian gauge theories [70,74,89,119,122]. Two possible explanations for this are: (i) that there are no perturbatively massless particles in first-order phase transitions in this model, and hence no true infrared divergences, unlike in non-Abelian gauge theories [47], (ii) that there is a tree-level barrier between the phases in this model, unlike in gauge-Higgs theories, and hence an -expansion is amenable. Although both explanations surely play a role, an understanding of their relative importance could inform theoretical studies attempting to reduce uncertainties in calculations of first-order phase transitions. Such an understanding could be achieved, for example, by performing a similar lattice study around the electroweak scale in the xSM, including couplings to the electroweak sector in the 3d EFT.

A Loop integrals
Our notation for momenta and loop integration follows Ref. [53,62]. In particular, thermal fourmomenta are denoted by uppercase letters, P = (p 0 , p), their components being the Matsubara frequencies, p 0 = 2πT n, and the spatial momenta, p. Their norms squared are P 2 = p 2 0 + p 2 , where p 2 ≡ p.p. The loop integration measure in d = 3 − 2ǫ dimensions is defined as, which includes powers of the MS renormalisation scale, Λ, to make the measure up to mass dimension 3. Sum-integrals over loop momenta are then defined as, The necessary one-loop integrals in the effective theory are, In fact this is true even if we vary the renormalisation scale over a factor of 10 larger range (i.e. a factor of 100).

-27 -
The only two-loop integral in the effective theory that we need is the sunset integral [52], In the full theory, at one-loop we have the following master sum-integral, We also note the result for the logarithmic integrand, Massless two-loop sum-integrals in the full theory can be most simply calculated using integrationby-parts techniques, which reduces them to products of the one-loop sum-integrals above [123,124]. The relevant sum-integrals for us are, 10 These sum-integrals have also been calculated directly in Refs. [53,125,126].

B Correlation functions
In this appendix, we give the results for the connected, 1PI correlation functions in the minimal real singlet scalar model. The one-and two-point functions are calculated to two-loop order, and the three-and four-point functions are calculated to one-loop order, which are input to the dimensional reduction matching relations discussed in Sec. 2. Results for all the sum-integrals that appear are given in Appendix A. As we are interested only in zero Matsubara mode external legs, we will explicitly show only the dependence on the spatial momenta. The two-point function is given, up to two-loop order, by the Feynman diagrams in Fig. 9. In 10 We thank Juuso Österman and Philipp Schicho for bringing this to our attention, and for verifying the results. Figure 9: Feynman diagrams contributing to the two-point function up to two-loop order. They are shown here in the same order that they appear in Eq. (B.1).
the full theory, these are where we have expanded assuming p ∼ √ λT . As regards the momentum dependence, this justifies the following expansion, In accordance with the philosophy of the strict perturbative expansion, this expansion projects out any non-analytic IR contributions. In our scheme in which the mass term is treated as a perturbation, the self-energy, Π, is equal to the two-point function up to the p 2 term, Figure 10: Feynman diagrams contributing to the three-point function up to one-loop order. They are shown here in the same order that they appear in Eq. (B.6).
(a) (b) (c) (d) Figure 11: Feynman diagrams contributing to the four-point function up to one-loop order. They are shown here in the same order that they appear in Eq. (B.9).
In the effective theory the corresponding calculation is trivial, due to the vanishing of scaleless integrals in dimensional regularisation, The vanishing of momentum-dependent loop contributions follows an expansion identical to Eq. (B.3) for Γ 3 . Just as for the one-point function, the scaleless integrations in the 3d effective theory exactly match scaleless integrals for the zero Matsubara modes in Eq. (B.1).
The three-point function is given, up to one-loop order, by the Feynman diagrams in Fig. 10. In the full theory, these are Corrections related to the soft external momenta arise at higher order than we consider. In the effective theory the corresponding calculation is again trivial, Γ The four-point function is given, up to one-loop order, by the Feynman diagrams in Fig. 11. In the full theory, these are Corrections related to the soft external momenta arise at higher order than we consider. In the effective theory the corresponding calculation is again trivial, Γ

C Some relations at zero temperature
To preserve the accuracy achieved in our dimensional reduction throughout the calculation, it is necessary to: (i) relate the MS parameters to physical quantities at some input scale at one-loop order [54,61] and (ii) run the MS couplings from the input scale up to some temperature-dependent scale chosen to minimise large logarithms. The details follow. The matching of physical quantities to MS parameters is carried out at zero temperature at some input renormalisation scale which we take to be equal to the physical (pole) mass of the particle µ input = m phys . The tadpole and three-and four-point MS couplings can be fixed by requiring at one-loop order. At zero temperature, the one-loop contribution to the MS -renormalised effective potential is given by where M 2 (φ) = m 2 + gφ + 1 2 λφ 2 . In principle one should in fact match the three-and four-point couplings on-shell, rather than at zero external momentum, but in lieu of a measurement of these couplings Eqs. (C.2) and (C.3) are reasonable renormalisation conditions.
For the MS mass parameter, one should match the pole mass to the physical mass, where, approximating m 2 ≈ m 2 phys within the one-loop self-energy, we have Π(m 2 ) = 2(4π) 2 λm 2 1 + log Here we have included only the 1PI part of the self-energy as the 1PR part cancels by virtue of Eq. (C.1). The necessary (zero temperature) counterterms in dimensional regularisation are By demanding that the bare coefficients are independent of the cut-off scale, we derive the beta functions, β X ≡ dX /d log Λ, up to the same order: If we assume the parametric scaling relation λ ∼ g/m, then the two-loop O( 2 ) parts of the tadpole and mass beta functions do not contribute at the order we work.

D Monte-Carlo update algorithm
A vanilla Monte-Carlo update strategy is not well suited to the study of first-order phase transitions as it is unable to efficiently sample both phases. For strong transitions, the significant barrier between the two phases will mean that any practical run will get stuck in one of the two phases. A standard solution to this problem is the multicanonical method [116,121], which takes advantage of the freedom to change the probability measure thuswise, Here W is a weight function, chosen in order to increase the probability of field configurations between the two phases. The weight function is taken to depend on a single variable, O[φ 3 ], and itself is calculated in a Monte-Carlo simulation. For the choice of O[φ 3 ], the key requisites are that it should distinguish between the two phases, be relatively quick to calculate, and its use should enhance as much as possible the rate at which the simulation tunnels between the two phases. We investigated a variety of field polynomials, and settled with the spatial average of (φ 3 + g 3 /λ 3 ) 3 . The specific iterative multicanonical algorithm that we adopted follows Ref. [87]. Algorithm (i) is standard and (ii) follows Ref. [70] closely. We used the checkerboard update order for (i) and (iii), with a random choice of whether the first sweep was over the odd or the even sites. The overrelaxation algorithm, (iii), for the update φ i → φ ′ i starts by solving the equation, 11 Here we drop the subscript 3 on the field, to avoid notational clutter.
-32 -  The continuum-extrapolated changes in the quadratic and quartic condensates across the transition, with statistical and systematic uncertainties in the final digits following in parentheses. The UV divergences of the quartic condensate has been removed, following Ref. [63]. The results should be exactly zero, by symmetry, thus providing a consistency check on the simulation and analysis pipeline, including the error estimates.
where S i is the part of the action which depends on the field at lattice site i. This quartic equation for φ ′ i can be reduced to a cubic equation by factorising the solution φ ′ i = φ i . The first step in the overrelaxation algorithm is to choose for φ ′ i one of the real solutions of the cubic equation. In the case where there are three real solutions to Eq. (D.2), we choose each solution with equal probability (though in practice this case is vanishingly rare). For this step we thus have that, where we have used that for this step the probability densities are equal, p(φ i → φ ′ i ) = p(φ ′ i → φ i ). To ensure detailed balance, we must then take a second step, which we perform as a Metropolis, accept/reject on the measure, The two-step algorithm then has equal probability of going forward or backwards and hence, due to Eq. (D.2), satisfies detailed balance. For one specific parameter point, we measured the autocorrelation times of Markov chains using combinations of these three algorithms. The autocorrelation functions are shown in Fig. 12. Note that neither the overrelaxation nor the global, radial update are ergodic, and hence they must be combined with the local Metropolis update. When combined with the overrelaxation and Metropolis, the global, radial update does not further diminish the autocorrelation time, and hence in the final runs we did not use it. Taking into account also the computational cost of the algorithms, our final update algorithm consisted of 1 local Metropolis update followed by 4 overrelaxation updates.
To validate the simulation code, we performed a range of tests. The implementations of the action and its various terms were shown to converge towards exact analytic results for specific smooth field configurations. Our final results at small loop expansion parameters α 3 ≪ 1 agree well with perturbation theory; see Fig. 7. Our results for the condensates of even powers of the field are close to zero, as required by Eq. (3.9), with all discrepancies being consistent with the estimated statistical and systematic errors. This is shown in Table 2.
At a small subset of parameter points we performed additional checks. The results of different combinations of algorithms, both multicanonical and canonical, were shown to agree within errors. Our chosen random number generator, the Tausworthe generator of L'Ecuyer [127,128] (implemented in GSL [129] as gsl_rng_taus2), was shown to produce results in agreement with two others: the Mersenne Twister algorithm [130] and the RANLUX algorithm of Luscher [131,132].  Figure 13: Taking the infinite volume, followed by the zero lattice spacing limits of lattice data for ∆ φ 3 at two parameters points. For the fits to the lattice spacing dependence in Figs. 13c and 13d, we perform linear fits to all the data. This gives reduced χ 2 ≈ 3 in both cases.
Finally, we would like to thank Kari Rummukainen, who shared his Monte-Carlo simulation code for the cubic anisotropy model [68], and with which we found agreement using a version of our code which implements that model.

E Additional continuum extrapolations
In this appendix we collect plots of the remaining lattice-continuum extrapolations which were not given in the body of the text, Figs. 13 and 14. In all cases the exponential ansatz for the volume dependence (see Sec. 5.3) appears to fit the data well. Further, for each value of a this extrapolation differs little from that of the largest volume. For the stronger phase transitions of Fig. 13, a linear fit to the a-dependence of all data points gives a reduced χ 2 ≈ 3 in both cases, appearing to fit the data reasonably well. However, for the weaker phase transitions of Fig. 14, there are significant deviations from linearity at the larger lattice spacings. This is to be expected because, according to Eqs. (5.10) and (3.8), for small |g 3 |/λ 3 the screening mass grows significantly larger than the treelevel mass, meaning a/ξ ≫ am 3 . To remedy this difficulty, we have simulated additional smaller -34 -  Figure 14: Taking the infinite volume, followed by the zero lattice spacing limits of lattice data for ∆ φ 3 at two parameters points. For the fits to the lattice spacing dependence here, we show both a linear fit to the three data points with smallest a, and a cubic fit including three additional points. In Fig. 14c the linear and cubic fits give reduced χ 2 ≈ 0.5 and 0.2 respectively. In Fig. 14d these are instead 0.02 and 2.
lattice spacings. In each of Figs. 14c and 14d we show two fits. The difference between the results of these two fits gives a measure of the systematic uncertainty in the extrapolations.