Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model

TL;DR

This study refines the estimation of primordial gravitational wave spectra by incorporating detailed thermodynamic effects from the Standard Model, revealing subtle damping features and higher amplitudes relevant for future detection efforts.

Contribution

It provides a comprehensive reevaluation of the effective degrees of freedom in the Standard Model, including interaction effects, and assesses their impact on gravitational wave spectra.

Findings

Additional damping features due to free-streaming particles.

Smoother spectrum across the QCD phase transition.

Primordial gravitational wave amplitude is about 1% larger than previous estimates.

Abstract

In this paper, we revisit the estimation of the spectrum of primordial gravitational waves originated from inflation, particularly focusing on the effect of thermodynamics in the Standard Model of particle physics. By collecting recent results of perturbative and non-perturbative analysis of thermodynamic quantities in the Standard Model, we obtain the effective degrees of freedom including the corrections due to non-trivial interaction properties of particles in the Standard Model for a wide temperature interval. The impact of such corrections on the spectrum of primordial gravitational waves as well as the damping effect due to free-streaming particles is investigated by numerically solving the evolution equation of tensor perturbations in the expanding universe. It is shown that the reevaluation of the effects of free-streaming photons and neutrinos gives rise to some additional damping features overlooked in previous studies. We also observe that the continuous nature of the QCD crossover results in a smooth spectrum for modes that reenter the horizon at around the epoch of the QCD phase transition. Furthermore, we explicitly show that the values of the effective degrees of freedom remain smaller than the commonly used value 106.75 even at temperature much higher than the critical temperature of the electroweak crossover, and that the amplitude of primordial gravitational waves at a frequency range relevant to direct detection experiments becomes $\mathcal{O}(1)\,\%$ larger than previous estimates that do not include such corrections. This effect can be relevant to future high-sensitivity gravitational wave experiments such as ultimate DECIGO. Our results on the temperature evolution of the effective degrees of freedom are made available as tabulated data and fitting functions, which can also be used in the analysis of other cosmological relics.

Figures

Tables

Table 1. Table 1: Coefficients for the fitting functions for 120 ​ MeV ≤ T ≤ 10 16 ​ GeV 120 MeV 𝑇 superscript 10 16 GeV 120\,\mathrm{MeV}\leq T\leq 10^{16}\,\mathrm{GeV} .

$i$	$a_i$	$b_i$	$c_i$	$d_i$
0	1	1.43382E$-$02	1	7.07388E+01
1	1.11724E$+$00	1.37559E$-$02	6.07869E$-$01	9.18011E+01
2	3.12672E$-$01	2.92108E$-$03	$-$1.54485E$-$01	3.31892E+01
3	$-$4.68049E$-$02	$-$5.38533E$-$04	$-$2.24034E$-$01	$-$1.39779E+00
4	$-$2.65004E$-$02	$-$1.62496E$-$04	$-$2.82147E$-$02	$-$1.52558E+00
5	$-$1.19760E$-$03	$-$2.87906E$-$05	2.90620E$-$02	$-$1.97857E$-$02
6	1.82812E$-$04	$-$3.84278E$-$06	6.86778E$-$03	$-$1.60146E$-$01
7	1.36436E$-$04	2.78776E$-$06	$-$1.00005E$-$03	8.22615E$-$05
8	8.55051E$-$05	7.40342E$-$07	$-$1.69104E$-$04	2.02651E$-$02
9	1.22840E$-$05	1.17210E$-$07	1.06301E$-$05	$-$1.82134E$-$05
10	3.82259E$-$07	3.72499E$-$09	1.69528E$-$06	7.83943E$-$05
11	$-$6.87035E$-$09	$-$6.74107E$-$11	$-$9.33311E$-$08	7.13518E$-$05

Equations

$$w(T)\equiv \frac{p(T)}{\rho (T)},$$

$$g_{\ast \rho }(T)\equiv \frac{\rho (T)}{\left[\frac{{\pi }^2{T}^4}{30}\right]},g_{\ast s}(T)\equiv \frac{s(T)}{\left[\frac{2{\pi }^2{T}^3}{45}\right]},$$

$$d{s}^2=-d{t}^2+a^2(t)({\delta }_{ij}+h_{ij})d{x}^{i}d{x}^j,$$

$${\ddot{h}}_{ij}+3H{\dot{h}}_{ij}-\frac{{\nabla }^2}{a^2}{h}_{ij}=16\pi G{\Pi }_{ij}^{\mathrm{TT}},$$

$$a^2{\Pi }_{ij}=T_{ij}-p{g}_{ij}.$$

$$h_{ij}(t,\mathbf{x})=\sum _{\lambda }\int \frac{d^{3}k}{(2\pi {)}^3}{h}^{\lambda }(t,\mathbf{k}){ϵ}_{ij}^{\lambda }(\mathbf{k})e^{i\mathbf{k}\cdot \mathbf{x}},$$

$${\ddot{h}}_{\mathbf{k}}^{\lambda }+3H{\dot{h}}_{\mathbf{k}}^{\lambda }+\frac{k^2}{a^2}{h}_{\mathbf{k}}^{\lambda }=16\pi G{\Pi }_{\mathbf{k}}^{\lambda },$$

$$h_{\mathbf{k}}^{\lambda }=h_{\mathbf{k},\mathrm{prim}}^{\lambda }\chi (t,k),$$

$${\rho }_{\mathrm{gw}}(t)=\frac{1}{32\pi G}⟨{\dot{h}}_{ij}(t,\mathbf{x}){\dot{h}}_{ij}(t,\mathbf{x}){⟩}_{\mathrm{av}}=\frac{1}{32\pi G}\sum _{\lambda }\int \frac{d^{3}k}{(2\pi {)}^3}2|{\dot{h}}_{\mathbf{k}}^{\lambda }{|}^2,$$

$${\Omega }_{\mathrm{gw}}(t,k)=\frac{1}{{\rho }_{\mathrm{crit}}(t)}\frac{d{\rho }_{\mathrm{gw}}(t,k)}{d\ln k},$$

$${\Omega }_{\mathrm{gw}}(\tau ,k)=\frac{1}{12a^2(\tau )H^2(\tau )}{\mathcal{P}}_T(k){\left[{\chi }^{\prime }(\tau ,k)\right]}^2,$$

$${\mathcal{P}}_T(k)\equiv \frac{k^3}{{\pi }^2}\sum _{\lambda }|h_{\mathbf{k},\mathrm{prim}}^{\lambda }{|}^2={\frac{2H^2}{{\pi }^2{M}_{\mathrm{Pl}}^2}|}_{k=aH},$$

$${\left[{\chi }^{\prime }(\tau ,k)\right]}^2\approx \frac{k^2}{2}{\left(\frac{a({\tau }_{\mathrm{hc}})}{a(\tau )}\right)}^2=\frac{a^4({\tau }_{\mathrm{hc}})H^2({\tau }_{\mathrm{hc}})}{2a^2(\tau )},$$

$$\frac{H_{\mathrm{hc}}^2}{H_0^2}=\frac{\frac{{\pi }^2}{30}{g}_{\ast \rho ,\mathrm{hc}}{T}_{\mathrm{hc}}^4}{{\rho }_{\mathrm{crit},0}}={\Omega }_{\gamma }\left(\frac{g_{\ast \rho ,\mathrm{hc}}}{2}\right){\left(\frac{g_{\ast s,\mathrm{hc}}}{g_{\ast s,\mathrm{fin}}}\right)}^{-\frac{4}{3}}{\left(\frac{a_0}{a_{\mathrm{hc}}}\right)}^4,$$

$${\Omega }_{\mathrm{gw}}({\tau }_0,k)\approx \frac{1}{24}{\Omega }_{\gamma }\left(\frac{g_{\ast \rho ,\mathrm{hc}}}{2}\right){\left(\frac{g_{\ast s,\mathrm{hc}}}{g_{\ast s,\mathrm{fin}}}\right)}^{-\frac{4}{3}}{\mathcal{P}}_T(k).$$

$${\Omega }_{\mathrm{gw}}({\tau }_0,k)h^2$$

$$\displaystyle\approx 1.29\times 10^{-17}\left(\frac{g_{*s,\mathrm{fin}}}{3.931}\right)^{\frac{4}{3}}\left(\frac{g_{*\rho,{\rm hc}}}{106.75}\right)\left(\frac{g_{*s,{\rm hc}}}{106.75}\right)^{-\frac{4}{3}}\bigg{(}\frac{V^{1/4}_{\rm inf}}{10^{16}\,\mathrm{GeV}}\bigg{)}^{4}.$$

$$f=\frac{k}{2\pi a_0}=\frac{H_{\mathrm{hc}}}{2\pi }\frac{a_{\mathrm{hc}}}{a_0}\approx 2.65\mathrm{Hz}{\left(\frac{g_{\ast s,\mathrm{fin}}}{3.931}\right)}^{\frac{1}{3}}{\left(\frac{g_{\ast \rho ,\mathrm{hc}}}{106.75}\right)}^{\frac{1}{2}}{\left(\frac{g_{\ast s,\mathrm{hc}}}{106.75}\right)}^{-\frac{1}{3}}\left(\frac{T_{\mathrm{hc}}}{10^8\mathrm{GeV}}\right).$$

$${\Pi }_{\mathbf{k}}^{\lambda }=-4{\rho }_{\nu }(\tau ){\int }_{{\tau }_{\nu \mathrm{dec}}}^{\tau }d{\tau }^{\prime }\frac{j_2[k(\tau -{\tau }^{\prime })]}{k^2(\tau -{\tau }^{\prime }{)}^2}\frac{d{h}_{\mathbf{k}}^{\lambda }({\tau }^{\prime })}{d{\tau }^{\prime }},$$

$${\Pi }_{\mathbf{k}}^{\lambda }=-4{\rho }_{\gamma }(\tau ){\int }_{{\tau }_{\mathrm{ls}}}^{\tau }d{\tau }^{\prime }\frac{j_2[k(\tau -{\tau }^{\prime })]}{k^2(\tau -{\tau }^{\prime }{)}^2}\frac{d{h}_{\mathbf{k}}^{\lambda }({\tau }^{\prime })}{d{\tau }^{\prime }},$$

$$\frac{d\ln {\Omega }_{\mathrm{gw}}(\tau ,k)}{d\ln a}=\frac{d\ln {\tilde{\rho }}_{\mathrm{gw}}(\tau ,k)}{d\ln a}-\frac{d\ln {\rho }_{\mathrm{crit}}(\tau )}{d\ln a},$$

$$\frac{d\ln {\rho }_{\mathrm{crit}}(\tau )}{d\ln a}=-3\left(1+w(a)\right).$$

$$\frac{d\ln {\tilde{\rho }}_{\mathrm{gw}}(\tau ,k)}{d\ln a}\approx -4.$$

$${\Omega }_{\mathrm{gw}}(\tau ,k)\approx \exp \left[{\int }_{a_{\mathrm{hc}}}^a\left(3w(a)-1\right)d\ln a\right]{\Omega }_{\mathrm{gw}}({\tau }_{\mathrm{hc}},k).$$

$$\frac{{\Omega }_{\mathrm{gw}}({\tau }_0,k_2)}{{\Omega }_{\mathrm{gw}}({\tau }_0,k_1)}$$

$$=\left(\frac{g_{\ast \rho ,\mathrm{hc},2}}{g_{\ast \rho ,\mathrm{hc},1}}\right){\left(\frac{g_{\ast s,\mathrm{hc},2}}{g_{\ast s,\mathrm{hc},1}}\right)}^{-\frac{4}{3}}\left[\frac{{\Omega }_{\mathrm{gw}}({\tau }_{\mathrm{hc},2},k_2)}{{\Omega }_{\mathrm{gw}}({\tau }_{\mathrm{hc},1},k_1)}\right].$$

$${\ddot{h}}_{ij}+\left(3H+16\pi G\eta \right){\dot{h}}_{ij}-\frac{{\nabla }^2}{a^2}{h}_{ij}=0.$$

$$p\to p-3H\zeta .$$

$$w-\frac{1}{3}=-\frac{3H\zeta }{\rho }=-\frac{8\pi G\zeta }{H}.$$

$$\Delta (T)\equiv \frac{\rho (T)-3p(T)}{T^4}=T\frac{d}{dT}\left\{\frac{p(T)}{T^4}\right\},$$

$$\rho (T)=T^4\left[\Delta (T)+\frac{3p(T)}{T^4}\right],s(T)=T^3\left[\Delta (T)+\frac{4p(T)}{T^4}\right],$$

Full text

1. 1 Introduction

2. 2 Inflationary gravitational waves

3. 2.1 Energy density of gravitational waves

4. 2.2 Damping effects

5. 3 Equation of state in the Standard Model

6. 3.1 Photons and leptons

7. 3.2 Strongly interacting particles

8. 3.3 Electroweak crossover

9. 3.4 Full equation of state and its uncertainty

10. 4 Spectrum of the gravitational waves

11. 5 Conclusions and discussion

12. A Supplementary material

13. B Notes on thermodynamic quantities

14. C Fitting functions for effective degrees of freedom

15. D Estimation of the sensitivity of gravitational wave detectors

1 Introduction

The recent detections of gravitational waves (GWs) in the Laser Interferometer Gravitational-Wave Observatory (LIGO) and Virgo [1, 2, 3, 4, 5, 6] opened up new possibilities to investigate various astrophysical phenomena which cannot be probed by the conventional methods based on electromagnetic waves. Although the ongoing detectors are only sensitive to strong transient events such as merging black holes, future GW experiments are expected to detect much weaker signatures of GWs produced in the early universe [see, e.g. Refs. [7, 8, 9] for reviews]. Several space-borne interferometers such as Laser Interferometer Space Antenna (LISA) [10] and Deci-hertz Interferometer Gravitational Wave Observatory (DECIGO) [11, 12] are planned to be launched in the future with the aim of detecting the primordial GW background. Given these perspectives, it will become more important to improve the precision of theoretical calculations of the primordial GWs in order to distinguish various models describing the history of the early universe.

The existence of the primordial GW background is one of the most crucial predictions of the inflationary scenario of the early universe [13, 14]. In order to estimate the spectrum of the inflationary GWs, basically we need to know two factors: One is the power spectrum of primordial tensor perturbations generated during inflation, and the other is the expansion rate of the universe from the end of inflation until today. The former defines the initial magnitude of the GW signature, and it is directly associated with the detailed properties of inflationary models [15, 16, 17]. On the other hand, the latter describes how the density of the primordial GWs has been diluted in subsequent stages of the cosmic expansion. In other words, the primordial GWs “know” how the universe has evolved after inflation. This fact implies that there is a possibility to obtain detailed information about the early history of the universe as well as the inflationary models by observing the spectrum of the primordial GWs.

It is known that in the standard slow-roll inflationary models the spectrum of the inflationary GW background is almost flat for a broad range of frequencies, since the amplitude of GWs gets a common dilution factor if the corresponding mode reenters the horizon during the radiation dominated era. The deviation from the flat spectrum can be associated with either the tilt of the primordial power spectrum or the fact that the hot plasma produced after inflation does not exactly behave like an ideal gas of massless particles, which can be caused by non-trivial interaction properties of elementary particles. In order to evaluate the latter effect correctly, we need to know the concrete theory of particle physics. The Standard Model (SM), which formulates all fundamental interactions of quarks and leptons in terms of electromagnetic, weak, and strong forces, can be regarded as a “benchmark” model to describe the properties of the primordial plasma. The non-trivial interactions of the SM plasma can affect the spectrum of the primordial GWs, which has to be quantitatively taken into account.

The intervening effects which alter the nature of the spectrum of GWs during their propagation are represented as a transfer function, and we have to follow several steps in order to calculate it correctly. First, Schwarz [18] explicitly pointed out that the deviation of the behavior of fluid components of the universe from that of a simple ideal gas of relativistic particles is imprinted in the spectrum of GWs, and that the observation of the spectrum of GWs can be used to probe the equation of state of the early universe [see also Ref. [19]]. Second, Weinberg [20] considered the effects of neutrinos, which decouple from the thermal bath and start to free-stream when the universe becomes cooled below $\sim\mathcal{O}(1)\,\mathrm{MeV}$. It was shown that such neutrinos damp the amplitude of GWs by 35.5 % in the frequency range between $\sim 10^{-16}\,\mathrm{Hz}$ and $\sim 10^{-10}\,\mathrm{Hz}$.111It appears that the damping effect due to the free-streaming particles was already pointed out in earlier works [see, e.g. Refs. [21, 22]]. Weinberg’s paper was the first one to point out that this effect might indeed be relevant to observations. This damping effect due to the free-streaming particles was further investigated by several authors [23, 24, 25, 26, 27, 28]. A more precise computational method was developed by Watanabe and Komatsu [29], and it was pointed out that the change of the effective degrees of freedom in the SM elementary particles modifies the spectrum of GWs. In addition to such an effect, Boyle and Steinhardt [30] explored general possibilities that the equation of state might be modified by some non-trivial features of interacting particles. Finally, Kuroyanagi, Chiba, and Sugiyama [31] performed full numerical calculations of the inflationary GW spectrum over broad frequency ranges by taking into account the scalar field dynamics during inflation, the process during the reheating after inflation, the damping due to the free-streaming neutrinos, and the change of the effective degrees of freedom.

It should be noted that the current estimate of the transfer function is still incomplete in spite of several developments described above. The key ingredients considered here are the equation of state of the primordial plasma,

[TABLE]

and the effective degrees of freedom,

[TABLE]

where $p(T)$, $\rho(T)$, and $s(T)$ are the pressure, energy density, and entropy density of the primordial plasma at temperature $T$. In most previous studies of primordial GWs (and also in those of other cosmological relics), the values of $g_{*\rho}(T)$ and $g_{*s}(T)$ were estimated based on the ideal gas approximation, and the equation of state was assumed to be $1/3$ in the radiation dominated universe. However, these approximations do not always hold once we take account of the effect of interactions of elementary particles. It is reasonable to expect that $g_{*\rho}$ and $g_{*s}$ as well as $w$ can deviate from commonly used values due to the non-trivial properties of particle interactions even in the SM of particle physics, and in principle such corrections can be imprinted on the spectrum of primordial GWs. The main question addressed in this paper is to what extent the shape of the spectrum of primordial GWs can be modified if we systematically include the corrections on the thermodynamic quantities of the primordial plasma arising from particle interactions in the SM.

In order to estimate the thermodynamic properties of interacting particles precisely, we must use state-of-the-art methods of quantum field theory. In the literature, there has been a lot of work on the estimation of thermodynamic quantities in the SM based on the finite temperature field theory, and we aim to collect all the results to reconstruct the relevant quantities such as the equation of state parameter and the effective degrees of freedom at arbitrary temperature. We show that the improvement in estimation of thermodynamic properties of high temperature plasma in the SM has several impacts on the calculation of the spectrum of primordial GWs, and some of them are relevant to future high-sensitivity GW experiments. As a by-product, we also obtain some detailed information on the temperature evolution of $g_{*\rho}(T)$ and $g_{*s}(T)$ and their uncertainty for a wide temperature interval, which is made available as tabulated data and fitting functions.

The organization of this paper is as follows. In Sec. 2, we describe some basic properties of inflationary GWs and discuss several damping effects that can take place during their evolution after inflation. In Sec. 3, we estimate the effective degrees of freedom $g_{*\rho}(T)$ and $g_{*s}(T)$ and the equation of state parameter $w(T)$ at arbitrary temperature by collecting recent developments including the analysis of the neutrino decoupling, quantum chromodynamics (QCD) phase transition, and electroweak phase transition. After obtaining the updated values of the effective degrees of freedom in the SM, in Sec. 4 we calculate the spectrum of GWs based on them and highlight difference from the results obtained in previous studies. Sec. 5 is devoted to conclusions and discussion. Some technical details are presented in appendixes. Appendix A deals with our tabulated data of $g_{*\rho}(T)$ and $g_{*s}(T)$ obtained in Sec. 3 and that of the transfer function of GWs obtained in Sec. 4. Appendix B is devoted to a short review of thermodynamic quantities used in the main text. In appendix C, we give fitting functions for $g_{*\rho}(T)$ and $g_{*s}(T)$ that remain consistent with the actual data obtained in Sec. 3 within a range of uncertainty. In appendix D, we describe some details about sensitivity of ultimate DECIGO mentioned in Sec. 4.

2 Inflationary gravitational waves

In this section, we describe how to estimate the spectrum of the inflationary GW background and enumerate several post-inflationary events that can affect its shape. First, we consider the evolution of GWs in the expanding universe and derive the relation among the energy density of GWs, the primordial tensor power spectrum and the transfer function in Sec. 2.1. After that, in Sec. 2.2 we discuss various damping effects of the primordial tensor perturbations encoded in the transfer function.

2.1 Energy density of gravitational waves

We work in a spatially flat Friedmann-Robertson-Walker (FRW) background with a metric,

[TABLE]

where GWs are represented by spatial metric perturbations ($|h_{ij}|\ll 1$) that satisfy the transverse-traceless conditions, $h^{i}_{i}=\partial^{i}h_{ij}=0$. The evolution of GWs is described by the linearized Einstein equation,

[TABLE]

where a dot represents a derivative with respect to cosmic time $t$, and $G$ is the Newton’s gravitational constant. $\Pi^{\rm TT}_{ij}$ are the transverse-traceless part of the anisotropic stress $\Pi_{ij}$, which are defined in terms of the spatial components of the stress-energy tensor $T_{ij}$, those of the metric tensor $g_{ij}$, and the background homogeneous pressure $p$,

[TABLE]

It is convenient to decompose $h_{ij}$ into their Fourier modes,

[TABLE]

where $\lambda=+,\times$ specifies two independent polarization states, and $\epsilon^{\lambda}_{ij}({\bf k})$ are the spin-2 polarization tensors satisfying the normalization conditions $\sum_{ij}\epsilon^{\lambda}_{ij}(\epsilon^{\lambda^{\prime}}_{ij})^{*}=2\delta^{\lambda\lambda^{\prime}}$. In terms of the Fourier modes, Eq. (2.2) can be rewritten as

[TABLE]

where $h^{\lambda}_{\bf k}(t)\equiv h^{\lambda}(t,{\bf k})$, and $\Pi^{\lambda}_{\bf k}$ are corresponding Fourier components of the source term $\Pi^{\rm TT}_{ij}$.

The GW background is originated from the quantum fluctuations during inflation, which are assumed to be classicalized after the corresponding modes cross outside the horizon ($k<aH$). For the modes far outside the horizon ($k\ll aH$), we can ignore the source term and the third term of the left-hand side of Eq. (2.5), which implies that $h^{\lambda}_{\bf k}\propto\text{const}$. At some point after inflation, these modes reenter the horizon ($k>aH$) and start to evolve according to Eq. (2.5). Based on these facts, we write the solution of Eq. (2.5) as

[TABLE]

where $h^{\lambda}_{{\bf k},\mathrm{prim}}$ represents the amplitude of the primordial tensor perturbations, and $\chi(t,k)$ is the transfer function. The magnitude of the transfer function is normalized such that $\chi(t,k)\to 1$ for $k\ll aH$. If we ignore the source therm ($\Pi^{\lambda}_{\bf k}=0$), the evolution of the modes deep inside the horizon ($k\gg aH$) is well approximated by the WKB solution, $\chi(\tau,k)\propto a^{-1}e^{\pm ik\tau}$, where $\tau$ is conformal time defined by $d\tau=dt/a$.

The energy density of the relic GWs is given by [32]

[TABLE]

where $\langle\dots\rangle_{\rm av}$ denotes an spatial average. The spectrum of GWs is described in terms of the fraction of their energy density per logarithmic frequency interval,

[TABLE]

where $\rho_{\rm crit}=3H^{2}/8\pi G$ is the critical energy density of the universe. Substituting Eq. (2.6), we obtain

[TABLE]

where a prime represents a derivative with respect to conformal time $\tau$, i.e. ${}^{\prime}=\frac{d}{d\tau}=a\frac{d}{dt}$. Here we introduced the primordial tensor power spectrum $\mathcal{P}_{T}(k)$, which is determined by the Hubble parameter at the time when the corresponding mode crosses the horizon during inflation ($k=aH$),

[TABLE]

where $M_{\rm Pl}\simeq 2.435\times 10^{18}\mathrm{GeV}$ is the reduced Planck mass.

Let us estimate the typical amplitude of GWs at the present time $\tau=\tau_{0}$ for the modes that reenter the horizon during the radiation dominated era. From the fact that the transfer function can be approximated by the WKB solution ($\chi\propto a^{-1}e^{\pm ik\tau}$) after the modes reenter the horizon, and that it should approach to $\chi\to 1$ at early times, we can make the following replacement

[TABLE]

where $\tau_{\rm hc}$ represents the conformal time of the horizon crossing, and we used $k=a(\tau_{\rm hc})H(\tau_{\rm hc})$ in the second equality. The factor $1/2$ arises from the time average of the rapidly oscillating function. If the horizon crossing happens during the radiation dominated era, we can also use the following relation,

[TABLE]

where the subscript “[math]” and “${\rm hc}$” represent the quantity at the time $\tau_{0}$ and $\tau_{\rm hc}$, respectively, and $g_{*\rho}$ and $g_{*s}$ are the effective degrees of freedom for the energy density and for the entropy density, respectively. $\Omega_{\gamma}=\rho_{\gamma,0}/\rho_{\rm crit,0}$ is the fraction of the energy density of photons $\rho_{\gamma,0}=(\pi^{2}/15)T_{\gamma,0}^{4}$ at the present time. Using the measured value $T_{\gamma,0}=2.72548(57)\,\mathrm{K}$ [33] for current cosmic microwave background (CMB) temperature, we obtain $\Omega_{\gamma}h^{2}=2.4728(21)\times 10^{-5}$, where $h$ is the renormalized Hubble parameter, $H_{0}=100\,h\,\mathrm{km}\,\mathrm{sec}^{-1}\mathrm{Mpc}^{-1}$. In the last equality in Eq. (2.12), we used the conservation of entropy, $g_{*s}a^{3}T^{3}=\text{const}$. Note that we have introduced the notation $g_{*s,\mathrm{fin}}$ rather than $g_{*s,0}$ in Eq. (2.12), where $g_{*s,\mathrm{fin}}$ is the value of $g_{*s}$ evaluated after the neutrino decoupling. We use this notation in order to remind that the value of $g_{*s}$ is not directly evaluated at the present time, but indirectly estimated at the epoch of the neutrino decoupling. Actually, this procedure gives rise to some ambiguities as we discuss in Sec. 3.1.

Substituting Eqs. (2.11) and (2.12) into Eq. (2.9), we obtain

[TABLE]

We see that the spectrum shown in Eq. (2.13) is almost flat except for the (weak) scale dependences in $\mathcal{P}_{T}(k)$, $g_{*\rho,{\rm hc}}$ and $g_{*s,{\rm hc}}$. It should be noted that the amplitude of GWs for a given mode is related to the Hubble parameter $H_{\rm inf}$ or the energy scale of inflation $V^{1/4}_{\rm inf}$ when the corresponding mode exits the horizon,

[TABLE]

Furthermore, the frequency $f$ of GWs can be related to the temperature $T_{\rm hc}$ at which the corresponding mode reenters the horizon,

[TABLE]

We emphasize that Eq. (2.13) is derived based on the assumption that the transfer function is exactly given by the WKB solution after the horizon crossing. Strictly speaking, this is not the case if the equation of state parameter $w$ of the background fluid deviates from the value for pure radiation $1/3$ around the time of the horizon crossing. Accordingly, the amplitude of GWs can be smaller (or larger) than the right-hand side of Eq. (2.13) under certain conditions. We further elaborate on this point in Sec. 4.

2.2 Damping effects

Various phenomena occurring after inflation can affect the transfer function, which results in a non-trivial shape of the spectrum of GWs at the present time. One important effect is the collisionless damping caused by free-streaming neutrinos.222In this paper, we only consider the damping effects due to the SM neutrinos and photons. If there exist additional free-streaming particles, they can cause similar damping effects. Such possibilities are explicitly analyzed in Refs. [34, 35]. We also note that the damping effect due to cold dark matter is negligibly small since the fraction of its kinetic energy is small or its contribution to the energy density is small compared to the total energy density of the universe [36]. Neutrinos decouple from the thermal bath when the temperature drops below $\sim 2\,\mathrm{MeV}$, and the subsequent evolution of their distribution function is described by the collisionless Boltzmann equation. The existence of primordial GWs causes a tensor type perturbation in the neutrino distribution function, which contributes to the transverse-traceless part of the anisotropic stress and results in the energy flow from neutrinos to GWs, or vice versa.

At linear order in perturbation theory, the contribution to the anisotropic stress is simply given by [20, 29]

[TABLE]

where $\rho_{\nu}(\tau)$ is the homogeneous energy density of neutrinos, $\tau_{\nu{\rm dec}}$ is the conformal time of the neutrino decoupling, and $j_{n}(z)$ is the spherical Bessel function of the first kind. The evolution equation of GWs [Eq. (2.5)] becomes an integro-differential equation, which can be solved numerically. It was shown that the effect of the neutrino free-streaming leads to the damping of the amplitude of primordial GWs by 35.5 % in the frequency range $10^{-16}\,\mathrm{Hz}\lesssim f\lesssim 10^{-10}\,\mathrm{Hz}$. Note that this damping effect does not work for the modes that reenter the horizon before the neutrino decoupling ($f\gtrsim 10^{-10}\,\mathrm{Hz}$), since these modes are rapidly oscillating already at the time of the neutrino decoupling and there is no net energy conversion between neutrinos and GWs. This effect also becomes less important for the modes that reenter the horizon after the time of mater-radiation equality, since the energy density of neutrinos becomes small compared to the total energy density of the universe during the matter dominated era.

In principle, the damping of primordial GWs can also be caused by photons, which decouple from the thermal bath at $T\sim 3000\,\mathrm{K}$. We can describe this contribution analogous to Eq. (2.16):

[TABLE]

where $\tau_{\rm ls}$ is the conformal time of the photon last scattering. Although this effect must exist in the standard cosmology, it was overlooked in most previous studies on the spectrum of primordial GWs. This is because the photon last scattering occurs after the time of matter-radiation equality and the damping effect soon becomes irrelevant as the energy density of photons becomes small compared to the total energy density of the universe. Nevertheless, we take account of this contribution in the analysis in this paper for the sake of completeness. We actually see that the inclusion of the contribution of photons to the anisotropic stress leads to an additional damping of the amplitude of GWs by $\lesssim 14\,\%$ at $f\sim 10^{-17}\,\mathrm{Hz}$ as shown in Sec. 4.

Another important effect is the change of the equation of state of the universe [18, 19]. To see this effect, let us consider the evolution of the energy density of GWs relative to that of the total energy density of the universe. From Eq. (2.8), we have

[TABLE]

where $\tilde{\rho}_{\rm gw}\equiv d\rho_{\rm gw}/d\ln k$. The evolution of the total energy density of the universe is determined by the equation of state [see Eq. (B.7)],

[TABLE]

On the other hand, if there is no energy flow due to the anisotropic stress, the energy density of GWs after the horizon crossing is simply diluted as

[TABLE]

Note that the above equation does not exactly hold around the time of the horizon crossing and that there would be some modifications according to the behavior of the equation of state parameter $w$ at that time. Although we stick to the approximation (2.20) in this section, we quantify the corresponding modifications via numerical methods in Sec. 4. Combining Eqs. (2.18)-(2.20), we obtain

[TABLE]

From Eqs. (2.9) and (2.11), we see that the scale dependence of $\Omega_{\rm gw}$ at the time of the horizon crossing is determined by the primordial tensor power spectrum: $\Omega_{\rm gw}(k,\tau_{\rm hc})\sim\mathcal{P}_{T}(k)$. Therefore, some feature of the present spectrum of GWs that is not caused by the primordial power spectrum can be characterized by the change of the equation of state parameter $w$.

Any non-trivial evolution of the equation of state after inflation can lead to the modification of the spectrum of GWs. For example, the inflaton field may undergo an oscillation around the minimum of its potential before it decays into radiations to reheat the universe after inflation. Such a dynamics of the inflaton field at the reheating stage causes some modification of the equation of state and results in a damping feature in the spectrum of GWs [19, 37, 38, 39, 40, 41].

The equation of state of the universe can vary even in the radiation dominated era, which leads to a deviation from the flat spectrum of GWs. We have already seen this fact in Eq. (2.13), which explicitly shows that there exists non-trivial frequency dependence according to the values of $g_{*\rho,\mathrm{hc}}$ and $g_{*s,\mathrm{hc}}$. Indeed, it is straightforward to show that the frequency dependence caused by $g_{*\rho,\mathrm{hc}}$ and $g_{*s,\mathrm{hc}}$ is equivalent to what shown in Eq. (2.21). For the modes $k_{1}$ and $k_{2}$ that reenter the horizon at $\tau_{\rm hc,1}$ and $\tau_{\rm hc,2}$, respectively, we obtain (see Appendix B for a more explicit proof)

[TABLE]

If the primordial tensor power spectrum is completely flat, we have $\Omega_{\rm gw}(\tau_{\rm hc,1},k_{1})=\Omega_{\rm gw}(\tau_{\rm hc,2},k_{2})$, and the variation of the amplitude of GWs at the interval between $k_{1}$ and $k_{2}$ is mostly determined by the integration of the equation of state parameter, or equivalently, the change of the effective degrees of freedom.

The change of the effective degrees of freedom can be understood by noting that the contribution of a particle with mass $m$ to radiations is exponentially suppressed when the temperature $T$ becomes smaller than $m$. For this reason, we expect that the values of $g_{*\rho,\mathrm{hc}}$ and $g_{*s,\mathrm{hc}}$ change at $T_{\rm hc}\sim m$, and that there exits a damping of the amplitude of GWs at the corresponding frequency. If the values of $g_{*\rho}$ and $g_{*s}$ are almost the same at high temperatures, which is the case in the SM, the damping factor can be approximated as $(g_{*\rho,\mathrm{hc}})(g_{*s,\mathrm{hc}})^{-4/3}\approx(g_{*\rho,\mathrm{hc}})^{-1/3}$.

In addition to the effect mentioned above, the change of the equation of state and effective degrees of freedom can also be caused due to the interaction properties of elementary particles composing radiations. If the interaction among particles is weak enough such that we can treat them as ideal gases and we ignore the change of the effective degrees of freedom caused by the finite particle masses, the equation of state is exactly given by $w=1/3$. However, this description breaks down if we take account of the interaction among elementary particles. In Sec. 3, we explicitly show that the inclusion of particle interactions causes shifts in the values of $g_{*\rho}$ and $g_{*s}$ as well as that of the equation of state parameter $w$ even if the number of particle species of the system remains unchanged. This effect leads to a further non-trivial change in the spectrum of GWs, and we discuss it more quantitatively in Sec. 4.

Once we take the interaction properties of elementary particles seriously, we may also speculate that the perfect fluid description of background radiations does not hold exactly. When we deal with such imperfect fluids, we have to include correction terms in the energy-momentum tensor characterized by three quantities, heat conduction, shear viscosity, and bulk viscosity. Although the heat conduction and bulk viscosity do not affect the evolution of GWs as far as we consider the linear order in perturbation theory, the existence of shear viscosity (parameterized by the shear viscosity coefficient $\eta$) modifies the propagation equation of GWs. In the flat FRW universe, we have [42]

[TABLE]

This fact implies that the amplitude of GWs decays in a viscous medium [43] at a rate $\Gamma_{g}=16\pi G\eta$. However, this damping is inefficient in the expanding universe for the following reasons. On dimensional grounds, the shear viscosity should behave as $\eta\sim\rho/\Gamma$ with $\Gamma$ being the interaction rate of constituent particles. Using this estimate in the damping rate and assuming the radiation dominated background ($H^{2}\sim G\rho$), we obtain $\Gamma_{g}/H\sim G\rho/\Gamma H\sim H/\Gamma$. Therefore, the damping rate due to the shear viscosity is negligible compared with the expansion rate as long as interactions of constituent particle proceed much faster than the cosmic expansion.333It was argued that macroscopic hydrodynamic fluctuations can produce GWs, whose production rate is proportional to the shear viscosity coefficient $\eta$ [44]. Although such a GW background may exist even in the SM, here we do not consider this effect, since the spectrum shows a peak at a GHz frequency range and it becomes negligibly small at lower frequencies considered in this paper.

In contrast to the shear viscosity, the bulk viscosity does not modify the propagation of GWs, but it modifies the equation for the background expansion [i.e. the evolution of the scale factor $a(t)$]. In the FRW universe, the effect of the bulk viscosity (parameterized by the bulk viscosity coefficient $\zeta$) is described by making the following replacement for the pressure [42],

[TABLE]

Hence, the whole effect results in the shift of the equation of state parameter [30],

[TABLE]

However, by using a rough estimate $\zeta\sim\rho/\Gamma$ we again see that the above correction is negligible as long as $\Gamma\gg H$. Furthermore, it is generically expected that the bulk viscosity coefficient $\zeta$ is much smaller than the shear viscosity coefficient $\eta\sim\rho/\Gamma$ for high temperature radiations [45], and the net effect can be further suppressed. For these reasons, we can safely neglect the correction due to the bulk viscosity.

Based on the above discussions, we treat the primordial plasma as a perfect fluid, but we still allow a deviation from the ideal gas description. In the next section, we explore the latter possibility in more detail and quantify its effect by taking account of the particle interactions in the SM.

3 Equation of state in the Standard Model

In order to estimate the equation of state of the early universe filled by the gases of interacting SM particles, we need to follow several steps according to the thermal history predicted by the SM. In particular, there are a couple of drastic changes according to two kinds of phase transitions: the QCD phase transition and electroweak phase transition. We note that the nature of these phase transitions might not be so drastic as expected: The result of lattice QCD simulation indicates that the nature of the QCD phase transition is not a sharp transition but a smooth crossover [46], and the observed value of the Higgs boson mass $m_{h}\approx 125\,\mathrm{GeV}$ [47] implies that the nature of the electroweak symmetry breaking is also a smooth crossover rather than a strong first order phase transition [48]. Hereinafter, we basically rely on these facts and follow an approach based on the equilibrium thermodynamics, referring to two kinds of phase transitions as the QCD crossover and electroweak crossover.

In the literature, there has been a lot of progress on the calculation of thermodynamic properties of the primordial plasma in the SM. At temperatures far above the QCD crossover, we can adopt the perturbative method to estimate the pressure of the quark-gluon plasma. This was investigated in Refs. [49, 50, 51, 52, 53, 54, 55, 56, 57], and the pressure of QCD with massless quarks has been estimated up to $\mathcal{O}(g_{s}^{6}\ln(1/g_{s}))$, where $g_{s}$ is the gauge coupling of the strong interactions. Furthermore, at temperatures above the electroweak crossover, it is also possible to estimate the pressure of the SM plasma via the perturbative method, and the detailed calculations were performed in Refs. [58, 59, 60]. These perturbative approaches break down when the temperature of the universe becomes comparable to the critical temperatures of the QCD crossover and electroweak crossover, and we need to adopt some non-perturbative methods to investigate the thermodynamics of interacting plasma. Recently, there has been some progress in estimation of the equation of state in such non-perturbative phases. These include the investigation of the equation of state of $2+1$ [61, 62, 63] and $2+1+1$ [64] flavor lattice QCD in the continuum limit, and that of the electroweak crossover with the physical Higgs mass of $125\,\mathrm{GeV}$ on the lattice [65, 66].

In this section, we collect the results obtained in the literature and combine them to estimate the effective degrees of freedom for the energy density $g_{*\rho}(T)$ and entropy density $g_{*s}(T)$ of the SM at arbitrary temperatures. Note that the energy density $\rho(T)$ and entropy density $s(T)$ can be derived from the pressure $p(T)$ by using thermodynamic equations [see Eqs. (B.3) and (B.4)]. Therefore, practically we need to know $p(T)$ and its derivative with respect to $T$ by taking account of the interaction among particles in the SM. To this end, it becomes convenient to introduce the following quantity:

[TABLE]

where the second equality follows from the thermodynamic equation [Eq. (B.4)]. The quantity $\Delta(T)$ is sometimes called the trace anomaly, since it is related to the trace of the stress-energy tensor $T^{\mu}_{\mu}=\rho-3p$, which vanishes at classical level for scale-invariant systems. Once we obtain $p(T)$ and $\Delta(T)$, the energy density and entropy density can be estimated as

[TABLE]

and they are translated into $g_{*\rho}(T)$ and $g_{*s}(T)$ via Eq. (1.2).

As mentioned above, we could obtain all relevant quantities from a single function $p(T)$ if we knew the form of $p(T)$ exactly as a differentiable function at all temperatures. However, in some cases it becomes difficult to obtain the value of $p(T)$ unambiguously, as we discuss in the following subsections. For this reason, we take $\Delta(T)$ as a fundamental quantity and estimate $p(T)$ by integrating $\Delta(T)$ with appropriate boundary conditions.

In the literature, there are pioneering works to apply the equation of state in the SM to the calculation of the Weakly Interacting Massive Particle (WIMP) dark matter abundance [67, 68] and that of the axion dark matter abundance [69, 64]. We basically follow a similar approach, but take different steps to estimate the effective degrees of freedom and their uncertainty at high temperatures up to those far above the electroweak crossover, as our main focus in this paper is to apply it to the calculation of the spectrum of primordial GWs.444See Refs. [18, 70] for some earlier attempts to apply the equation of state for strongly interacting particles to the calculation of the spectrum of primordial GWs.

3.1 Photons and leptons

Let us start by considering the contribution of photons and leptons. First of all, we need to take account of the effects of neutrino decoupling and $e^{+}e^{-}$ annihilation, which are relevant at temperatures below a few $\mathrm{MeV}$. Neutrinos decouple from the thermal bath when their weak interaction rate becomes smaller than the expansion rate of the universe at around $T\sim 2\,\mathrm{MeV}$, and after that $e^{+}e^{-}$ annihilation takes place when the temperature of the universe becomes comparable to the electron mass $m_{e}\simeq 0.511\,\mathrm{MeV}$.

The energy density of radiations for $T\ll 1\,\mathrm{MeV}$ is conventionally parameterized by the effective neutrino degrees of freedom $N_{\rm eff}$,

[TABLE]

where $\rho_{\gamma}=(\pi^{2}/15)T^{4}$ is the energy density of photons. If we assume that the neutrino decoupling has happened instantaneously and that they do not couple with the electromagnetic plasma when $e^{+}e^{-}$ are annihilated, it is straightforward to obtain $N_{\rm eff}=3$. However, neutrinos are still slightly interacting at that epoch, and the value of $N_{\rm eff}$ deviates from $3$. Detailed calculations [71, 72, 73] showed that the residual interactions lead to a distortion in the neutrino spectra, and that neutrinos (photons) have higher (lower) energy density with respect to those obtained in the instantaneous decoupling approximation. To evaluate the value of $N_{\rm eff}$ in the non-instantaneous decoupling scenario, it is convenient to recast Eq. (3.3) in the following form [71],

[TABLE]

where $\rho_{\nu}^{0}$ denotes the energy density of a single specie of neutrinos with an equilibrium distribution, $\delta\rho_{\nu_{i}}$ represent deviations from $\rho_{\nu}^{0}$ for three neutrino species $i=e,\mu,\tau$, $z\equiv Ta$ is the dimensionless photon temperature with a scale factor $a$ being normalized such that $a\to 1/T$ at large temperatures, and $z_{0}$ and $z_{\rm fin}$ represent asymptotic values for the case of instantaneous decoupling and that of full numerical calculation, respectively. Recent results of numerical calculations including the effect of neutrino oscillations [73] showed that $\delta\rho_{\nu_{i}}/\rho_{\nu}^{0}$ increase by $0.5\textendash 0.7\%$ and that the photon temperature becomes $z_{\rm fin}=1.39779$, which is smaller than $z_{0}=(11/4)^{1/3}\simeq 1.40102$ obtained in the instantaneous decoupling approximation. These values result in the following estimate for the effective neutrino degrees of freedom:

[TABLE]

This corresponds to the asymptotic value of the effective degrees of freedom for the energy density $g_{*\rho,\mathrm{fin}}\simeq 3.383$.

In principle, it should also be possible to estimate an asymptotic value of the effective degrees of freedom for the entropy density, $g_{*s,\mathrm{fin}}$. To the best of our knowledge, however, there is no clear consensus as to the value of $g_{*s,\mathrm{fin}}$ compared to the value of $g_{*\rho,\mathrm{fin}}$ (or $N_{\rm eff}$) in the literature. Here we adopt the following naive argument based on the conservation of entropy. In the numerical calculations in Refs. [71, 72, 73], the initial time was set to be $a_{\rm in}^{-1}=10\,\mathrm{MeV}$, at which the dimensionless photon temperature was given by $z_{\rm in}=1.00003$ [74]. At that time, neutrinos are regarded to be in equilibrium with the electromagnetic plasma, and hence we can at least evaluate the entropy density at the initial time straightforwardly. Then, with the definition of the dimensionless variable $z=Ta$, the conservation of entropy of the system implies

[TABLE]

where $g_{*s,\mathrm{EM}}(T_{\rm in})$ parameterizes the contributions of the electromagnetic plasma to the entropy density at the temperature $T_{\rm in}=z_{\rm in}\cdot 10\,\mathrm{MeV}$. This quantity can be estimated as

[TABLE]

where

[TABLE]

with $i=e$ is the contribution of electrons/positrons to the entropy density, and $p^{0}_{F,i}(T)$ and $\Delta_{F,i}^{0}(T)$ represent the contributions of fermion species $i$ with mass $m_{i}$ and internal degrees of freedom $g_{i}$ to the pressure and trace anomaly, respectively, in the ideal gas approximation:

[TABLE]

The functions $I_{F1}(y)$ and $I_{F2}(y)$ are given by

[TABLE]

Note that $I_{F1}(y)$ and $I_{F2}(y)$ exhibit asymptotic behaviors,

[TABLE]

In Eq. (3.7), we also introduced corrections from quantum electrodynamics (QED),

[TABLE]

Here we adopt leading order QED corrections suggested in Refs. [75, 76, 71]:

[TABLE]

where $E_{k}\equiv\sqrt{k^{2}+m_{e}^{2}}$,

[TABLE]

are finite temperature contributions to the electron mass and photon mass, respectively, and $\alpha=e^{2}/(4\pi)$ is the fine structure constant. In Eq. (3.16), we have neglected the momentum-dependent term as suggested in Ref. [71]. Substituting the values $z_{\rm in}=1.00003$ and $z_{\rm fin}=1.39779$ to Eq. (3.6), we obtain $g_{\rm*s,\mathrm{fin}}\simeq 3.931$.

Although the results of the numerical calculations mentioned above appear to be reasonable, we warn that there is a subtlety in estimating accuracy of them. Actually, we find that the value of $g_{*s,\mathrm{fin}}$ can increase by $0.1\%$ when we add the contribution of muons with a mass $m_{\mu}\simeq 105.7\,\mathrm{MeV}$ to $g_{*s,\mathrm{EM}}(T_{\rm in})$. Since the effect of muons was neglected in Refs. [74, 71, 72, 73], it is possible that all numerical results are subjected to similar uncertainty. Regarding this fact, here we adopt the above estimate with $0.1\%$ uncertainty:

[TABLE]

For the same reason, we use the following value for the effective degrees of freedom for the energy density and their uncertainty:

[TABLE]

If we assume that neutrinos have decoupled completely from the thermal bath at the epoch of $e^{+}e^{-}$ annihilation and that the electromagnetic plasma behaves as an ideal gas, the temperature dependence of $\rho(T)$ during this process can be described by the following formula [see, e.g. Ref. [77]]:

[TABLE]

where

[TABLE]

and

[TABLE]

with $i=e$ is the contribution of electrons/positrons to the energy density in the ideal gas approximation. The function $\mathcal{S}(m_{e}/T)$ describes the temperature dependence of the ratio between the effective neutrino temperature $T_{\nu}$ and photon temperature,

[TABLE]

and this relation follows from entropy conservation of the electromagnetic plasma and the scaling behavior of the effective neutrino temperature $T_{\nu}\propto 1/a$. Similarly, the temperature dependence of $s(T)$ in the instantaneous decoupling scenario is given by

[TABLE]

Note that the function $\mathcal{S}(y)$ takes the following asymptotic values:

[TABLE]

Equations (3.20) and (3.24) should be modified if we take account of the relic interactions between neutrinos and electromagnetic plasma. In principle, it is possible to follow the evolution of $\rho(T)$ and $s(T)$ including all necessary physics by solving coupled Boltzmann equations numerically, which is out of the scope of this paper. Instead, here we modify Eqs. (3.20) and (3.24) such that $g_{*\rho}(T)$ and $g_{*s}(T)$ exhibit appropriate asymptotic behavior:

[TABLE]

where the coefficients $a_{r1}$, $a_{r2}$, $a_{s1}$, and $a_{s2}$ parameterize deviations from the instantaneous decoupling results, and

[TABLE]

is the QED corrections to the energy density. The values of $a_{r1}$ and $a_{r2}$ can be fixed by requiring that the first two terms in the right-hand side of Eq. (3.26) approach asymptotically to $g_{*\rho,\mathrm{fin}}$ for $T\ll m_{e}$ and to $2+\frac{7}{8}\cdot 6$ for $T\gg m_{e}$. Adopting these conditions, we obtain

[TABLE]

where we have substituted $g_{*\rho,\mathrm{fin}}=3.383$ in the last equality. Similarly, requiring that the first two terms in the right-hand side of Eq. (3.27) approach asymptotically to $g_{*s,\mathrm{fin}}$ for $T\ll m_{e}$ and to $2+\frac{7}{8}\cdot 6$ for $T\gg m_{e}$, we obtain

[TABLE]

where we have substituted $g_{*s,\mathrm{fin}}=3.931$ in the last equality. In Eqs. (3.26) and (3.27), we also added the contributions of muons [given by Eqs. (3.22) and (3.8) with $i=\mu$] and those of hadrons ($\rho_{\rm hadrons}$ and $s_{\rm hadrons}$) in order to accomplish a smooth connection to the calculations at higher temperatures. We discuss more details about the contributions of hadrons in the next subsection.

In Fig. 1, we compare the simplistic results (3.20) and (3.24) with modified functions (3.26) and (3.27). Although the corrections look small at low temperatures, they become relevant for $T\gtrsim 10\,\mathrm{MeV}$. Actually, we see that the values of $g_{*\rho}$ and $g_{*s}$ increase by $\mathcal{O}(1)\%$ at $T=20\,\mathrm{MeV}$ due to the contributions of muons and light hadrons.

For temperatures below $10\,\mathrm{MeV}$, we evaluate the effective degrees of freedom by using Eqs. (3.26) and (3.27) with coefficients specified by Eqs. (3.29), (3.30), (3.31), and (3.32), and take their uncertainty by simply multiplying the functions (3.26) and (3.27) by a factor 1.001 or 0.999. Here we neglect the effect of the neutrino masses, though they can lead to further modifications at $T\lesssim 0.1\,\mathrm{eV}$ if some of them are close to the current cosmological bound on the sum $\sum_{j}m_{j}<0.17\,\mathrm{eV}$ reported by the Planck Collaboration [78, 79] (including data on the Baryon Acoustic Oscillations). For $T>10\,\mathrm{MeV}$, on the other hand, neutrinos are kept in equilibrium with the thermal bath via weak interactions, and hence in such higher temperatures we evaluate the energy density and entropy density in terms of thermodynamic quantities $p(T)$ and $\Delta(T)$ [i.e. via Eq. (3.2)] instead of using Eqs. (3.26) and (3.27). In particular, we estimate the pressure $p(T)$ by integrating the trace anomaly $\Delta(T)$ as we mentioned in the beginning of this section. Note that the contribution of taus with a mass $m_{\tau}\simeq 1.777\,\mathrm{GeV}$ becomes relevant for $T\gtrsim\mathcal{O}(0.1)\,\mathrm{GeV}$. We include this contribution and calculate $p(T)$ and $\Delta(T)$ at higher temperatures in Sec. 3.4 after we discuss the effects of the QCD crossover and electroweak crossover in the next subsections.

3.2 Strongly interacting particles

The most drastic change happens at the epoch of the QCD crossover. Here we can consider three different phases: the hadronic phase [$T\lesssim\mathcal{O}(100)\,\mathrm{MeV}$], the non-perturbative phase [$T\sim\mathcal{O}(100)\,\mathrm{MeV}$], and the perturbative phase [$T\gg\mathcal{O}(100)\,\mathrm{MeV}$]. In this subsection, we discuss how to estimate the contribution of strongly interacting particles in each phase. Interpolation between different phases should be considered carefully, and we quantify the uncertainty arising from the interpolation procedure in Sec. 3.4.

The thermodynamic properties of hadronic plasma at temperatures below the critical temperature of the QCD crossover can be described by the hadron resonance gas model, in which the system is approximated by free hadrons and resonances [80, 81, 82]. The hadron resonance gas model has shown a good agreement with lattice QCD results [62, 83] and successfully reproduced the experimental results observed in heavy ion collisions [see, e.g. Refs. [84, 85] and references therein]. In this model, thermodynamic quantities are expressed as sums over single free particle contributions with mass $m_{i}$ and internal degrees of freedom $g_{i}$. The pressure and trace anomaly read

[TABLE]

where

[TABLE]

[TABLE]

$I_{F1}(y)$ and $I_{F2}(y)$ are defined in Eq. (3.11), and $I_{B1}(y)$ and $I_{B2}(y)$ are given by

[TABLE]

These functions exhibit the following asymptotic behaviors,

[TABLE]

In the sums in Eqs. (3.33) and (3.34), we include all known baryons and mesons whose masses are less than $2.5\,\mathrm{GeV}$ as enumerated in the Particle Data Group’s review [86]. The energy density and entropy density can be estimated as $\rho_{\rm hadrons}(T)=T^{4}[\Delta_{\rm hadrons}(T)+3p_{\rm hadrons}(T)/T^{4}]$ and $s_{\rm hadrons}(T)=T^{3}[\Delta_{\rm hadrons}(T)+4p_{\rm hadrons}(T)/T^{4}]$, respectively, and they were used in Eqs. (3.26) and (3.27).

The description based on the hadron resonance gas model breaks down at around the critical temperature of the QCD crossover. At that stage, interactions among constituent particles become very strong, and we must use the non-perturbative method, i.e. lattice QCD simulations, to estimate the thermodynamic properties of quark-gluon plasma. The equation of state in $2+1$ flavor lattice QCD with physical quark masses was calculated by the “Budapest-Wuppertal” collaboration [61, 62], and their result was confirmed by the “HotQCD” collaboration in Ref. [63]. Recently, the contribution of the charm quark was investigated in $2+1+1$ flavor lattice QCD [64], which enables us to estimate the equation of state for higher temperatures up to $1\,\mathrm{GeV}$. In this paper, we adopt the result of Ref. [64] to calculate the effective degrees of freedom in the non-perturbative phase.555 For $2+1$ flavor lattice QCD results plotted in Fig. 2, we used the tabulated data attached to Ref. [62] [dubbed as “Budapest-Wuppertal (2013)”] and table I in Ref. [63] [dubbed as “HotQCD (2014)”]. For $2+1+1$ flavor lattice QCD results, we extract the data from Fig. S9 (or Fig. S7 in the arXiv version) of Ref. [64] [dubbed as “Budapest-Wuppertal (2016)”].

For temperatures much higher than the critical temperature of the QCD crossover, we can use the perturbative method of finite temperature QCD. The pressure for QCD with massless quarks at high temperature $p_{\rm QCD}(T)$ was computed by using effective field theory approach [87, 88, 56], which provides a systematic framework for combining perturbative and non-perturbative results in a consistent way. The coefficients of the expansion in the strong gauge coupling $g_{s}$ have been computed at orders $g_{s}^{2}$ [49, 50], $g_{s}^{3}$ [51], $g_{s}^{4}\ln(1/g_{s})$ [52], $g_{s}^{4}$ [53, 54], $g_{s}^{5}$ [55, 56], and $g_{s}^{6}\ln(1/g_{s})$ [57]. Here we quote the results obtained in Ref. [57]:

[TABLE]

where $\alpha_{s}=g_{s}^{2}/4\pi$,

[TABLE]

$N_{f}$ is the number of massless quarks, and $\mu$ is the renormalization scale in the modified minimal subtraction ($\overline{\text{MS}}$) scheme. For the evaluation of the renormalization group running of gauge coupling $g_{s}$, we adopt four-loop beta function coefficients obtained in Ref. [89]. The contribution to the trace anomaly can be evaluated as

[TABLE]

The coefficient $q_{c}(N_{f})$ in Eq. (3.49), which might depend on $N_{f}$, is unknown, and it is related to the serious infrared problem of finite temperature field theory (known as “Linde’s problem”) [90, 91]. A power counting argument implies that the contributions of $(l+1)$-loop gluon diagrams to the free energy (or pressure) contain infrared divergences of the form $\sim g_{s}^{6}T^{4}\ln(T/m)$ for $l=3$ and $\sim g_{s}^{6}T^{4}(g_{s}^{2}T/m)^{l-3}$ for $l>3$, where $m$ is an infrared cutoff scale. If the system were described by QED, the cutoff scale could be identified as the screening mass $m\sim g_{s}T$, which would regularize all the infrared divergences and make the perturbative expansion in terms of $g_{s}$ totally controllable. However, this is not the case in QCD due to the self-coupling of the spacelike gluons, which give rise to the magnetic mass $m\sim g_{s}^{2}T$ for the cutoff scale. Consequently, all contributions with $l>3$ become of $\mathcal{O}(g_{s}^{6})$, independently of the order of perturbation theory. In other words, we cannot determine the contribution of $\mathcal{O}(g_{s}^{6})$ in perturbative approach, and such a contribution is parameterized by the coefficient $q_{c}(N_{f})$ in Eq. (3.49).

In Fig. 2, we show the contributions of strongly interacting particles to the pressure and trace anomaly. We see that the results based on the hadron resonance gas model [Eqs. (3.33) and (3.34)] agree with those obtained from lattice QCD simulations at $T\gtrsim 0.1\,\mathrm{GeV}$. Both $p(T)$ and $\Delta(T)$ blow up at higher temperatures in the hadron resonance gas model, indicating that the system should be described as interacting quark-gluon plasma rather than a collection of free hadrons and resonances and that we have to adopt lattice data. In Fig. 2, we plot three different sets of lattice data, in which the analysis was performed with physical quark masses and the results were shown in the continuum limit. These include the results of $2+1$ flavor lattice QCD in Refs. [62, 63] and those of $2+1+1$ flavor lattice QCD in Ref. [64]. The results of $2+1$ flavor lattice QCD are consistent with those of $2+1+1$ flavor lattice QCD at low temperatures, while they disagree at higher temperatures as the contribution of charm quarks (mass $m_{c}\simeq 1.3\,\mathrm{GeV}$) becomes important. This $2+1+1$ flavor lattice QCD results should be connected with the perturbative results [Eqs. (3.43) and (3.50)] at some higher temperature.666It should be noted that the perturbative results for $\Delta_{\rm QCD}(T)$ shown in Fig. 2 and used in the analysis in Sec. 3.4 do not represent a genuine perturbative expansion of the trace anomaly but a direct differentiation of pressure $p_{\rm QCD}(T)$ obtained in the perturbative expansion. Namely, “$\mathcal{O}(g_{s}^{n})$” result of $\Delta_{\rm QCD}(T)$ shown in the lower panel of Fig. 2 is obtained by performing the temperature derivative in Eq. (3.50), where $p_{\rm QCD}(T)$ includes terms up to $\mathcal{O}(g_{s}^{n})$ and the renormalization group running of $g_{s}$ is also treated as a function of $T$. Strictly speaking, the latter should be regarded as higher order effects, but we include them in $\Delta_{\rm QCD}(T)$ since their effects are not negligible in the temperature range considered here. However, the convergence of the perturbative expansion turns out to be very poor at $T\sim\mathcal{O}(1\textendash 10)\,\mathrm{GeV}$. In fact, the prediction of the perturbation theory changes significantly at each order of weak coupling expansion, as shown in Fig. 2. Furthermore, we must adjust the value of the unknown constant $q_{c}(N_{f})$ at order $g_{s}^{6}$. Although the choice $q_{c}(N_{f}=4)=-3000$ [64] appears to coincide with the $2+1+1$ flavor lattice QCD result as shown in Fig. 2, large uncertainty of the perturbative results around $T\sim\mathcal{O}(1)\,\mathrm{GeV}$ leaves the possibility of other choices, which we investigate more carefully when we match the perturbative results to those of the $2+1+1$ flavor lattice QCD simulations in Sec. 3.4.

Before going to the matching procedure between the perturbative and non-perturbative phases in QCD, we must consider the effect of quark mass thresholds, which should be included in order to describe the contributions of heavier quarks correctly. The perturbative results shown in Eqs. (3.43)-(3.49) are valid only for massless quarks, and we cannot use them when the temperature becomes comparable to the masses of heavier quarks, since these heavier states start to contribute to the pressure. The corrections to the QCD pressure with massive quarks were estimated up to $\mathcal{O}(g_{s}^{2})$ in Ref. [92], and it was shown that higher order corrections to the quark mass dependence are insignificant compared with those to the pressure as a whole. Based on this fact, in Ref. [64] it was argued that the quark mass threshold can be sufficiently described by using a tree-level correction factor. Here we follow the approach of Ref. [64] and describe the effect of quark masses as follows. First, we write the pressure of 4-flavor QCD including the effect of the charm quark mass as

[TABLE]

where

[TABLE]

is the pressure in the Stefan-Boltzmann limit of the $N_{f}$-flavor theory, $p_{F,c}^{0}(T)$ is the contribution of free charm quarks [see Eq. (3.9)], and $p_{\rm QCD}(T)|_{N_{f}\,=\,4}$ corresponds to Eq. (3.43) evaluated at $N_{f}=4$. Next, we describe the effect of the bottom quark mass $m_{b}\approx 4.2\,\mathrm{GeV}$ as [64]

[TABLE]

where $p_{F,b}^{0}(T)$ is the contribution of free bottom quarks evaluated by using Eq. (3.9). Finally, the contribution to the trace anomaly can be estimated as

[TABLE]

Equations (3.53) and (3.54) allow us to model the variation of the pressure and trace anomaly due to the charm and bottom quark masses, and they can be used to extrapolate the results of lattice QCD simulations to higher temperatures.

In addition to the contributions of charm and bottom quarks discussed above, we also have to include the contribution of top quarks (mass $m_{t}\gtrsim 170\,\mathrm{GeV}$). Here we do not add a top quark mass threshold effect just by extrapolating the 4-flavor result as in Eq. (3.53), since at temperature much higher than $m_{t}$ the thermodynamic quantities of QCD can be exactly described by using the result with massless quarks [Eqs. (3.43) and (3.50)]. For this reason, we model the trace anomaly such that it approaches Eq. (3.54) for $T\ll m_{t}$ and Eq. (3.50) with $N_{f}=6$ for $T\gg m_{t}$,

[TABLE]

where

[TABLE]

and $p_{F,t}^{0}(T)$ is the contribution of free top quarks evaluated by using Eq. (3.9).777If we consider a trace anomaly given by a derivative of the following combination,

$p_{\rm QCD}^{(u+d+s+c+b+t)}(T)=c_{t}(T)\left.p_{\rm QCD}(T)\right|_{N_{f}\,=\,6}+\left(1-c_{t}(T)\right)p_{\rm QCD}^{(u+d+s+c+b)}(T),$

(3.57)

we have

$\displaystyle\Delta_{\rm QCD}^{(u+d+s+c+b+t)}(T)$ $\displaystyle=c_{t}(T)\left.\Delta_{\rm QCD}(T)\right|_{N_{f}\,=\,6}+\left(1-c_{t}(T)\right)\Delta_{\rm QCD}^{(u+d+s+c+b)}(T)$

$\displaystyle\quad+T\left(\frac{dc_{t}(T)}{dT}\right)\left(\frac{\left.p_{\rm QCD}(T)\right|_{N_{f}\,=\,6}-p_{\rm QCD}^{(u+d+s+c+b)}(T)}{T^{4}}\right).$

(3.58)

For a tree-level approximation, $\left.p_{\rm QCD}(T)\right|_{N_{f}\,=\,6}-p_{\rm QCD}^{(u+d+s+c+b)}(T)\approx\frac{7\pi^{2}}{60}T^{4}$, we see that the last term in the right-hand side of the above equation represents the contribution of free top quarks, $\Delta_{F,t}^{0}(T)=T\frac{d}{dT}\{p_{F,t}^{0}(T)T^{-4}\}$. Since such a contribution is automatically included once we evaluate the Higgs potential through the electroweak crossover [see Eq. (3.63)], we dropped the term proportional to $dc_{t}(T)/dT$ in Eq. (3.58) in order to avoid double counting. The correction factor $c_{t}(T)$ approaches a limit $c_{t}(T)\to 0$ for $T\ll m_{t}$ and $c_{t}(T)\to 1$ for $T\gg m_{t}$ and smoothly interpolates between $\Delta_{\rm QCD}^{(u+d+s+c+b)}(T)$ and $\left.\Delta_{\rm QCD}(T)\right|_{N_{f}\,=\,6}$. We also note that the top quark mass $m_{t}(v_{T})$ in Eq. (3.56) is evaluated by taking account of the temperature variation of a Higgs expectation value $v_{T}$ [see Eq. (3.64)]. The evolution of the expectation value $v_{T}$ can be traced together with the evaluation of the Higgs potential, which we discuss in the next subsection.

3.3 Electroweak crossover

Across the electroweak crossover, the universe evolved from high temperature symmetric phase to low temperature broken phase, in which the Higgs field has a non-vanishing expectation value. In this subsection, we discuss the additional contributions to thermodynamic quantities arising through the electroweak crossover.

Following Ref. [92], we estimate the thermodynamic quantities in the low temperature broken phase by using the finite temperature one-loop effective potential of the Higgs field [hereinafter, referred to as “one-loop Coleman-Weinberg (CW) method”]. Namely, we evaluate the free energy density $f(\varphi,T)$ in the presence of a Higgs expectation value $\varphi$ at temperature $T$ as

[TABLE]

where $\nu^{2}$ and $\lambda$ are the quadratic Higgs coupling and quartic Higgs coupling, respectively,

[TABLE]

are one-loop contributions of Higgs bosons (subscript $h$), $W$ bosons ($W$), $Z$ bosons ($Z$), and top quarks ($t$), and $I_{F1}(y)$ and $I_{B1}(y)$ are integrals defined in Eq. (3.11) and Eq. (3.41), respectively. The $\varphi$-dependent masses in Eqs. (3.60)-(3.63) are given by

[TABLE]

where $g$ is the SU(2)L gauge coupling, $g^{\prime}$ is the U(1)Y gauge coupling, and $y_{t}$ is the top Yukawa coupling. Let $v_{T}$ and $v_{0}$ be values of $\varphi$ that minimize $f(\varphi,T)$ and $f(\varphi,0)$, respectively. The pressure in the low temperature broken phase can be estimated as

[TABLE]

The contribution to the trace anomaly is given by

[TABLE]

We evaluate $\Delta_{\rm electroweak,low}(T)$ by minimizing $f(\varphi,T)$ numerically for a given value of $T$. For coupling parameters $\nu(\mu)$, $\lambda(\mu)$, $g^{2}(\mu)$, $g^{\prime 2}(\mu)$, $y_{t}^{2}(\mu)$ defined in the $\overline{\text{MS}}$ scheme, we adopt three-loop order running evaluated in Ref. [93]. The renormalization scale $\mu$ is simply taken to be a fixed value $\mu=m_{Z}$, where $m_{Z}\simeq 91.19\,\mathrm{GeV}$ is the $Z$ boson mass, and this scale is varied by a factor $0.5\dots 2$ in order to see the dependence on the choice of the renormalization scale, which can be regarded as theoretical uncertainty.

The pressure and trace anomaly in the high temperature crossover/symmetric phase were also estimated in Refs. [58, 59, 60] based on the dimensional reduction to the effective 3d theory [94], which is essentially the same with the effective field theory approach used to estimate the pressure of QCD at high temperature [87, 88, 56]. The trace anomaly in this regime can be split into three terms [60]:

[TABLE]

where888 Equation (3.68) is slightly different from the expression of $\Delta_{1}(T)$ given in Ref. [60]. Here, we dropped the terms proportional to $g_{s}^{4}$ and $g_{s}^{5}$ in $\Delta_{1}(T)$, since they are included in the pure QCD contribution $\Delta_{\rm QCD}(T)$, which we treat separately.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$n_{G}=3$ is the number of generations, $\gamma_{E}\simeq 0.577216$ is the Euler-Mascheroni constant, $\zeta(s)$ denotes the Riemann zeta function, $g_{3}^{2}=g^{2}T$ is the coupling parameter in the 3d theory,

[TABLE]

are the Debye mass parameters, and $p_{0R}$ is a renormalized expression of the zero-temperature pressure [associated with $f(v_{0},0)$ in Eq. (3.65)]. The three terms in the right-hand side of Eq. (3.67) have different origins and hence represent different effects. According to the terminology of Ref. [60], $\Delta_{1}(T)$ represents breaking of scale invariance by quantum corrections, $\Delta_{2}(T)$ represents the temperature evolution of the Higgs condensate, and $\Delta_{3}(T)$ represents the vacuum term. Mathematically, this distinction arises due to the fact that $T$ can appear in $p(T)/T^{4}$ as $\mu/T$, $\nu^{2}/T^{2}$, or $p_{0R}/T^{4}$ because of dimensional reasons. $\Delta_{1}(T)$ and $\Delta_{2}(T)$ correspond to the derivative of $p(T)/T^{4}$ with respect to $\mu/T$ and $\nu^{2}/T^{2}$, respectively. In Eq. (3.69), we also split $\Delta_{2}(T)$ into three terms for practical reasons discussed below. These three terms correspond to three different momentum scales appearing in the effective 3d theory [56, 94, 58, 59, 60]. $\Delta_{2E}(T)$ contains contributions from the “superheavy” scale $k\sim T$, $\Delta_{2M}(T)$ those from the “heavy” scale $k\sim gT$,999Hereafter we use an abbreviation $g^{2}$ to represent the expansion parameter of the perturbative analysis in the symmetric phase, i.e. $g^{2}\in\{g^{2},g^{\prime 2},g_{s}^{2},y_{t}^{2},\lambda,\nu^{2}/T^{2}\}$. and $\Delta_{2G}(T)$ those from “light” scale $k\sim g^{2}T$. The vacuum contribution $\Delta_{3}(T)$