Standard Model thermodynamics across the electroweak crossover

Even though the Standard Model with a Higgs mass mH = 125 GeV possesses no bulk phase transition, its thermodynamics still experiences a"soft point"at temperatures around T = 160 GeV, with a deviation from ideal gas thermodynamics. Such a deviation may have an effect on precision computations of weakly interacting dark matter relic abundances if their mass is in the few TeV range, or on leptogenesis scenarios operating in this temperature range. By making use of results from lattice simulations based on a dimensionally reduced effective field theory, we estimate the relevant thermodynamic functions across the crossover. The results are tabulated in a numerical form permitting for their insertion as a background equation of state into cosmological particle production/decoupling codes. We find that Higgs dynamics induces a non-trivial"structure"visible e.g. in the heat capacity, but that in general the largest radiative corrections originate from QCD effects, reducing the energy density by a couple of percent from the free value even at T>160 GeV.


Introduction
The excellent performance of the Standard Model (SM) in describing LHC data suggests that the SM represents a precise description of nature up to energy scales of several hundred GeV. If this observation continues to be confirmed by future runs, one consequence is that the Higgs phenomenon set in smoothly, i.e. without a phase transition, as the Universe cooled down to temperatures below 160 GeV [1]- [4] (cf. ref. [5] for a review of refined investigations). This would mean that the bulk motion of the plasma did not depart from thermal equilibrium and therefore did not generate cosmological relics.
Nevertheless, the thermodynamics of the SM could still have left a "background" imprint on other non-equilibrium physics that just happened to be going on. One example is that the B+L violating rate switched off rapidly around the crossover [6], and therefore determined which fraction of a given lepton number being produced around these temperatures could be converted into baryons (cf. e.g. ref. [7]). Another is that if Dark Matter particles decoupled at around this period, even small features in the equation of state could have had an impact [8].
(Similar considerations can be carried out for the QCD crossover, cf. refs. [9,10] and references therein.) A rough estimate for the decoupling temperature of weakly interacting particles of mass M is T ∼ M/25, so that we may expect the largest sensitivity for M ∼ 4 TeV or so.
The purpose of our study is to estimate the equation of state of the SM around the electroweak crossover. We do this through perturbative computations extending up to 3-loop level, as well as by making use of existing lattice simulations within a dimensionally reduced effective theory. Only the temperature is treated as a non-trivial canonical variable, with chemical potentials associated with conserved charges such as the baryon minus lepton number or the hypercharge magnetic flux set to zero.
The paper is organized as follows. After defining the basic observables in sec. 2, we derive a "master equation" in sec. 3 which relates the trace of the expectation value of the energymomentum tensor to three ingredients. The first ingredient, scale violations through quantum corrections, is addressed in sec. 4. The second ingredient, the expectation value of the Higgs condensate, needs to be evaluated up to the non-perturbative level, and this is achieved in sec. 5. The third ingredient is related to vacuum renormalization and is discussed in sec. 6. All ingredients are put together in sec. 7, where we also present phenomenological results. Section 8 contains some conclusions as well as a discussion of the theoretical uncertainties of the current analysis. Readers not interested in technical details may start from sec. 7.

Basic setup
We consider the Standard Model with a Higgs potential of the form where ν 2 B , λ B are bare parameters (the corresponding renormalized parameters are denoted by ν 2 , λ), and L E is a Euclidean (imaginary-time) Lagrangian. Denoting by g 2 ∈ {λ, h 2 t , g 2 1 , g 2 2 , g 2 3 } a generic coupling constant 1 , we concentrate on the parametric temperature range in the following. At the lower edge of this range, the effective Higgs mass parameter of the dimensionally reduced theory [11,12],m 2 3 [13], is assumed to satisfy but it can have either sign. If it is negative, we may expect to find ourselves in a "Higgs phase". Within a gauge-fixed perturbative treatment, this would imply that the Higgs field has an expectation value v 2 ∼ −m 2 3 /λ > 0. Within the range of eq. (2.3), only relatively small expectation values v < ∼ g 1 2 T can be considered, where we counted λ ∼ g 2 . It is well known that when momenta in the range |m 2 3 | ∼ (g 2 T /π) 2 are considered, which is a special case of eq. (2.3), then the dynamics of the system is non-perturbative [14,15]. Therefore the dynamics needs to be treated with lattice methods. However, non-perturbativity is associated with particular modes only, and can be captured with a dimensionally reduced effective description [11,12]. Even though the construction of this theory is perturbative, we expect to obtain a description accurate on the percent level within a weakly coupled theory such as the SM [16,17].
The basic observable that we consider is the thermodynamic pressure of the SM. The pressure can formally be defined through the grand canonical partition function Z as where the thermodynamic limit V → ∞ is implied, and p B denotes the bare result. The computation of p B has previously been considered in ref. [18] for the casem 2 3 ∼ g 2 T 2 , and in ref. [19] form 2 3 ∼ g 3 T 2 /π. Following refs. [18,19] 2 (which were inspired by ref. [20]) we write where p E collects contributions from the momentum scale k ∼ πT , p M those from k ∼ gT , p G those from k ∼ g 2 T /π. Rephrasing terminology often used in the QCD context, we refer to the effective theory contributing to p M as ESM ("Electrostatic Standard Model") and to that contributing to p G as MSM ("Magnetostatic Standard Model").
As defined by eq. (2.4) the pressure is ultraviolet divergent. We renormalize it by assuming that the pressure (and energy density) vanish at T = 0. The renormalized pressure can then be written as Our goal is to determine the dimensionless functions p(T )/T 4 and T ∂ T { p(T )/T 4 } up to O(g 5 ). For the latter quantity, it turns out that the contribution of the softest momentum modes needs to be treated non-perturbatively in order to reach this precision.

Master equation
Consider a normalized form of eq. (2.6), Hereμ denotes the MS renormalization scale. We envisage that at some temperature (T = T 0 ) on either side of the crossover, say T 0 ≪ 160 GeV or T 0 ≫ 160 GeV so that |m 2 3 | > ∼ g 3 T 2 /π, p/T 4 can be determined by a direct perturbative computation. The task then is to integrate p/T 4 across the electroweak crossover: So, we need to compute the logarithmic temperature derivative of the dimensionless ratio p/T 4 . The result is non-zero because of the breaking of scale invariance, either by the explicit mass term ν 2 , or by quantum corrections related to running couplings. Making use of standard thermodynamic relations, we note that where e denotes the energy density. In the context of QCD this quantity is referred to as the trace anomaly, and has been studied with lattice methods. Were it not that chiral fermions (in particular the top quark) are essential for the physics considered, the problem could in principle be studied with full 4-dimensional lattice simulations here as well. In the present paper, we circumvent the problem of chiral fermions by using lattice input only for the Boseenhanced infrared degrees of freedom, treating fermions within a (resummed) weak-coupling expansion.
Now, because of dimensional reasons, we may rewrite eq. (3.1) as We note from the Euclidean path integral representation, viz.
that the bare pressure obeys ( where we wrote ν 2 B = Z m ν 2 (μ). 3 It is now straightforward to obtain the following relation: Here eq. (3.6) has been replaced by its renormalized version. It can be observed from eq. (3.7) that three ingredients are needed: the determination of explicit logarithms appearing inp R ("breaking of scale invariance by quantum corrections"); the temperature evolution of the Higgs condensate ("explicit breaking of scale invariance"); as well as the vacuum term needed for renormalization. We discuss the first of these in the next section, the condensate in sec. 5, and the vacuum term in sec. 6. The results are collected together and evaluated numerically in sec. 7.

Effects from scale violation by quantum corrections
The first ingredient needed in eq. (3.7) are logarithms ofμ/T that are induced by loop corrections. There are two ways to extract ∂p R ∂ ln(μ/T ) from ref. [18]: either by reading logarithms from the explicit expressions for the various coefficients given there, or by deducing them from the running of the couplings. Making use of the notation in eq. (2.5), withp ≡ p/T 4 , the partp G contributes at O(g 6 ). Fromp E +p M we need terms up to and including O(g 5 ). The expansion readŝ 3 The renormalization factor reads Zm = 1 + 3 (4π) 2 ǫ 2λ + h 2 t − 1 4 (g 2 1 + 3g 2 2 ) + O(g 4 ) but is not separately needed, because it always appears in the combination Zm φ † φ .
where the leading-order Debye masses read (n G = 3 denotes the number of generations) All coefficients in eq. (4.1) apart from α Eνν are scale independent; its value differs from that in ref. [18] because of our different renormalization condition: where d F , n S are specified in eq. (B.1). Given that the pressure as a whole is scale independent, running from the couplings must cancel against explicit logarithms; therefore, The runnings can be read from eqs. (7)-(12) of ref. [18].
Putting together, we obtain This expression is renormalization group (RG) invariant up to O(g 6 ). The numerically largest corrections originate from terms involving the strong gauge coupling g 2 3 .

Outline
The Higgs condensate can be obtained from eq. (3.6). For a non-perturbative evaluation, we make use of the lattice simulations in refs. [21,6]. This means that we split the pressure into contributions from various momentum scales like in eq. (2.5): Given that Z m φ † φ is multiplied by ν 2 in eq. (3.7) and that we assume That said, it turns out that computations going beyond those in refs. [18,19] are needed.

Perturbative contributions
The first term of eq. (5.1), the contribution from the "hard" scales k ∼ πT , can be directly extracted from the results of ref. [18]: Inserting coefficients from ref. [18] and from eq. (4.3), this contains 1/ǫ-divergences in the terms proportional to g 2 1 , g 2 2 and ν 2 : The second term of eq. (5.1), the contribution from the "soft" scales k ∼ gT , originates at O(g 3 ) from the dependence of the effective mass parameters on ν 2 : The first two terms give contributions that can be extracted from [19]: the Debye mass parameters depend on ν 2 as where β ′ E1 and β E1 are as in eq. (4.2). Differentiating the term on the second line of eq. (4.1) and inserting the values of the coefficients from ref. [19]  The last term of eq. (5.4) also contributes, but this contribution cannot be extracted from ref. [19]. The reason is that if we consider p M without a derivative, then contributions involving m 2 3 are of O(g 6 ) for |m 2 3 | < ∼ g 3 T 2 /π. However, the derivative of m 2 3 is larger than m 2 3 itself: Therefore terms that are of higher order in p M do contribute to the trace anomaly.
The computation of ∂p M /∂m 2 3 is simplified by the fact that since the softest physics is captured by a study of p G , the Higgs and gauge fields can be treated as massless in the computation of the matching coefficient p M . Then most diagrams, in particular all diagrams which do not contain at least one adjoint scalar propagator, vanish in dimensional regularization. The leading non-zero contributions, which are the only terms needed at O(g 3 ), are those given in fig. 1. Reducing the diagrams to master integrals and evaluating them with standard techniques (cf. e.g. ref. [22]), we obtain

Non-perturbative contribution
The third term of eq. (5.1), the contribution from the "ultrasoft" scales k ∼ g 2 T /π, can be expressed as wherem 2 3 is the Higgs mass parameter within MSM and the error comes from partial derivatives with respect to the other effective couplings of MSM. The dependence ofm 2 3 on ν 2 is known up to O(g 2 ) [13], with the leading term reading ∂m 2 3 /∂ν 2 = −1. The condensate can be expressed in dimensional regularization as where the argument of φ † φ 3d indicates the renormalization scale used within MSM. As discussed below eq. (2.3) the condensate is parametrically < ∼ gT 2 in our power counting, and therefore we keep the correction of O(g 2 ) in the coefficient of φ † φ 3d , even though we do not keep it in the ultraviolet part where it is unambiguously of O(g 4 ). For future reference we express the ultraviolet part in terms of the couplings of the full theory, by inserting 4 Thereby the ultrasoft contribution from k ∼ g 2 T /π becomes (5.14) The condensate φ † φ 3d (g 2 M2 ) can be measured non-perturbatively by subtracting proper counterterms [24] from lattice measurements extrapolated to the infinite-volume limit, and extrapolating subsequently to the continuum limit. A continuum extrapolation has only been carried out at an unphysically small Higgs mass [21] but cutoff effects have been seen to be modest, as long as we are not in the broken phase. Therefore we make use of lattice results at the physical Higgs mass [6] only for −0.1 < ∼ y < ∼ 0.2 in terms of the parameters in eq. (A.2). In order to extrapolate to higher or lower temperatures, perturbative expressions are employed; their values are discussed in appendix A. The procedure is illustrated in fig. 2.

Combined expression
Summing together the contributions from eqs. (5.3), (5.6), (5.7) and (5.14), most of the 1/ǫ-divergences cancel, and we get 4 For g 2 M2 corrections are known up to 2-loop order [23]. The other parameters of MSM read:  . The lattice data corresponds to a fixed lattice spacing a, which induces errors particularly at low temperatures [21]. For the results of sec. 7 the phenomenological interpolation indicated with the thin dashed line is employed.
where the finite part reads Apart from the last term the result isμ-independent up to O(g 6 ). Thisμ-dependence as well as the divergence in eq. (5.15) cancel against terms from the vacuum subtraction, cf. sec. 6.

Vacuum subtraction and renormalization
Suppose that we compute the bare vacuum pressure, p 0B , in a perturbative loop expansion: . In a gauge-fixed computation, the result is a function of the Higgs expectation value, which can likewise be determined order by order: v = ℓ v (ℓ) . Even though v is gauge-dependent, p 0B is gauge-independent order by order in perturbation theory. Inserting v into p 0B , we can re-expand the result as where we made use of p The terms not shown are of O(g 6 ) in our power counting. Therefore, it is sufficient for our purposes to compute p 0B up to 1-loop level and insert the tree-level Higgs vacuum expectation value v 2 = ν 2 /λ into the expression.
Dropping the superscripts, we get The counterterms read whereas at the minimum the masses take the values Most divergences cancel, and the vacuum part of the result becomes In practice, it is not useful to evaluate eq. (6.7) directly, because our renormalization scale will beμ ∼ πT , which would introduce large logarithms if inserted into eq. (6.7). Rather, all vacuum parameters are first evaluated at a scaleμ 0 ≡ m Z , and then evolved fromμ 0 to the thermal scale through RG equations. As an illustration, the 1-loop RG equations for the two dimensionful parameters appearing in our analysis read from where we get ∆ 3 (T ;μ) = 4p 0R /T 4 . At the scaleμ =μ 0 , the running couplings are expressed in terms of physical parameters through 1-loop relations specified e.g. in ref. [13].
(For many parameters a higher loop order could be extracted from more recent literature, but such corrections are smaller than thermal uncertainties, so for the sake of simplicity and reproducibility we make use of explicit 1-loop expressions.)

Phenomenological results
Summing together the terms from eqs. (4.6), (5.15), (6.6) as dictated by eq. (3.7), all 1/ǫdivergences cancel, and we obtain Because of a cancellation between ∆ 2 and ∆ 3 the result is formally independent of the renormalization scaleμ, even though in practice a residualμ-dependence is left over, as a reflection of unknown higher order corrections. We writeμ = απT , and vary α in the range α ∈ (0.5...2.0) in order to get one impression on the corresponding uncertainty. The result of eq. (7.1) is accurate up to and including the order O(g 5 ). It is well known from studies of the pressure of QCD, however, that certain odd orders show anomalously poor convergence. In particular, whereas the O(g 2 ) correction to the pressure provides for a reasonable approximation, the O(g 3 ) correction is far off. For ∆, the order O(g 4 ) is related to the O(g 2 ) correction to the pressure, and the order O(g 5 ) to the O(g 3 ) correction (cf. sec. 4). For this reason, the numerically most accurate estimate for ∆ can probably be obtained by restricting to O(g 4 ). A result corresponding to this accuracy is shown for the observable of eq. (7.1) in fig. 3. (For the pressure we display also the O(g 5 ) result in fig. 4(left).) In order to obtain other thermodynamic functions, the boundary value needed for eq. (3.2) should be fixed. We do this on the low-temperature side, making use of the results of ref. [25] 5 and thereby setting p(T 0 )/T 4 0 ≃ 10.91 at T 0 = 100 GeV. On the high-temperature side, we expect to match to the results of ref. [18], after changing the overall renormalization condition to our eq. (2.6) and by making the changes listed in appendix B. Like for ∆, we make use of the result of O(g 4 ), with the negative O(g 5 ) QCD contribution expected to lead to an underestimate [26,27]. (We have also experimented with an approximate O(g 6 ) QCD contribution as estimated in ref. [25], finding a result which is in between the O(g 4 ) and O(g 5 ) ones, confirming that the O(g 5 ) result is most likely on the low side.) After choosing an initial condition, other thermodynamic functions are obtained as follows: the pressure p/T 4 from eq. (3.2); the energy density from e/T 4 = ∆ + 3p/T 4 ; the entropy density s = p ′ from s/T 3 = ∆ + 4p/T 4 ; the heat capacity c = e ′ from c/T 3 = T ∆ ′ + 7∆ + 12p/T 4 ; the equation-of-state parameter from w = p/e = 1/(3 + ∆T 4 /p); and the speed of sound squared from c 2 s = p ′ /e ′ = s/c. Some of these functions are conveniently parametrized through Results for all of these functions are shown in fig. 4. (The heat capacity and the speed of sound squared are the most "difficult" quantities, because they require taking a numerical temperature derivative from ∆.) It may be wondered why our new results do not match exactly the previous ones for high temperatures [18], even though both have been computed up to the same order in g. The main reason is that they involve different sets of higher-order corrections. In particular, eq. (7.1) has been evaluated "strictly" to O(g 4 ) or O(g 5 ), because this has led to fairly simple expressions. In contrast, the high-temperature result is more cumbersome, including more terms (with non-zero corrections of O(1), O(g 2 ), O(g 3 )) and complicated "soft" contributions, because the thermal mass parameters associated with the Higgs and gauge fields are of the same order in this regime. For evaluating such contributions the experience from QCD suggests that it is not worth re-expanding them in terms of the original couplings but rather to evaluate p M "as is". We have followed separate procedures for the two regimes, because their difference permits for us to get another impression on the magnitude of unknown higherorder corrections. Indeed, the mismatches in fig. 4 are on the percent level which we expect to represent the theoretical uncertainty of our analysis.
With the help of the functions in eq. (7.2), the relationship of time and temperature in the Early Universe (assuming a flat geometry) can be expressed as It is seen that a peak in heat capacity (i.e. i eff ), visible in fig. 4(middle), leads to a short period of slower temperature change [28], and correspondingly a mildly increased abundance of produced particles if a particle production process is under way, or a reduced density of weakly interacting relics disappearing through a co-annihilation process.

Conclusions
Drawing on existing multiloop computations [18] and lattice simulations within a dimensionally reduced effective theory [6], and complementing these through new 1-loop, 2-loop and 3-loop computations needed for determining the "trace anomaly" of the electroweak theory, we have estimated the basic thermodynamic functions that play a role for Standard Model thermodynamics at temperatures between 100 GeV and 300 GeV (cf. fig. 4). These results can be matched, for most observables with modest (∼ 1%) discrepancies, to perturbative computations at low [25] or high [18] temperatures.
One finding of our study is that despite the high temperatures considered, radiative corrections are larger than might naively be anticipated. Larger corrections also imply that uncertainties are less well under control. On the low-temperature side, we believe that the results of ref. [25] do contain uncertainties because the analysis of the electroweak sector was based on a low loop order and because QCD corrections, whose convergence is slow at finite temperatures, are substantial. On the high-temperature side, the QCD corrections continue to be an issue (cf. fig. 4(left)). As an example of the physical significance of these uncertainties, it may be noted from fig. 4 that the effective numbers of degrees of freedom remain below the canonical value 106.75 in the whole temperature range considered, and even decrease modestly as the temperature increases above 300 GeV. Even though there could be valid physics reasons for the decrease (such as that the effective Higgs mass parameter is very small across the crossover but then increases again), a similar tendency also originates from the O(g 4 ) QCD contribution and is then an artifact of the truncation. (The O(g 5 ) correction decreases the effective numbers of degrees of freedom further but turns the results into slowly increasing functions). Another peculiar feature of the O(g 4 ) result is that the speed of sound squared c 2 s is ∼ 0.1% above 1/3 at T > 300 GeV, however this is again reversed by the O(g 5 ) correction. In general, the estimated theoretical uncertainty of our results is on the 1% level, and no conclusions should be drawn from features finer than this.
In comparison with ref. [25], we find results for the effective numbers of degrees of freedom that are about 1% lower at T > 300 GeV. The reason is mostly due to negative electroweak radiative corrections proportional to g 2 2 and h 2 t , which were omitted in ref. [25]. Apart from radiative corrections, another important ingredient in our results is the 3d condensate φ † φ 3d (cf. eq. (5.14)), measured non-perturbatively on the lattice and subsequently converted to the MS scheme. It is interesting to note that apart from the known strong correlation of φ † φ 3d with the anomalous rate of baryon plus lepton number violation (cf. refs. [30,6] and references therein), the same quantity also plays a role for other cosmologically relevant observables, such as the production rate of non-relativistic right-handed neutrinos [31] or the relationship between lepton and baryon number densities [32]. Therefore, it seems well motivated to improve on the existing measurements [6] by including the U(1) subgroup and by taking the continuum limit for a physical Higgs mass. In addition, it would be helpful to measure the susceptibility related to this condensate, so that taking a numerical temperature derivative, needed for estimating the heat capacity or the speed of sound squared, could be avoided.
Our interpolations of all thermodynamic functions shown in fig. 4 can be downloaded, for a wide temperature interval, from www.laine.itp.unibe.ch/eos15/. (In these results we switch from one regime to another with a temperature gap of 15-20 GeV in between.) Despite the remaining uncertainties, we hope that these results turn out to be helpful as a background equation of state in Dark Matter or Leptogenesis computations operating in this temperature range.
In the broken symmetry phase, i.e. form 2 3 ≪ 0, an explicit loopwise result is available up to 2-loop level, however it shows poor convergence [29]. Therefore we make use of a numerically determined value (referred to as the "Coleman-Weinberg (CW) method") which includes a subset of higher-order corrections and has been tested against lattice simulations in refs. [29,21]. The result is only available in the approximation g M1 = 0 but the corrections originating from g M1 are expected to be small. A comparison with lattice data is shown in fig. 2 (long-dashed line), with the low-temperature deviation expected to be reduced by a continuum extrapolation [21].