Freezing-In Gravitational Waves

Jacopo Ghiglieri, Jan Schütte-Engel, Enrico Speranza

► To cite this version:

Jacopo Ghiglieri, Jan Schütte-Engel, Enrico Speranza. Freezing-In Gravitational Waves. Physical Review D, 2024, 109 (2), pp.023538. 10.1103/PhysRevD.109.023538. hal-03895048

HAL Id: hal-03895048

https://hal.science/hal-03895048v1

Submitted on 26 May 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Freezing-in gravitational waves

Jacopo Ghiglieri<sup>1</sup>, Jan Schütte-Engel<sup>2,3</sup> and Enrico Speranza<sup>2,3</sup>

<sup>1</sup>SUBATECH, Université de Nantes, IMT Atlantique, IN2P3/CNRS

4 rue Alfred Kastler, La Chantrerie BP 20722, 44307 Nantes, France

<sup>2</sup>Department of Physics, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA

<sup>3</sup>Illinois Center for Advanced Studies of the Universe,
University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA

<figure><img src="image_2.png" /></figure>

(Received 23 December 2022; accepted 9 January 2024; published 31 January 2024)

The thermal plasma in the early Universe produced a stochastic gravitational wave (GW) background, which peaks today in the microwave regime and was dubbed the cosmic gravitational microwave background (CGMB). In previous works, only single graviton production processes that contribute to the CGMB have been considered. Here, we also investigate graviton pair production processes and show that these can lead to a significant contribution if the ratio between the maximum temperature and the Planck mass, $T_{max}/m_p$, divided by the internal coupling in the heat bath is large enough. As the dark matter freeze-in production mechanism is conceptually very similar to the GW production mechanism from the primordial thermal plasma, we refer to the latter as “GW freeze-in production.” We show that quantum gravity effects appear in single graviton production and are smaller by a factor $(T_{max}/m_p)^2$ than the leading order contribution. In our work, we explicitly compute the CGMB spectrum within a scalar model with quartic interaction.

DOI: 10.1103/PhysRevD.109.023538

# I. INTRODUCTION

The first detection of gravitational waves (GWs) from back hole and neutron star mergers [1,2] opened up a new window to explore our Universe. While the GWs that have been detected so far were emitted in the late-time Universe, GWs can also be produced in the early Universe. These GWs are stochastic in nature, and their detection would yield unprecedented information about early Universe cosmology as well as high energy particle physics. To give a few examples, GWs in the early Universe can be produced from inflation [3–6], preheating [7,8], inflaton annihilation into gravitons [9–11], first-order phase transitions [12,13], cosmic defects such as cosmic strings [14,15], noisy turbulent motion [16–20] and equilibrated gravitons [21,22]. For a review on early Universe GW sources, see Ref. [23]. The full GW spectrum for a specific particle physics model that can describe the entire cosmological history was worked out in Ref. [24].

In this paper, we consider GWs that were produced from the thermal plasma [25–28] in the early Universe. Every

plasma, even in thermal equilibrium, produces GWs due to microscopic particle collisions and macroscopic hydrodynamic fluctuations, cf. Ref. [26]. In the former case, the GW momenta are on the order of the temperature, $k \sim T$, and in the latter case, they are much smaller, $k \ll T$. Here, we focus on GWs produced by microscopic particle collisions since it enables us to probe elementary particle physics theories at high energies. Furthermore, the GW contribution from microscopic particle collisions to the final spectrum is larger compared to the contribution from hydrodynamic fluctuations [26].

Our main assumption is that after the hot big bang, a thermal plasma of particles in thermal equilibrium at a maximum temperature $T_{max}$ was present. In addition, we assume that at this time no GWs are present. In an expanding Universe, GWs from the thermal plasma are continuously produced as the temperature decreases. The spectrum of the produced GWs peaks at a frequency on the order of the temperature at the time of production. If the redshift-temperature relation is linear, the GW spectra that are produced at different temperatures add up such that the observed GW spectrum today is enhanced. The spectrum of the produced GWs today peaks in the microwave regime and is hence dubbed the cosmic gravitational microwave background (CGMB).

In principle, the maximum temperature, $T_{max}$, of the thermal plasma can be as high as the Planck mass

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP<sup>3</sup>.
$m_p \approx 1.2 \times 10^{19}$ GeV [29]. However, in slow roll inflationary cosmology, it cannot be much higher than $10^{-3} m_p$. This bound follows by first inferring the energy scale of inflation from the amplitude of scalar perturbations and the tensor-to-scalar ratio and then assuming an instantaneous and a maximally efficient reheating [30] to a radiation dominated Universe, cf. Ref. [28] and Ref. [31] for a review. Note that $T_{max} > 10^{-3} m_p$ can be achieved in noninflationary scenarios. One particular example are bouncing cosmology scenarios which can lead to $T_{max}$ which goes up to the Planck scale: $T_{max} < m_p$; cf. Refs. [32–34]. The maximum temperature of the thermal plasma is also bounded from below. The most conservative estimates set a lower limit around a few MeV [35–39], shortly before big bang nucleosynthesis took place. However, most scenarios require temperatures reaching well above the electroweak scale such that, e.g., sphalerons can be active in leptogenesis scenarios [40].

In previous works, the CGMB spectrum has been calculated within the Standard Model (SM) [26,27] and for beyond Standard Model (BSM) theories [28,41–43]. Those works considered only single graviton production processes. In this case, the resulting GW energy density per logarithmic momentum interval, $\Omega_{gw}$, is proportional to $g^2 \times (T_{max}/m_p)$, where g is the internal coupling in the thermal bath. Here, we extend previous works by also including GW production processes with two gravitons in the final state. These give a contribution to $\Omega_{gw}$ that is proportional to $(T_{max}/m_p)^3$. We refer to this GW production channel as graviton pair production. Depending on the values of $T_{max}/m_p$ and g, the graviton pair production channel can be the dominating contribution to the CGMB spectrum. In analogy to dark matter production from the thermal plasma, we dub the GW production from the thermal plasma GW freeze-in production.

We also identify at which order quantum gravity and backreaction effects would appear in the CGMB spectrum. Observing these effects in the CGMB spectrum would therefore probe the quantization of gravity and reveal fundamental information about particle physics if the GW production occurs at high energy scales that cannot be probed with particle colliders on Earth.

Throughout this paper, we work with a complex scalar field with quartic coupling $\lambda$ that is the internal coupling in the thermal bath; i.e., for our model, the previously mentioned generic coupling g is the quartic coupling $\lambda$. In previous works [26–28], such a coupling has not been considered even though the SM has such a coupling in the Higgs sector. That is because Refs. [26,27] worked under the assumption that in the SM, the three gauge couplings

and the top Yukawa are of order of the square root of the Higgs self-coupling. Note that in BSM theories, this is not necessarily the case.

This paper is organized as follows: In Sec. II, we introduce our model, which is a complex scalar field coupled to gravity. This is then followed by Sec. III, where we introduce the full evolution equations for the two distribution functions $f_\phi$ and $f_h$, which describe the scalars and gravitons, respectively. Furthermore, we perturbatively expand the distribution functions around their initial states. This enables us to find a solution of the coupled nonlinear integral-differential equations for the distribution functions. In Sec. IV, we calculate the Matrix elements squared for the graviton production processes. We then compute the GW spectrum in Sec. V. Finally, conclusions are given in Sec. VI. Throughout this paper, we use natural units with $\hbar = c = k_B = 1$, where $k_B$ is the Boltzmann constant.

## II. SCALAR MODEL

The action for a complex scalar field on curved space-time is

$$S_{\phi} = \int d^4 x \mathcal{L}_{\phi} = \int d^4 x \sqrt{-g} (-g^{\mu\nu} (\nabla_{\mu} \phi)^{\dagger} \nabla_{\nu} \phi - U), \quad (1)$$

where $\nabla_\mu$ is a covariant derivative, $g_{\mu\nu}$ the metric tensor, $g = \text{Det}[g_{\mu\nu}]$ and U is the potential. The flat space-time metric is defined as $\eta_{\mu\nu} \equiv \text{diag}(-1, 1, 1, 1)$. Note that for a scalar field, the covariant derivative reduces to a partial derivative: $\nabla_\mu \phi = \partial_\mu \phi$. The considered complex scalar is not charged under a local transformation, and we consider a quartic potential $U = -\frac{\lambda}{4}|\phi|^4$. We assume that the scalar field is massless, which is justified if the considered temperature in the thermal plasma is larger than the mass of the scalar field.

On top of the action in Eq. (1) we need the Einstein-Hilbert action

$$S_{\text{EH}} = \int d^4 x \mathcal{L}_{\text{EH}} = \int d^4 x \frac{1}{16\pi G} \sqrt{-g} R, \quad (2)$$

where R is the Ricci scalar, and $G \equiv 1/m_p^2$ is the gravitational constant. In the following, we expand the metric around flat space-time: $g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu}$, with $h_{\mu\nu} \ll 1$. A detailed expansion has been worked out before in Ref. [44] and yields
<figure><img src="image_3.png" /><figcaption>FIG. 1. Lowest order vertices that arise from the expansion of the scalar field and Einstein-Hilbert Lagrangian. Scalars are represented by dashed lines and gravitons by double lines.</figcaption></figure>

$$ \mathcal{L} = \mathcal{L}_{\text{EH}} + \mathcal{L}_{\phi} = -\eta^{\mu\nu}(\partial_{\mu}\phi)^{\dagger}\partial_{\nu}\phi - U + \frac{1}{2}\partial_{\mu}h^{\sigma\nu}\partial^{\sigma}h_{\sigma\nu} + \kappa h^{\mu\nu}(\partial_{\mu}\phi)^{\dagger}\partial_{\nu}\phi \\ + \kappa\left(-\frac{1}{2}h^{\alpha}_{\beta}\partial_{\alpha}h^{\mu}_{\nu}\partial^{\beta}h^{\nu}_{\mu} - h^{\alpha}_{\beta}\partial_{\mu}h^{\nu}_{\alpha}\partial^{\mu}h^{\beta}_{\nu} + h^{\beta}_{\mu}\partial_{\nu}h^{\alpha}_{\beta}\partial^{\mu}h^{\nu}_{\alpha}\right) \\ + \kappa^2\left[-(h^{\mu\lambda}h^{\nu}_{\lambda} + \frac{1}{4}\eta^{\mu\nu}h^{\alpha}_{\rho}h^{\rho}_{\alpha})(\partial_{\mu}\phi)^{\dagger}\partial_{\nu}\phi + \frac{1}{4}h^{\alpha}_{\rho}h^{\rho}_{\alpha}U\right] + \mathcal{O}(h^3). \quad (3) $$

Note that the h fields in Eq. (3) have been rescaled with a factor $\kappa = \sqrt{32\pi G}$ and have now mass dimension one. Furthermore, we have adopted the so-called transverse-traceless (TT) gauge, which includes the De Donder gauge: $\partial_{\alpha}h^{\alpha}_{\mu} = \frac{1}{2}\partial_{\mu}h$ together with the requirement that the trace $h = h^{\mu}_{\mu}$ is zero. In the first and second line of Eq. (3), we wrote down the zeroth and first order terms coming from $\mathcal{L}_{\text{EH}}$ and $\mathcal{L}_{\phi}$. In the second line, we only write down the second order term coming from $\mathcal{L}_{\phi}$ since the second order term from $\mathcal{L}_{\text{EH}}$ will not be needed for our calculations. The lowest order Feynman vertices for our theory are shown in Fig. 1.

### III. EVOLUTION EQUATIONS FOR THE DISTRIBUTION FUNCTIONS

We describe the thermal plasma of $\phi$ particles and the produced gravitons with two distribution functions defined as

$$ f_{\phi}(t, k) \equiv \frac{N_{\phi}^{k}}{V d^{3} k / (2\pi)^{3}}, \quad f_{h}(t, k) \equiv \frac{N_{h}^{k}}{V d^{3} k / (2\pi)^{3}}. \quad (4) $$

$V$ is the considered volume, and $N_{\phi}^k$ and $N_h^k$ are the numbers of $\phi$-states and gravitons with momentum $k = |k|$ in the interval $d^3k$. Since we shall expand $f_{\phi}$ around an isotropic equilibrium state, $f_h$ is understood to be the polarization-averaged distribution function. Also we do not introduce a distribution function for $\phi^\dagger$ since our model and initial conditions are CP symmetric, and therefore, it would always be equal to $f_{\phi}$.

In the regime where the momentum $k$ is on the order of the hard scale that in equilibrium corresponds to the temperature, i.e., $k \sim T$, kinetic theory is expected to be a good approximation for our system. The evolution equations for the $\phi$ and graviton distribution functions can thus be written in the following Boltzmann-like form:

$$ \dot{f}_{\phi}(t, k) = G_{\phi}(t, k) - L_{\phi}(t, k), \quad (5) $$

$$ \dot{f}_{h}(t, k) = G_{h}(t, k) - L_{h}(t, k), \quad (6) $$

where the $G$ and $L$ terms describe the gain and loss terms of particle states. A generic expression for the graviton production term $G_h$ is given by

$$ G_h(t, k) = \frac{1}{4k} \sum_{\substack{\text{all processes } r \\ \text{with at least one} \\ \text{final state graviton}}} S_r \int d\Omega_r |\mathcal{M}_r|^2 \times f_{\phi}(p'_1) \cdots f_{\phi}(p'_m) f_h(k'_1) \cdots f_h(k'_n) \\ \times (1 + f_{\phi}(p_1)) \cdots (1 + f_{\phi}(p_i))(1 + f_h(k)) \cdots (1 + f_h(k_j)), \quad (7) $$

where the index $r$ labels all possible processes. We call the momenta of the incoming $\phi$ and graviton states $p'_1, \dots, p'_m$ and $k'_1, \dots, k'_n$, respectively. The momenta of the outgoing $\phi$'s and gravitons are $p_1, \dots, p_i$ and $k_1 = k, k_2, \dots, k_j$. In our notation, the $m$ incoming and $i$ outgoing $\phi$ states can be

$\phi$ or $\phi^\dagger$ states. In Eq. (7), the symmetry factor $S_r$ has to be included if two or more indistinguishable particles appear in the initial or final state. We will make the symmetry factor explicit in section Sec. IV where we calculate the graviton rate. The sum in Eq. (7) runs over combinations of
all processes with at least one graviton with momentum $k$ in the final state. The prefactor $1/(4k)$ is a combination of $1/(2k)$ from the phase space measure and $1/2$ from the graviton polarization degeneracy. We need the factor 2 from the polarization degeneracy since the matrix element squared is summed over polarizations and the distribution function is defined to be averaged over both polarizations.

The loss term is analogous to the gain term in Eq. (7) with the difference that one sums over all processes with at least one graviton in the initial state. The Boltzmann-like Eqs. (5), (6), and (7) come with important caveats on their validity beyond leading order, which we shall discuss later. The integral that appears in $G_h$ is the phase space integral that has to be performed over all momenta, except $k$:

$$ \int d\Omega_r = \int \frac{d^3 p'_1}{(2\pi)^3 2p'_1} \cdots \int \frac{d^3 p'_m}{(2\pi)^3 2p'_m} \int \frac{d^3 k'_1}{(2\pi)^3 2k'_1} \cdots \int \frac{d^3 k'_n}{(2\pi)^3 2k'_n} \int \frac{d^3 p_1}{(2\pi)^3 2p_1} \cdots \int \frac{d^3 p_i}{(2\pi)^3 2p_i} \\ \times \int \frac{d^3 k_2}{(2\pi)^3 2k_2} \cdots \int \frac{d^3 k_j}{(2\pi)^3 2k_j} \times (2\pi)^4 \delta^{(4)}(P'_1 \cdots + P'_m + K'_1 \cdots + K'_n - P_1 \cdots - P_i - K \cdots - K_j), \quad (8) $$

where we use capital letters to denote four-vectors. For further use, we introduce the shorthand notation:

$$ \int_r := \frac{S_r}{4k} \int d\Omega_r |M_r|^2. \quad (9) $$

We have written Eqs. (5), (6), (7), and (8) in a rather generic form which includes all possible processes. In our specific model of a massless complex scalar field, $1 \leftrightarrow 2$ processes are only allowed in the collinear limit, i.e., when the three-momenta of all three particles are exactly parallel. However, in this case, the thermal and vacuum masses of the scalars have to be taken into account, which leads to the fact that $1 \leftrightarrow 2$ graviton production processes are not even allowed in the collinear limit. The first kinematically allowed processes for graviton production are $2 \to 2$ and $2 \leftrightarrow 3$ processes, which have a lowest order matrix element squared of $O(\kappa^4)$ and $O(\lambda^2\kappa^2)$, respectively.

We are interested in tracking the evolution of $f_h$ under the assumption that it starts from an initially vanishing value in a bath of equilibrated scalars.<fnref n="1" /> In our stated freeze-in scenario, we further assume that throughout the entire evolution $f_h \ll 1$ and $|n_B - f_\phi| \ll 1$, where $n_B$ is the Bose-Einstein distribution $n_B(k) \equiv 1/(e^{k/T} - 1)$. We can thus expand the distribution functions up to fourth order in $\lambda$ and $\kappa$:

$$ f_{\phi}(k) = n_{B}(k) + f_{\phi}^{(2,2)}(k), \quad (10) $$

$$ f_h(k) = 0 + f_h^{(2,2)}(k) + f_h^{(0,4)}(k), \quad (11) $$

where the superscript stands for the order of $\lambda$ and $\kappa$ that is considered, i.e., $O(\lambda^i \kappa^j) = (i, j)$. Note that in our expansion of the distribution functions, we have also implicitly expanded the matrix element squared. As $\kappa$ is

dimensionful, the expansion in $\kappa \sim 1/m_p$ has to be understood as an expansion in $T\kappa$, the corresponding dimensionless quantity. The zeroth order term of $f_\phi$ is set to be $n_B$, and the zeroth order term of the $f_h$ distribution function is zero. We have suppressed the time arguments in Eqs. (10) and (11). Note that in the expansion, we treat $\lambda$ and $T\kappa$ on equal footing, and for the distribution functions, we have only written out the nonvanishing terms up to fourth order. Terms of $O(\lambda^n \kappa^2)$ with $n > 2$ are also expected, which—depending on the temperature and the value for $\lambda$—can be larger than $O(\kappa^4)$. However, as we discuss later in this section, these terms cannot in general be included in a straightforward manner into the Boltzmann-like ansatz, cf. Eqs. (5), (6), and (7).

In the following, we discuss based on three examples why the terms shown in Eqs. (10) and (11) are the only nonzero contributions. The evolution Eqs. (5) and (6) yield $\dot{f}_h^{(2,0)} = 0$ since, in order to produce or annihilate a graviton, one has to go at least to second order in $\kappa$. Furthermore, $\dot{f}_h^{(0,2)} = 0 = \dot{f}_\phi^{(0,2)}$, since massless $1 \leftrightarrow 2$ processes are kinematically forbidden. Finally, $\dot{f}_\phi^{(2,0)}$ vanishes because of detailed balance arguments. As an example, consider the two terms $\int_{\phi\phi \to \phi\phi}^{(2,0)} n_B n_B(1+n_B)(1+n_B(k)) - \int_{\phi\phi \to \phi\phi}^{(2,0)} n_B(k)n_B(1+n_B)(1+n_B)$ that appear in $\dot{f}_\phi^{(2,0)}$. After a redefinition of variables in the phase space integral, one can show that both terms cancel each other. Similar arguments hold for the other terms in $\dot{f}_\phi^{(2,0)}$ such that overall $\dot{f}_\phi^{(2,0)} = 0$. Note that if $\dot{f}^{(i,j)} = 0$, then $f^{(i,j)} = 0$ for all times for $(i, j) \neq (0, 0)$, which follows from the initial conditions: $f_h(t=0) = 0$ and $f_\phi(t=0) = n_B$. In complete analogy to the discussed examples, one can show with kinematic and detailed balance arguments that the other terms that are not shown in Eqs. (10) and (11) also vanish.

Next we discuss the nonzero fourth order terms for the graviton distribution function. The $\dot{f}_h^{(2,2)}$ rate is nonzero, and it is given by

<footnote n="1">¹Although we consider the case of an initial vanishing distribution function of gravitons, our framework is more generic, as it is valid as long as $f_h \ll 1$.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Lowest order Feynman diagrams for the $\phi\phi \to \phi\phi h$ process. The matrix element squared is of the order $O(\lambda^2\kappa^2)$. In addition to the diagrams shown, there are two more diagrams that can be obtained by crossing the final state graviton.</figcaption></figure>

$$ \dot{f}_h^{(2,2)} = \int_{\phi\phi \to \phi\phi h}^{(2,2)} n_B n_B (1 + n_B)(1 + n_B) + \dots, \quad (12) $$

where the dots stand for other processes of the same order, e.g., $\phi\phi^\dagger\phi \to \phi h$. All possible processes are written down in Sec. IV. In Fig. 2, we show the Feynman diagrams for the single graviton production process $\phi\phi \to \phi\phi h$ at order $O(\lambda^2\kappa^2)$.

The contributions to $\dot{f}_h$ at order $O(\kappa^4)$ are sourced by processes that have two gravitons in the final state. These processes can be relevant in high-temperature early Universe scenarios since the dimensionless expansion parameter is $T/m_p$, which can be relatively large if the temperature is close to the Planck mass. The explicit form of $\dot{f}_h^{(0,4)}$ is

$$ \dot{f}_h^{(0,4)} = \int_{\phi^\dagger\phi \to hh}^{(0,4)} n_B n_B. \quad (13) $$

Processes with a graviton in the initial state do not contribute since the initial state graviton always comes with a factor $f_h$, which makes the whole term of higher order. The same holds for final state $1 + f_h$ amplification factors. The corresponding Feynman diagrams that contribute to two graviton production are shown in Fig. 3. In Sec. IV, we calculate these diagrams.

<figure><img src="image_5.png" /><figcaption>FIG. 3. Lowest order Feynman diagrams for the process $\phi^\dagger\phi \to hh$. The matrix element squared is of the order $O(\kappa^4)$. There is one additional diagram that is not shown since it can be obtained by crossing the final state gravitons in the second diagram.</figcaption></figure>

In the following, we discuss how one could extend our calculation of $f_h$ to higher orders. In particular, we point out the limitations and challenges that one would face. The evolution Eq. (5) for the graviton distribution function can definitively be used without problems at lowest order in perturbation theory; i.e., in our case, these are the $2 \to 2$ and $2 \leftrightarrow 3$ processes that are of the order $O(\kappa^4)$ and $O(\lambda^2\kappa^2)$, respectively. These lowest order terms have real corrections and virtual corrections. Virtual corrections are loop correction, while real corrections come from tree-level processes with extra initial- or final-state particles. If they are finite, the latter are easily incorporated into the Boltzmann equation formalism, i.e., Eq. (7). Incorporating the former, on the other hand, is not straightforward as the matrix elements squared contain not only standard vacuum fluctuations but also statistical fluctuations, which, in turn, depend on the distribution functions themselves [45]. Furthermore, while real and renormalized virtual corrections might separately be finite in a standalone scalar theory, in more complex systems such as the scalar theory coupled to gravity or gauge theories, they are in general not finite, with infrared (IR) divergences canceling between the two, as in the case of the Kinoshita-Lee-Nauenberg theorem [46,47]. In conclusion, it is a challenging task to incorporate higher order effects with the Boltzmann-like approach since it is only possible to incorporate the finite higher order effects.

Virtual gravitons arise already at order $O(\kappa^4)$ in the graviton distribution function, cf. Fig. 3. Quantum gravity effects start to play a role at order $O(\lambda^2\kappa^4)$ since at this order, diagrams with graviton loops exist. The three diagrams that we show in Fig. 4 are of order $O(\lambda^1\kappa^1)$, $O(\lambda^1\kappa^3)$ and $O(\lambda^1\kappa^2)$, respectively. The interference term of the first two diagrams is of order $O(\lambda^2\kappa^4)$ and is the first virtual correction involving loops of gravitons. Conversely, the square of the third diagram is $O(\lambda^2\kappa^4)$ and is part of the real corrections at that order. This further exemplifies the challenge in going beyond leading order: The real corrections can be dealt with in a Boltzmann form in a rather straightforward way, while the virtual corrections cannot, as their matrix element squared will depend in nontrivial ways on the statistical factors—see Ref. [48] for a recent work on this problem in a nongravitational setting.

In an alternative way, one could systematically study quantum effects in kinetic theory from first principles using the Wigner-function formalism. By performing an expansion in $\hbar$ of the Wigner function, one can in principle derive quantum corrections to the classical Boltzmann equation, see, e.g., Refs. [49–53]. The development of such quantum kinetic theory is left for future studies.

Along the lines that we discussed before, we want to mention that it is possible to calculate $\dot{f}_h^{(\text{all},2)}$ in a full quantum picture. As shown in Ref. [27] for the specific case of GWs and more generally in Ref. [54] for any state that is feebly coupled to a thermal bath, the Boltzmann-equation
<figure><img src="image_6.png" /><figcaption>FIG. 4. The interference of the left and center diagrams yields a term $O(\lambda^2\kappa^4)$ in the matrix element squared. The right diagram denotes a real correction that is on the same order when it's squared.</figcaption></figure>

based approach that we use here agrees, at $O(\lambda^2\kappa^2)$, with the thermal-field-theoretical approach of production and equilibration rates. Namely, for GWs, one has $\dot{f}_h^{(\text{all},2)}(k) = \Gamma(k)[n_B(k) - f_h^{(\text{all},2)}(k)]$, [26] where the production/equilibration rate $\Gamma(k)$ is proportional to $\kappa^2$ times the imaginary part of the retarded two-point function of the $T_{12}$ component of the energy-momentum tensor of the equilibrium particles, i.e., in our case, the scalars. This formalism defines the single-graviton production rate to all orders in $\lambda$. Within this formalism, higher orders in $\lambda$ naturally incorporate both real and virtual corrections, without the issues that would plague direct attempts in the Boltzmann-like approach. However note that, in principle, we do not want to go to higher orders in $\lambda$ but to higher orders in $\kappa$ to identify quantum gravity effects. The discussed thermal-field-theoretical approach is not suited for this, and new strategies have to be developed for a full quantum treatment of the graviton production rate.

Finally let us discuss backreaction effects, which can be incorporated into the Boltzmann-like formalism. If we stay at order $O(\kappa^4)$ and go to nonzero order in $\lambda$, we can identify backreaction effects. These appear at lowest order at $O(\kappa^4\lambda^4)$. The graviton production rate at this order contains the following backreaction terms:

$$ \dot{f}_h^{(4,4)}(k) = \int_{\phi\phi\to\phi\phi h}^{(2,2)} f_{\phi}^{(2,2)} n_B(1+n_B)(1+n_B) + \int_{\phi\phi^{\dagger}\to hh}^{(0,4)} f_{\phi}^{(4,0)} n_B + \dots, \quad (14) $$

where the dots above stand for other terms that we have omitted here. We call the terms in Eq. (14) backreaction terms since the small corrections on top of the Bose Einstein distribution, $f_{\phi}^{(2,2)}$ and $f_{\phi}^{(0,4)}$, appear in the phase space integral. Backreaction effects also appear in the $\dot{f}_{\phi}$ rate.

We further note that the RHS in Eqs. (12) and (13) are time independent. Therefore, $f_h^{(2,2)}$ and $f_h^{(0,4)}$ are linear in time, and the backreaction rates, i.e., $\dot{f}_h^{(4,4)}$, are linear in time. From this follows that $f_h^{(4,4)}$ has a quadratic time dependence.

## IV. MATRIX ELEMENTS AND PHASE SPACE INTEGRALS

In this section, we calculate the matrix elements squared for graviton production at order $O(\lambda^2\kappa^2)$ and $O(\kappa^4)$. Let us start with the $O(\lambda^2\kappa^2)$ component. As argued previously, it arises from $2 \to 3$ and $3 \to 2$ processes. The corresponding expressions for the distribution functions are

$$ \dot{f}_h^{2 \to 3}(k) = \frac{1}{16k} \int d\Omega_{2 \to 3} \sum_{abcd} |M_{cdh}^{ab}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{p}_1, \mathbf{p}_2, \mathbf{k})|^2 \\ \times n_B(p'_1) n_B(p'_2) [1 + n_B(p_1)][1 + n_B(p_2)], \quad (15) $$

$$ \dot{f}_h^{3 \to 2}(k) = \frac{1}{24k} \int d\Omega_{3 \to 2} \sum_{abcd} |M_{dh}^{abc}(\mathbf{p}'_1, \mathbf{p}'_2, \mathbf{p}'_3; \mathbf{p}_1, \mathbf{k})|^2 \\ \times n_B(p'_1) n_B(p'_2) n_B(p'_3) [1 + n_B(p_1)], \quad (16) $$

$$ \dot{f}_h^{(2,2)}(k) = \dot{f}_h^{2 \to 3}(k) + \dot{f}_h^{3 \to 2}(k), \quad (17) $$

where at order $O(\lambda^2\kappa^2)$, there is no $f_h$ on the right-hand side. The sums run over all $abcd$ scalar and antiscalar degrees of freedom and thus over all $ab \to cdh$ and $abc \to dh$ processes, with $h$ denoting the graviton. The quantities $|M_{cdh}^{ab}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{p}_1, \mathbf{p}_2, \mathbf{k})|^2$ and $|M_{dh}^{abc}(\mathbf{p}'_1, \mathbf{p}'_2, \mathbf{p}'_3; \mathbf{p}_1, \mathbf{k})|^2$ are the corresponding matrix elements squared summed over the graviton polarizations. For $k \sim T$, the contribution of the thermal mass $m_{\phi T} = \sqrt{\lambda/12}T$ is suppressed, so the external states can be considered massless. The prefactor $1/(16k)$ is a combination of $1/(2k)$ from the phase space measure, $1/2$ for the graviton polarization degeneracy, and $1/(2!)^2$ for the symmetry factors for identical initial and final state particles. In the cases where $a \neq b$ or $c \neq d$, the sum over $abcd$ counts the process two times and compensates for this factor. Similarly, $1/(24k)$ is a combination of $1/(2k)$ from the phase space measure, $1/2$ for the graviton polarization degeneracy, and $1/(3!)^2$ for the symmetry factors for identical initial state particles.

The phase spaces can be read off from Eq. (8). For the matrix element squared, we used the automated pipeline introduced in Ref. [27]. We first used FeynRules [55] to derive Feynman rules for the Lagrangian in Eq. (3). Using the appropriate interface [56], FeynRules generates a model file
for FeynArts [57]. This package and its companion FormCalc [58] were then used to generate, evaluate and square all amplitudes. Tensor boson polarization sums had to be implemented following the method discussed in Ref. [27]. For $\phi\phi^\dagger \to \phi\phi^\dagger h$, the result is

$$ |\mathcal{M}_{\phi\phi^\dagger h}^{\phi\phi^\dagger}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{p}_1, \mathbf{p}_2, \mathbf{k})|^2 = \frac{\kappa^2 \lambda^2}{2} \left[ \frac{(P'_1 \cdot P'_2)^2}{P'_1 \cdot KP'_2 \cdot K} + \frac{(P_1 \cdot P_2)^2}{P_1 \cdot KP_2 \cdot K} - \frac{(P'_1 \cdot P_2)^2}{P'_1 \cdot KP_2 \cdot K} - \frac{(P_1 \cdot P'_2)^2}{P_1 \cdot KP'_2 \cdot K} - \frac{(P'_1 \cdot P_1)^2}{P'_1 \cdot KP_1 \cdot K} - \frac{(P'_2 \cdot P_2)^2}{P'_2 \cdot KP_2 \cdot K} - 2 \right]. \quad (18) $$

Equation (18) arises from diagrams such as the ones in Fig. 2. The four structures in the denominator, e.g., $P'_1 \cdot K$, correspond to the propagator of the virtual, intermediate scalar connecting the $\phi\phi h$ vertex with the $\phi^4$ one. It is reassuring to see that this matrix element squared is finite even when one of the scalar products in the denominators vanishes. For example, the term $P'_1 \cdot K$ can vanish either for $p'_1 \to 0$ or when $\mathbf{p}'_1$ is parallel to $\mathbf{k}$. In the former case, the powers of $p'_1$ at the numerator immediately remove the divergence, and similarly the phase space is free of endpoint divergences for $p'_1 \to 0$, even in the presence of Bose enhancement ($n_B(p'_1 \to 0) \approx T/p'_1$). In the collinear case, $\mathbf{p}'_1 \| \mathbf{k}$, one can show that the sum of the divergent terms is finite.

Let us now discuss the terms that are obtained by crossing. By crossing an initial state $\phi$ ($\phi^\dagger$) to the final state, one obtains a $\phi^\dagger$ ($\phi$). Hence, one finds

$$ \begin{aligned} |\mathcal{M}_{\phi\phi^\dagger h}^{\phi\phi^\dagger}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{p}_1, \mathbf{p}_2, \mathbf{k})|^2 &= |\mathcal{M}_{\phi\phi h}^{\phi\phi}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{p}_1, \mathbf{p}_2, \mathbf{k})|^2 \\ &= |\mathcal{M}_{\phi^\dagger\phi^\dagger h}^{\phi^\dagger\phi^\dagger}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{p}_1, \mathbf{p}_2, \mathbf{k})|^2. \end{aligned} \quad (19) $$

Furthermore, $|\mathcal{M}_{\phi\phi^\dagger h}^{\phi\phi^\dagger}|^2$ is symmetric under permutations within the initial and final states of the $\phi$ and $\phi^\dagger$. Hence, in the sum over $abcd$ of Eq. (15), it is counted four times, whereas the $\phi$-only or $\phi^\dagger$-only processes are counted once.

The $3 \to 2$ matrix elements squared can be obtained by crossing Eq. (18), too. For instance, if we cross the final-state $\phi^\dagger$ with momentum $P_2$ into an initial-state $\phi$ with momentum $P'_3$, we have

$$ |\mathcal{M}_{\phi h}^{\phi\phi^\dagger\phi}(\mathbf{p}'_1, \mathbf{p}'_2, \mathbf{p}'_3; \mathbf{p}_1, \mathbf{k})|^2 = \frac{\kappa^2 \lambda^2}{2} \left[ \frac{(P'_1 \cdot P'_2)^2}{P'_1 \cdot KP'_2 \cdot K} + \frac{(P'_1 \cdot P'_3)^2}{P'_1 \cdot KP'_3 \cdot K} + \frac{(P'_2 \cdot P'_3)^2}{P'_2 \cdot KP'_3 \cdot K} - \frac{(P_1 \cdot P'_1)^2}{P_1 \cdot KP'_1 \cdot K} - \frac{(P_1 \cdot P'_2)^2}{P_1 \cdot KP'_2 \cdot K} - \frac{(P_1 \cdot P'_3)^2}{P_1 \cdot KP'_3 \cdot K} - 2 \right]. \quad (20) $$

By further crossing, one finds

$$ |\mathcal{M}_{\phi^\dagger h}^{\phi^\dagger\phi^\dagger}(\mathbf{p}'_1, \mathbf{p}'_2, \mathbf{p}'_3; \mathbf{p}_1, \mathbf{k})|^2 = |\mathcal{M}_{\phi h}^{\phi\phi^\dagger\phi}(\mathbf{p}'_1, \mathbf{p}'_2, \mathbf{p}'_3; \mathbf{p}_1, \mathbf{k})|^2. \quad (21) $$

The processes with two $\phi$ and with two $\phi^\dagger$ in the initial state are counted three times each in the sum over $abcd$. Putting everything together, we then find

$$ \dot{f}_h^{2 \to 3}(k) = \frac{3}{8k} \int d\Omega_{2 \to 3} |\mathcal{M}_{\phi\phi h}^{\phi\phi}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{p}_1, \mathbf{p}_2, \mathbf{k})|^2 n_B(p'_1) n_B(p'_2) [1 + n_B(p_1)] [1 + n_B(p_2)], \quad (22) $$

$$ \dot{f}_h^{3 \to 2}(k) = \frac{1}{4k} \int d\Omega_{3 \to 2} |\mathcal{M}_{\phi h}^{\phi\phi^\dagger\phi}(\mathbf{p}'_1, \mathbf{p}'_2, \mathbf{p}'_3; \mathbf{p}_1, \mathbf{k})|^2 n_B(p'_1) n_B(p'_2) n_B(p'_3) [1 + n_B(p_1)]. \quad (23) $$

Details on the seven-dimensional numerical integration of the phase space can be found in the Appendix. The $O(\kappa^4)$ component is sourced by $2 \to 2$ processes:

$$ \dot{f}_h^{(0,4)}(k) = \frac{1}{8k} \int d\Omega_{2 \to 2} \sum_{ab} |\mathcal{M}_{hh}^{ab}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{k}, \mathbf{k}_2)|^2 n_B(p'_1) n_B(p'_2), \quad (24) $$
where 1/(8k) is the product of the 1/(2k) from the Lorentz phase space measure, a factor of 1/2 for the two polarizations of the graviton and a factor of 1/(2!) for possible identical initial state particles.

The matrix element squared arises from four diagrams. Three of them are shown in Fig. 3, and the fourth comes from the u-channel analog of the second diagram. We have computed the matrix element squared by using the previously described FeynRules, FeynArts and FormCalc machinery:

$$ |\mathcal{M}_{hh}^{\phi\phi^\dagger}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{k}, \mathbf{k}_2)|^2 = \frac{\kappa^4 t^2 u^2}{8 s^2}, \quad (25) $$

where s, t and u are the standard Mandelstam invariants. Note that the matrix element squared can also be extracted from Refs. [59,60]<fnref n="2" /> and we have checked that our results agree. Accounting for a factor of 2 from the sum in Eq. (24) yields:

$$ \dot{f}_h^{(0,4)}(k) = \frac{1}{4k} \int d\Omega_{2\to} |\mathcal{M}_{hh}^{\phi\phi^\dagger}(\mathbf{p}'_1, \mathbf{p}'_2; \mathbf{k}, \mathbf{k}_2)|^2 \times n_B(\mathbf{p}'_1) n_B(\mathbf{p}'_2). \quad (26) $$

We refer again to the Appendix for details on the reduction of the phase space integration to a two-dimensional integral that we evaluate numerically.

The production rates for $(i, j) = (2, 2)$ and $(i, j) = (0, 4)$ can be written in a compact form as

$$ \dot{f}_h^{(i,j)}(k) = \frac{1}{2k} T^2 \left(\frac{T}{m_p}\right)^j \lambda^i n_B\left(\frac{k}{T}\right) \psi^{(i,j)}\left(\frac{k}{T}\right), \quad (27) $$

where we have defined dimensionless $\psi$ functions and used the convention that $n_B(x) = 1/(e^x - 1)$ if the argument of $n_B$ is dimensionless. The $\psi$ functions are shown in Fig. 5 for the $2 \to 3$ (dotted black) and $3 \to 2$ (dashed black) processes. We also show the sum of the $2 \to 3$ and $3 \to 2$ processes, which is labeled as $2 \leftrightarrow 3$ and shown as a red dot-dashed line. The $2 \to 2$ processes are shown as a solid blue line. Note that while the $\psi$ function for the $2 \leftrightarrow 3$ processes has only a relatively mild $k/T$ dependence around $k/T \simeq 1$, this is not the case for the $2 \to 2$ processes. In the Appendix, we have derived an asymptotic form for $\dot{f}_h^{(0,4)}$, i.e., $\psi^{(0,4)}$ in the limit $k > T$. We find

<figure><img src="image_7.png" /><figcaption>FIG. 5. We plot the $\psi$ functions that have been obtained by numerically integrating the phase space integrals. The results are shown for the $2 \to 3$, $3 \to 2$ and $2 \to 2$ processes, respectively. The sum of the $2 \to 3$ and $3 \to 2$ contributions is denoted as $2 \leftrightarrow 3$. $2 \leftrightarrow 3$ processes come from single graviton production processes that have a lowest order matrix element squared of the order $\mathcal{O}(\lambda^2\kappa^2) = (2, 2)$. $2 \to 2$ processes arise in graviton pair production processes with a lowest order matrix element squared of the order $\mathcal{O}(\kappa^4)$.</figcaption></figure>

$\psi^{(0,4)} = 32/(15\pi)(k/T)^2$ for $k > T$ which is in agreement with Fig. 5. $\psi^{(2,2)}$ and $\psi^{(0,4)}$ differ in their functional form due to the fact that the underlying phase space structure and matrix elements are different for single and graviton pair production. The single gravitons are produced in $2 \leftrightarrow 3$ processes, while graviton pair production is a $2 \to 2$ process. In addition, the nonrenormalizable nature of gravity leads to a matrix element squared for single graviton production that is proportional to $\kappa^2 \sim \frac{1}{m_p^2}$ times an expression with dimensions (energy)$^2$. Similarly the matrix element squared for graviton pair production is proportional to $\kappa^4 \sim \frac{1}{m_p^4}$ times an expression with dimensions (energy)$^4$. In the case of graviton pair production, the $2 \to 2$ phase space integration for large $k$ translates this into a quadratic dependence on momentum, which explains the form of $\psi^{(0,4)}$ for large $k$.

As Fig. 5 shows, the $2 \to 3$ contribution to $\dot{f}_h^{(2,2)}$ leads to $\dot{f}_h^{(2,2)} \propto k^{-3}$ for $k \ll T$, which implies a naively IR-divergent contribution to the number density of gravitons ($n_{gw} \propto \int dk k^2 f_h$) and a finite but enhanced contribution to the energy density $\rho_{gw} \propto \int dk k^3 f_h$ from the IR domain $k \ll T$. This IR contribution is an artifact of treating the external and intermediate scalar states as massless. If we would include their thermal mass $m_{\phi T}$, then such a behavior would go away, as the matrix element would no longer behave like $1/k^2$ at small $k$. Scalars have the nice property that their thermal mass behaves like an ordinary local mass term in the Lagrangian, unlike gauge fields and fermions.

<footnote n="2">²These papers show that the amplitudes for graviton production factorize into simple products of photon amplitudes times kinematic factors. The massless limit of Eq. (21) of Ref. [59] gives the $\phi\gamma \to \phi\gamma$ scalar Compton amplitudes in the helicity basis. Equation (62) of the same paper then gives the matrix elements squared as the fourth power of these amplitudes, multiplied by the second power of the kinematic factor $F$, given in Eq. (61). The matrix elements squared for the two polarizations are identical and thus trivially summed. Finally, crossing symmetry relates the $\phi h \to \phi h$ matrix element squared to the $\phi^\dagger h \to hh$ matrix element squared.</footnote>
This in turn would make thermal mass resummation relatively straightforward. For this paper, we limit ourselves to unresummed (massless) results and consider the result for $\dot{f}_h$ to be a proper leading order determination in the regime $k \gtrsim m_{\phi T} = \sqrt{\lambda/12T}$. Furthermore, for smaller $k$, $k \ll \lambda^2 T$, the quasiparticle description breaks down completely, and gravitational waves are sourced from hydrodynamic fluctuations [26].

## V. GRAVITATIONAL WAVE SPECTRUM

In this section, we embed the graviton production rate into cosmological evolution. Our main assumption is that, after the hot big bang, a thermal plasma of $\phi$ particles with temperature $T_{max}$ is present. Throughout the expansion of the Universe, the thermal plasma produces GWs. From the definition of Eq. (4), it follows that the GW differential energy density is $d\rho_{gw}(t, k) = 2kf_h(t, k) \frac{d^3k}{(2\pi)^3}$, where a flat space-time metric has been assumed, and the factor of 2 takes the two polarization states into account which contribute to the energy density. We can rewrite the equation for the energy density as $\frac{d\rho_{gw}}{dt d\ln k} = \frac{k^4}{\pi^2} \dot{f}_h$. Generalizing to an expanding Universe, the GW energy density evolves as [26]

$$ (\partial_t + 4H)\rho_{gw}(t) = \int \frac{d^3k}{(2\pi)^3} R(t, k), \quad (28) $$

where $H$ is the Hubble parameter, and we have defined $R(t, k) \equiv 2kf_h(t, k)$. Note that for the $O(\lambda^2 k^2)$ and $O(k^4)$ contributions, which we shall discuss here, $\dot{f}_h$ has no explicit time dependence. We will therefore treat $\dot{f}_h$ without explicit time dependence in the following derivation. Now that we consider a thermal plasma in an expanding Universe, the temperature decreases over time. Therefore, we have an implicit time/temperature dependence.

We further note that, in a radiation-dominated Universe,

$$ H = \sqrt{\frac{8\pi\rho(T)}{3}} \frac{1}{m_p}, \quad (29) $$

with $\rho(T) = g_{*\rho}(T)\pi^2 T^4/30$, and $g_{*\rho}(T) = 2$ is the effective number of energy density degrees of freedom. The scalars remain in thermal equilibrium<fnref n="3" /> as long as their interaction rate, which is on the order $\lambda^2 T$, is at least as fast as the Hubble rate $H \sim T^2/m_p$. Therefore, we obtain the equilibrium condition, $\lambda^2 \gtrsim T/m_p$, which has to be understood as an order of magnitude estimate.

We can integrate Eq. (28):

$$ \frac{\rho_{gw}(t_1)}{s^{4/3}(t_1)} - \frac{\rho_{gw}(t_0)}{s^{4/3}(t_0)} = \int_{t_0}^{t_1} dt \frac{1}{s^{4/3}(t)} \int \frac{d^3 k}{(2\pi)^3} R(T(t), k), \quad (30) $$

where we have used that the entropy density $s$ fulfills $\dot{s} + 3Hs = 0$, and we have made the temperature/time dependence explicit. We assume that at the beginning of the GW production, no GWs are present, i.e., $\rho_{gw}(t_0) = 0$. At $t_0$, the thermal plasma was first in thermal equilibrium and had its maximum temperature $T_{max}$. The time $t_1$ is the time when the mass of the scalar field cannot be neglected anymore. The temperature corresponding to the time $t_1$ is referred to as $T_\phi$ in the following. In Eq. (30), we can only integrate to $t_1$ since the production rates that we have calculated are only valid for temperatures above $T_\phi$. The time integral in Eq. (30) can be transformed into an integral over the temperature by using the relation [61]:

$$ \frac{dT}{dt} = -\sqrt{\frac{4\pi^3}{45}} g_{*\rho}(T) \frac{1}{2} \frac{g_{*s}(T)}{g_{*c}(T)} \frac{T^3}{m_p}, \quad (31) $$

where $g_{*s}$ and $g_{*c}$ are the effective degrees of freedom for the entropy density and heat capacity, which are defined as $s(T) \equiv g_{*s}(T) \frac{2\pi^2}{45} T^3$ and $c(T) \equiv g_{*c}(T) \frac{2\pi^2}{15} T^3$. In the following, we use the assumption of isotropy under which we can simplify the $d^3k$ integral: $\int d^3k R(T, k) = 4\pi \int dk k^2 R(T, k) = 4\pi \int d\ln(k) k^3 R(T, k)$. From Eq. (30), we can then read off the GW energy density per logarithmic momentum interval at $T_\phi$, normalized to the total energy density: $\Omega_{gw} = \frac{1}{\rho} \frac{d\rho_{gw}}{d\ln k}$. Redshifting all corresponding quantities to today [28] yields

$$ h_0^2 \Omega_{gw} (f_g) = \frac{15\sqrt{45}}{4\pi^{11/2}} m_p g_{*s} (T_{\text{today}})^{1/3} h_0^2 \Omega_\gamma \left( \frac{2\pi f_g}{T_{\text{today}}} \right)^3 \\ \times \int_{T_\phi}^{T_{\text{max}}} dT \frac{1}{T^4} \frac{g_{*c}(T)}{g_{*\rho}(T)^{1/2} g_{*s}(T)^{4/3}} \\ \times R \left( T, T \frac{2\pi f_g}{T_{\text{today}}} \left( \frac{g_{*s}(T)}{g_{*s}(T_{\text{today}})} \right)^{1/3} \right), \quad (32) $$

where $f_g$ is the current day GW frequency, $h_0^2 \Omega_\gamma = 2.473 \times 10^{-5}$ is the present fractional photon energy density, $h_0^2$ a factor that eliminates the experimental uncertainty that is coming from measurements of the Hubble constant, $T_{\text{today}} = 2.7254$ K is the current day temperature [62] and $g_{*s}(T_{\text{today}}) = 3.931$ are the effective entropy degrees of freedom today [63].

With the parametric form of $\dot{f}_h$ from Eq. (27), we can write $R$ as

<footnote n="3">³By this, we mean that the zeroth order term of $f_\phi$ is a massless Bose-Einstein distribution with the current temperature.</footnote>
$$ 
\begin{aligned}
R(T, k) &= 2k \dot{f}_h(T, k) \\
&= T^2 n_B \left(\frac{k}{T}\right) \frac{T^2}{m_p^2} \left(\lambda^2 \psi^{(2,2)}\left(\frac{k}{T}\right) + \frac{T^2}{m_p^2} \psi^{(0,4)}\left(\frac{k}{T}\right) + \dots \right),
\end{aligned} 
\quad (33) $$

where the dots denote higher order terms. We plug Eq. (33) in Eq. (32) and, in order to get an analytical result, we approximate $g_{*s}(T) = g_{*\rho}(T) = g_{*c}(T) = g_{*}(T_{max})$ in the region of temperatures above $T_{\phi}$. We thus obtain

$$ h_0^2 \Omega_{gw}(f_g) = 5.54 \times 10^{-12} \left(\frac{f_g}{10^{10} \text{ Hz}}\right)^3 \times \left(\frac{2}{g_*(T_{\text{max}})}\right)^{\frac{5}{6}} \left(\frac{T_{\text{max}}/m_p}{10^{-3}}\right) n_B(y_{\text{max}}) \times \left(\lambda^2 \psi^{(2,2)}(y_{\text{max}}) + \frac{1}{3} \left(\frac{T_{\text{max}}}{m_p}\right)^2 \times \psi^{(0,4)}(y_{\text{max}}) + \dots \right), \quad (34) $$

where we have assumed $T_{max} \gg T_{\phi}$, and we have defined $y_{max} = \frac{2\pi f_g}{T_{\text{today}}} \left(\frac{g_{*s}(T_{max})}{g_{*s}(T_{\text{today}})}\right)^{1/3} = 0.14(\frac{f_g}{10^{10} \text{ Hz}})\left(\frac{g_{*s}(T_{max})}{2}\right)^{1/3}$. Note that models that can describe the entire thermal history of the early Universe have $g_{*}(T_{max}) > g_{*s}(T_{\text{today}})$. We work with a model that includes only one complex scalar field, and therefore, $g_{*}(T_{max}) < g_{*s}(T_{\text{today}})$. Nonetheless the

features in the GW spectrum that we work out here will hold in general even with a more realistic model that can describe the thermal history of the Universe consistently. We comment on this aspect further below. The terms that are represented by the dots in Eq. (34) include processes which encode quantum gravity effects, cf. Fig. 4. These effects arise at the order $O(\lambda^2\kappa^4)$ and would appear as a term $\psi^{(2,4)}\lambda^2(T_{max}/m_p)^2$ in the parentheses in the second line in Eq. (34).

The single graviton production processes have a matrix element squared of the order $O(\lambda^2\kappa^2)$, i.e., it is proportional to $1/m_p^2$. Since the GW production happens on a timescale that is comparable with $m_p$, the single graviton production processes are proportional to $T_{max}/m_p$ in the GW spectrum. A similar argument applies to the contribution from the graviton pair production. The matrix element squared is of the order $O(\kappa^4)$, i.e., proportional to $1/m_p^4$, and hence, the contribution to the GW spectrum is of the order $(T_{max}/m_p)^3$. We refer to this production mechanism as GW freeze-in, since the GWs are produced from the thermal plasma throughout the expansion of the Universe. Note that while the GW production mechanism is conceptually very similar to the production of dark matter from the thermal plasma, the final GW spectrum is very ultraviolet sensitive in the sense that it depends on the maximum temperature.

In Fig. 6, we plot the GW spectrum coming from single graviton production processes (red dot-dashed lines) and from graviton pair production processes (blue dashed

<figure><img src="image_8.png" /><figcaption>FIG. 6. Gravitational wave spectrum with respect to the present day GW frequency. The single graviton production processes that are of order $O(\lambda^2\kappa^2) = (2, 2)$ are shown as a red dot-dashed line. The graviton pair production processes are shown as a dashed blue line and are of order $O(\kappa^4) = (0, 4)$. The total GW spectrum is shown as a solid black line. On the left, we show a scenario where the maximum temperature is limited to $10^{-3} m_p$. In the right figure, we set the maximum temperature to $T_{max} = 10^{-2} m_p$. In this case, the graviton pair production processes yield an even larger contribution compared to the single graviton production processes. The GW spectra are calculated for a complex scalar model with $g_{*}(T_{max}) = 2$.</figcaption></figure>
lines). The sum of both contributions, i.e., the total GW spectrum, is shown as a black solid line. The quartic coupling is always set to $\lambda = 10^{-1}$.$^4$ When we evaluate the GW spectrum, we also have to evaluate the $\psi$ functions. Since these functions have only been calculated reliably for arguments that are larger than $\sqrt{\lambda} = 0.31$, cf. Sec. IV, we only show the GW spectrum in the corresponding frequency regime. Note that the GW spectrum from the single graviton production mimics to a good approximation a black body spectrum since the function $\psi^{(2,2)}$ is very flat in the regime $T/m_p > \sqrt{\lambda}$. The graviton pair production contribution has a significantly different shape since the function $\psi^{(0,4)}$ is not constant in the frequency interval of interest.

Figure 6 (left) shows a scenario where the maximum temperature is set to $T_{max}/m_p = 10^{-3}$. In this case, the single graviton production contribution dominates over the graviton pair production contribution. The total spectrum has therefore mostly the form of the single graviton production spectrum.

Figure 6 (right) shows the GW spectrum for a slightly larger maximum temperature, $T_{max}/m_p = 10^{-2}$, which requires a noninflationary scenario. Note that the chosen maximum temperature is consistent with the equilibrium condition for the scalar fields, i.e., $\lambda^2 \gtrsim T/m_p$. In the case at hand, both contributions are parametrically equally important, and around the peak frequency, the graviton pair production contribution is even substantially larger than the single graviton production contribution. This can be seen explicitly from Eq. (34) by comparing the two terms in the second line. The first term corresponds to single graviton production processes, while the second one describes the contribution due to graviton pair production. Comparing both terms, we find as an order of magnitude estimate that the contribution from graviton pair production processes are equally important or even larger than the contribution from the single graviton production processes if $10T_{max}/m_p \gtrsim \lambda$. Therefore, the relevance of the graviton pair production processes depends crucially on the size of the coupling and the maximum temperature. The GW spectra associated with single graviton and graviton pair production processes peak at slightly different frequencies and have a distinct functional form. As a result, the total spectrum takes a very characteristic form that is substantially different from the single graviton production spectrum, i.e., an approximate black body spectrum.

In the following, we derive analytic expressions for the peak frequencies of the single graviton and graviton pair production GW spectra. For the single graviton contribution, we find from Eq. (34) that the GW spectrum peaks at

$2 \times 10^{11}$ Hz$(2/g_{*s}(T_{max}))^{1/3}$, where we have assumed that $\psi^{(2,2)} = \text{const}$, which is a good approximation in the frequency interval of interest. The peak frequency of the graviton pair production curve lies at a slightly higher frequency: $3.5 \times 10^{11}$ Hz$(2/g_{*s}(T_{max}))^{1/3}$, where we have used the asymptotic form $\psi^{(0,4)} = 32/(15\pi)(k/T)^2$ that was derived in the Appendix.

The higher order terms that are depicted by the dots in Eq. (34) are suppressed further by powers of $T_{max}/m_p$ and $\lambda$. For values of $T_{max}/m_p$ and $\lambda$ that are relatively close to unity, one has to check in detail that the higher order $\psi$ functions are not larger than the leading order $\psi$ functions such that the suppression from the additional powers of $T_{max}/m_p$ and $\lambda$ is not spoiled. The values of $T_{max}$ and $\lambda$ that we consider in Fig. 6 are much smaller than unity, and therefore, we do not expect such a scenario to happen. For example, the $\psi^{(0,6)}$ function would have to be 4 orders of magnitude larger than the $\psi^{(0,4)}$ function for $T_{max}/m_p = 10^{-2}$. The higher order $\psi$ functions will be either phase space suppressed or have the same phase space as the leading order processes. In both cases, we do not expect an enhancement of the higher order $\psi$ functions by orders of magnitude.

The maximum of the CGMB spectrum is bounded from above, $h_0^2\Omega_{gw} \lesssim 10^{-6}$ [28], due to constraints on the additionally allowed amount of dark radiation. Both scenarios that are shown in Fig. 6 do not saturate this bound and are therefore not excluded. In the complex scalar model, the dark radiation bound is saturated for $T_{max}/m_p \approx 0.5$. Note that in this case, the main contribution to the GW spectrum is coming from graviton pair production processes, which illustrates the importance to include these processes at high temperatures. The SM predictions for the GW production from the thermal plasma includes currently only single graviton production processes, cf. Refs. [26–28], and yields $h_0^2\Omega_{gw} = 5 \times 10^{-7}$ for $T_{max} = m_p$. Therefore, adding graviton pair production processes to the SM calculation can already lead to a violation of the dark radiation bound for $T_{max} < m_p$ because the current prediction, which includes only single graviton production processes, is already very close to the dark radiation bound. In conclusion, this might then be used to constrain the maximum temperature of the Universe. We plan to work out the details in a follow-up study where we will also address the aspect of thermal equilibrium at ultra high temperatures in the SM and beyond the SM theories.

## VI. CONCLUSIONS AND OUTLOOK

The thermal plasma in the early Universe produced a guaranteed stochastic GW background through thermal fluctuations. At each time, the emitted GW spectrum peaks at the respective temperature. Due to the temperature-redshift relation, the peak frequencies of the GW spectra are all redshifted to the same frequency today and therefore add up. Conceptually, the GW production from the thermal

<footnote>$^4$In a more realistic model that can describe the entire history of our Universe, one would have to use renormalization group equations to run the parameters up to very high energy scales.</footnote>
plasma has many similarities with the so-called dark matter freeze-in production from the thermal plasma. The GWs are produced out of equilibrium, and their distribution function is small at all times, $f_h \ll 1$. Furthermore, the $f_h$ distribution function evolves much slower than the Hubble rate. We therefore dubbed the GW production from the thermal plasma *GW freeze-in production*. The GW freeze-in scenario is ultraviolet dominated in the sense that it depends on the maximum temperature of the Universe, as expected from a nonrenormalizable coupling.

In this paper, we use a Boltzmann-like formalism to study the microscopic particle collision processes that contribute to the CGMB spectrum. We have done all calculations in a model with a complex scalar field and quartic self-interaction. Our basic assumption is that after the hot big bang, a plasma of scalars with temperature $T_{max}$ is present, and this plasma produced the CGMB spectrum. First, we considered the contribution of single graviton production processes to the CGMB spectrum. In a scalar theory with quartic interaction, single graviton production processes are $2 \leftrightarrow 3$ processes, which have not been calculated before. Our calculation is motivated by the fact that a quartic coupling exists in the Higgs sector of the SM and in many BSM theories. The second class of processes that we investigate are graviton pair production processes. These are $2 \rightarrow 2$ processes and have not been considered before in the context of GWs from the thermal plasma. We show that their contribution to the CGMB spectrum can be larger than the contribution from the single graviton production processes. As an order of magnitude estimate, graviton pair production processes dominate the GW spectrum if $10 T_{max}/m_p \gtrsim \lambda$. Note however that the maximum temperature is also bounded from above by an equilibrium requirement for the scalar particles: $\lambda^2 \gtrsim T_{max}/m_p$, which has to be seen as a parametric estimate. Therefore, the degree to which graviton pair production processes contribute significantly to the CGMB spectrum depends on the values of the coupling coefficient and the maximum temperature. As an example, we show the two different scenarios in Fig. 6. On the left, single graviton production processes dominate ($T_{max}/m_p = 10^{-3}$ and $\lambda = 10^{-1}$). When increasing $T_{max}$ by one order of magnitude (Fig. 6, right), graviton pair production processes yield a significant contribution to the total GW spectrum.

The single graviton and graviton pair production processes are the lowest order contributions, which can easily be incorporated into our Boltzmann-like formalism. We have also discussed the first steps and problems that would arise if one would add real and virtual quantum gravity corrections to the presented results. While finite real corrections can be incorporated in our formalism, virtual corrections depend on the distribution functions themselves, and this complicates their inclusion into the Boltzmann-like approach that we use here. A possible

future direction is to derive a quantum Boltzmann equation from the Wigner function by performing a systematic $\hbar$ expansion. This would allow one to explicitly identify the quantum corrections.

The results that we have worked out for a scalar model are qualitatively also valid for more general theories. In this case, the coupling coefficient $\lambda$ would have to be replaced with the heat bath couplings in the more general theory, which we generically refer to as $g$. Then the contribution from graviton pair production processes to the GW spectrum dominates over the single graviton contribution if $X \times T_{max}/m_p \gtrsim g$, where $X$ is a model dependent constant. A value $X < 10$ in the SM would indicate that graviton pair production dominates at higher temperatures compared to our scalar model. In a follow-up study, we will answer this point with a full SM and BSM calculation. Confronting the graviton pair production calculation in the SM and BSM theories with existing dark radiation constraints can therefore already lead to constraints on $T_{max}$. Constraints on $T_{max}$ can be used to test different models of our Universe, cf. Ref. [28]. For example, nonstandard inflationary cosmological models would be required if a $T_{max} \gtrsim 10^{-3} m_p$ would be inferred from the GW spectrum. Furthermore, the CGMB can be used to constrain nonstandard cosmological histories, cf. Ref. [64]. The authors of Ref. [64] have only considered single graviton production processes. It would be interesting to add graviton pair production processes to their calculation since it could lead to stronger constraints on nonstandard cosmological histories.

A detection of the CGMB with Earth-based detectors will be challenging, cf. Ref. [65]; however, there exist detector proposals with sensitivities comparable to the dark radiation bound, cf. Ref. [66]. Future experimental work will have to show if the proposed detectors can be realized and if their foreseen sensitivity can even be improved. Our results motivate further work on high-frequency GW detection since a detection of the CGMB in the future would pave the way to probe our understanding of particle physics and cosmology at ultra high energies.

## ACKNOWLEDGMENTS

We thank Alex Dima, Yoni Kahn, Mikko Laine, Kaloian Lozanov, Daniel Meuser, Patrick Peter and Andreas Ringwald for useful discussions. We also thank Nick Abboud, Rachel Nguyen and Michael Wentzel for linguistic corrections. The work of J. S. E. is supported in part by DOE Grant No. DE-SC0015655. J. G. acknowledges support by a PULSAR grant from the Région Pays de la Loire and by the Agence Nationale de la Recherche (France), under Grant No. ANR-22-CE31-0018 (AUTOTHERM). J. S. E. would like to express special thanks to the Mainz Institute for Theoretical Physics (MITP) of the Cluster of Excellence PRISMA+ (Project ID 39083149), for its hospitality and support.
APPENDIX: DETAILS ON THE EVALUATION OF THE PHASE SPACE INTEGRALS

In this appendix, we provide some extra details on the phase space integrals of Sec. IV. Let us start from Eqs. (22) and (23). In order to carry out the integrations numerically, we can rewrite the phase space as

$$ \int d\Omega_{2\to3} = -\frac{1}{16(2\pi)^7} \int_0^\infty dp'_1 p'_1 \int_0^\infty dp'_2 p'_2 \int_{-1}^1 dc_{p'_1} \int_{-1}^1 dc_{p'_2} \int_{-1}^1 dc_{p_1} \int_0^{2\pi} d\phi_1 \int_0^{2\pi} d\phi_2 \frac{p_1^2 \theta(p_1) \theta(p'_1 + p'_2 - p_1 - k)}{P_1 \cdot P_2}, \quad (A1) $$

where $k$ is chosen to point in the $z$ direction, and the $c$ variables are the cosines of the angles between $k$ and the respective momenta: $c_p \equiv \cos\theta_{k,p}$. The $\phi$'s are two azimuthal angles, where the third one was integrated out. $p_1$ is fixed to

$$ p_1 = [P'_1 \cdot P'_2 - P'_1 \cdot K - P'_2 \cdot K][k(1 - c_{p_1}) - p'_1(1 - (c_{p'_1} c_{p_1} + \sqrt{(1 - c_{p'_1}^2)(1 - c_{p_1}^2)} \cos(\phi_1 + \phi_2))) - p'_2(1 - (c_{p'_2} c_{p_1} + \sqrt{(1 - c_{p'_2}^2)(1 - c_{p_1}^2)} \cos(\phi_2)))], \quad (A2) $$

and $P'_i \cdot K = p'_i k(c_{p'_i} - 1)$, $P'_1 \cdot P'_2 = p'_1 p'_2 (c_{p'_1} c_{p'_2} + \sqrt{(1 - c_{p'_1}^2)(1 - c_{p'_2}^2)} \cos(\phi_1) - 1)$. Similarly, $P'_i \cdot P_1 = p'_i p_1 (c_{p'_1} c_{p_1} + \sqrt{(1 - c_{p'_1}^2)(1 - c_{p'_1}^2)} \cos(\phi_1 + \phi_2) - 1)$. The other inner products, including $P_1 \cdot P_2$, follow from $P_2 = P'_1 + P'_2 - K - P_1$. The $3 \to 2$ analog of Eq. (A1) follows from simple crossings. These seven-dimensional integrals are carried out numerically using the Monte Carlo algorithm Vegas+ [67]. The results are shown in Fig. 5.

Let us now consider the $O(\kappa^4)$ contribution. Equation (26) can be evaluated using the standard “s-channel” parametrization of Refs. [68,69]. We can arrange the phase space integral as

$$ \int d\Omega_{2\to3} = \frac{1}{(4\pi)^3 k} \int_k^\infty dq^0 \int_{|2k-q^0|}^{q^0} dq \int_{q_-}^{q_+} dp'_2 \int_0^{2\pi} \frac{d\phi_{p'_2 k}}{2\pi}, \quad (A3) $$

where we defined $q_\pm \equiv (q^0 \pm q)/2$, and we chose $q^0, q$ such that $p'_1 = q^0 - p'_2$ and $\mathbf{p}'_1 + \mathbf{p}'_2 = \mathbf{q}$. $\phi_{p'_2 k}$ is the azimuthal angle between the $\mathbf{p}'_2, \mathbf{q}$ and $\mathbf{k}, \mathbf{q}$ planes. This corresponds to

$$ s = q_0^2 - q^2, \quad t = -\frac{s}{2q^2} \left[ (2k - q^0)(q^0 - 2p'_2) + q^2 - \cos(\phi_{p'_2 k}) \sqrt{(q^2 - (q^0 - 2k)^2)(q^2 - (q^0 - 2p'_2)^2)} \right]. \quad (A4) $$

We can perform the angular average $\langle \dots \rangle_{\phi_{p'_2 k}} \equiv \int_0^{2\pi} \frac{d\phi_{p'_2 k}}{2\pi} \dots$ to get rid of odd powers of the cosine and find

$$ \dot{f}_h^{(0,4)} = \frac{\kappa^4}{32k} \frac{1}{(4\pi)^3 k} \int_k^\infty dq^0 \int_{|2k-q^0|}^{q^0} dq \int_{q_-}^{q_+} dp'_2 \left\langle t^2 \left[ \frac{t^2}{s^2} + 2\frac{t}{s} + 1 \right] \right\rangle_{\phi_{p'_2 k}} n_B(q^0) [1 + n_B(p'_2) + n_B(q^0 - p'_2)]. \quad (A5) $$

We have used the identity $n_B(p'_2)n_B(q^0 - p'_2) = n_B(q^0)[1 + n_B(p'_2) + n_B(q^0 - p'_2)]$, which is useful for treating the $p'_2$ integration analytically.<fnref n="5" /> Carrying out the integration, we find

<footnote n="5">⁵The matrix element squared depends on $p'_2$ as function of $q^0 - 2p'_2$. It thus has a reflection symmetry around $p'_2 = q^0/2$, which is the midpoint of the $p'_2$ integration range. We can then reflect the $n_B(q^0 - p'_2)$ into an $n_B(p'_2)$ leaving the matrix element unchanged, further simplifying the integration.</footnote>
$$ \dot{f}_h^{(0,4)} = \frac{\kappa^4}{4(8\pi)^3 k^2} \int_k^\infty dq^0 \int_{|2k-q^0|}^{q^0} dq n_B(q^0)(q_0^2 - q^0)^2 \left\{ -\frac{11q^4 - 30q^2(2k-q^0)^2 + 15(2k-q^0)^4}{120q^3} \right. \\ + T \frac{(q^2 - (q^0 - 2k)^2)^2}{8q^4} \ln \frac{e^{\frac{q_+}{T}} - 1}{e^{\frac{q_+}{T}} - 1} - T^2 \frac{(q^2 - (q^0 - 2k)^2)(q^2 - 5(q^0 - 2k)^2)}{2q^5} \left( \mathrm{Li}_2\left(e^{-\frac{q_+}{T}}\right) + \mathrm{Li}_2\left(e^{-\frac{q_-}{T}}\right) \right) \\ - T^3 \frac{45(2k-q^0)^4 + q^2(5q^2 - 42(q^0-2k)^2)}{2q^6} \left( \mathrm{Li}_3\left(e^{-\frac{q_+}{T}}\right) - \mathrm{Li}_3\left(e^{-\frac{q_-}{T}}\right) \right) \\ - T^4 \frac{9q^4 - 90q^2(2k-q^0)^2 + 105(2k-q^0)^4}{q^7} \left[ \mathrm{Li}_4\left(e^{-\frac{q_+}{T}}\right) + \mathrm{Li}_4\left(e^{-\frac{q_-}{T}}\right) + \frac{2T}{q} \left( \mathrm{Li}_5\left(e^{-\frac{q_+}{T}}\right) - \mathrm{Li}_5\left(e^{-\frac{q_-}{T}}\right) \right) \right] \left. \right\} \quad (A6) $$

Note that the result of the integration is positive, though this might not appear obvious from this expression. For large momenta, $k \gg T$, Eq. (A6) asymptotes to

$$ \dot{f}_h^{(0,4)}|_{k \gg T} = \frac{\kappa^4 T^4}{15(4\pi)^3} e^{-k/T} \left( k + \mathcal{O}\left(\frac{1}{k^2}\right) \right). \quad (A7) $$

This result can be extracted by noting that in this asymptotic regime, $n_B(q^0) \approx e^{-q^0/T}$. This sharp exponential cutoff ensures that only the $q^0 \approx k$ and $2k-q^0 < q < q^0$ ($q \approx k$) ranges dominate the integral. Expanding the integrand for $q^0 - k \ll k$ and $q - k \ll k$ and then performing the integral, we recover Eq. (A7). We note that the form given in Eq. (A7), while valid for $k \gg T$, approximates the numerical results shown in Fig. 5 at better than 30% accuracy for $k > T$.

[1] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), Observation of gravitational waves from a binary black hole merger, Phys. Rev. Lett. **116**, 061102 (2016).

[2] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. **119**, 161101 (2017).

[3] L. P. Grishchuk, Amplification of gravitational waves in an isotropic Universe, J. Exp. Theor. Phys. **40**, 409 (1975) [Zh. Eksp. Teor. Fiz. **67**, 825 (1974)].

[4] A. A. Starobinskii, Spectrum of relict gravitational radiation and the early state of the Universe, JETP Lett. **30**, 682 (1979) [Pis'ma Zh. Eksp. Teor. Fiz. **11**, 719 (1979)].

[5] V. A. Rubakov, M. V. Sazhin, and A. V. Veryaskin, Graviton creation in the inflationary Universe and the grand unification scale, Phys. Lett. **115B**, 189 (1982).

[6] R. Fabbri and M. D. Pollock, The effect of primordially produced gravitons upon the anisotropy of the cosmological microwave background radiation, Phys. Lett. **125B**, 445 (1983).

[7] S. Y. Khlebnikov and I. I. Tkachev, Relic gravitational waves produced after preheating, Phys. Rev. D **56**, 653 (1997).

[8] K. Lozanov, *Reheating After Inflation* (Springer, New York, 2020).

[9] Yohei Ema, Ryusuke Jinno, Kyohei Mukaida, and Kazunori Nakayama, Gravitational effects on inflaton decay, J. Cosmol. Astropart. Phys. 05 (2015) 038.

[10] Yohei Ema, Ryusuke Jinno, Kyohei Mukaida, and Kazunori Nakayama, Gravitational particle production in oscillating backgrounds and its cosmological implications, Phys. Rev. D **94**, 063517 (2016).

[11] Yohei Ema, Ryusuke Jinno, and Kazunori Nakayama, High-frequency graviton from inflaton oscillation, J. Cosmol. Astropart. Phys. 09 (2020) 015.

[12] Edward Witten, Cosmic separation of phases, Phys. Rev. D **30**, 272 (1984).

[13] C. J. Hogan, Gravitational radiation from cosmological phase transitions, Mon. Not. R. Astron. Soc. **218**, 629 (1986).

[14] Thibault Damour and Alexander Vilenkin, Gravitational wave bursts from cosmic strings, Phys. Rev. Lett. **85**, 3761 (2000).

[15] Thibault Damour and Alexander Vilenkin, Gravitational wave bursts from cusps and kinks on cosmic strings, Phys. Rev. D **64**, 064008 (2001).

[16] Arthur Kosowsky, Andrew Mack, and Tinatin Kahniashvili, Gravitational radiation from cosmological turbulence, Phys. Rev. D **66**, 024030 (2002).

[17] Alberto Nicolis, Relic gravitational waves from colliding bubbles and cosmic turbulence, Classical Quantum Gravity **21**, L27 (2004).

[18] Chiara Caprini and Ruth Durrer, Gravitational waves from stochastic relativistic sources: Primordial turbulence and magnetic fields, Phys. Rev. D **74**, 063521 (2006).
[19] Grigol Gogoberidze, Tina Kahniashvili, and Arthur Kosowsky, The spectrum of gravitational radiation from primordial turbulence, Phys. Rev. D **76**, 083002 (2007).

[20] Tigran Kalaydzhyan and Edward Shuryak, Gravity waves generated by sounds from big bang phase transitions, Phys. Rev. D **91**, 083502 (2015).

[21] Edward W Kolb and Michael Stanley Turner, *The Early Universe*, Frontiers in Physics (Westview Press, Boulder, CO, 1990).

[22] Sunny Vagnozzi and Abraham Loeb, The challenge of ruling out inflation via the primordial graviton background, Astrophys. J. Lett. **939**, L22 (2022).

[23] Chiara Caprini and Daniel G. Figueroa, Cosmological backgrounds of gravitational waves, *Classical Quantum Gravity* **35**, 163001 (2018).

[24] Andreas Ringwald and Carlos Tamarit, Revealing the cosmic history with gravitational waves, Phys. Rev. D **106**, 063027 (2022).

[25] Steven Weinberg, *Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity* (John Wiley and Sons, New York, 1972).

[26] J. Ghiglieri and M. Laine, Gravitational wave background from standard model physics: Qualitative features, J. Cosmol. Astropart. Phys. 07 (2015) 022.

[27] J. Ghiglieri, G. Jackson, M. Laine, and Y. Zhu, Gravitational wave background from standard model physics: Complete leading order, J. High Energy Phys. 07 (2020) 092.

[28] Andreas Ringwald, Jan Schütte-Engel, and Carlos Tamarit, Gravitational waves as a big bang thermometer, J. Cosmol. Astropart. Phys. 03 (2021) 054.

[29] Andrei Sakharov, Maximum temperature of thermal radiation, Pis'ma Zh. Eksp. Teor. Fiz. **3**, 439 (1966).

[30] Gary N. Felder, Lev Kofman, and Andrei D. Linde, Instant preheating, Phys. Rev. D **59**, 123523 (1999).

[31] Mustafa A. Amin, Mark P. Hertzberg, David I. Kaiser, and Johanna Karouby, Nonperturbative dynamics of reheating after inflation: A review, Int. J. Mod. Phys. D **24**, 1530003 (2014).

[32] Robert Brandenberger and Patrick Peter, Bouncing cosmologies: Progress and problems, Found. Phys. **47**, 797 (2017).

[33] Robert Brandenberger and Ziwei Wang, Nonsingular Ekyrotic cosmology with a nearly scale-invariant spectrum of cosmological perturbations and gravitational waves, Phys. Rev. D **101**, 063522 (2020).

[34] Betty X. Hu and Abraham Loeb, An upper limit on the initial temperature of the radiation-dominated Universe, J. Cosmol. Astropart. Phys. 01 (2021) 041.

[35] M. Kawasaki, Kazunori Kohri, and Naoshi Sugiyama, Cosmological constraints on late time entropy production, Phys. Rev. Lett. **82**, 4168 (1999).

[36] M. Kawasaki, Kazunori Kohri, and Naoshi Sugiyama, MeV scale reheating temperature and thermalization of neutrino background, Phys. Rev. D **62**, 023506 (2000).

[37] Gian Francesco Giudice, Edward W. Kolb, and Antonio Riotto, Largest temperature of the radiation era and its cosmological implications, Phys. Rev. D **64**, 023508 (2001).

[38] Steen Hannestad, What is the lowest possible reheating temperature?, Phys. Rev. D **70**, 043506 (2004).

[39] Takuya Hasegawa, Nagisa Hiroshima, Kazunori Kohri, Rasmus S. L. Hansen, Thomas Tram, and Steen Hannestad,

MeV-scale reheating temperature and thermalization of oscillating neutrinos by radiative and hadronic decays of massive particles, J. Cosmol. Astropart. Phys. 12 (2019) 012.

[40] M. Fukugita and T. Yanagida, Baryogenesis without grand unification, Phys. Lett. B **174**, 45 (1986).

[41] Lucía Castells-Tiestos and Jorge Casalderrey-Solana, Thermal emission of gravitational waves from weak to strong coupling, J. High Energy Phys. 10 (2022) 049.

[42] P. Klose, M. Laine, and S. Procacci, Gravitational wave background from non-Abelian reheating after axion-like inflation, J. Cosmol. Astropart. Phys. 05 (2022) 021.

[43] P. Klose, M. Laine, and S. Procacci, Gravitational wave background from vacuum and thermal fluctuations during axion-like inflation, J. Cosmol. Astropart. Phys. 12 (2022) 020.

[44] S. Y. Choi, J. S. Shim, and H. S. Song, Factorization and polarization in linearized gravity, Phys. Rev. D **51**, 2751 (1995).

[45] Mikko Laine and Aleksi Vuorinen, Basics of thermal field theory, Lect. Notes Phys. **925**, 1 (2016).

[46] T. Kinoshita, Mass singularities of Feynman amplitudes, J. Math. Phys. (N.Y.) **3**, 650 (1962).

[47] T. D. Lee and M. Nauenberg, Degenerate systems and mass singularities, Phys. Rev. **133**, B1549 (1964).

[48] M. Laine, Resonant *s*-channel dark matter annihilation at NLO, J. High Energy Phys. 01 (2023) 157.

[49] Nora Weickgenannt, Xin-Li Sheng, Enrico Speranza, Qun Wang, and Dirk H. Rischke, Kinetic theory for massive spin-1/2 particles from the Wigner-function formalism, Phys. Rev. D **100**, 056018 (2019).

[50] Nora Weickgenannt, Enrico Speranza, Xin-li Sheng, Qun Wang, and Dirk H. Rischke, Generating spin polarization from vorticity through nonlocal collisions, Phys. Rev. Lett. **127**, 052301 (2021).

[51] Di-Lun Yang, Koichi Hattori, and Yoshimasa Hidaka, Effective quantum kinetic theory for spin transport of fermions with collsional effects, J. High Energy Phys. 07 (2020) 070.

[52] Nora Weickgenannt, Enrico Speranza, Xin-li Sheng, Qun Wang, and Dirk H. Rischke, Derivation of the nonlocal collision term in the relativistic Boltzmann equation for massive spin-1/2 particles from quantum field theory, Phys. Rev. D **104**, 016022 (2021).

[53] Xin-Li Sheng, Nora Weickgenannt, Enrico Speranza, Dirk H. Rischke, and Qun Wang, From Kadanoff-Baym to Boltzmann equations for massive spin-1/2 fermions, Phys. Rev. D **104**, 016029 (2021).

[54] D. Bödeker, M. Sangel, and M. Wörmann, Equilibration, particle production, and self-energy, Phys. Rev. D **93**, 045028 (2016).

[55] Adam Alloul, Neil D. Christensen, Céline Degrande, Claude Duhr, and Benjamin Fuks, FeynRules 2.0 - A complete toolbox for tree-level phenomenology, Comput. Phys. Commun. **185**, 2250 (2014).

[56] Neil D. Christensen, Priscila de Aquino, Celine Degrande, Claude Duhr, Benjamin Fuks, Michel Herquet, Fabio Maltoni, and Steffen Schumann, A comprehensive approach to new physics simulations, Eur. Phys. J. C **71**, 1541 (2011).
[57] Thomas Hahn, Generating Feynman diagrams and amplitudes with FeynArts 3, Comput. Phys. Commun. **140**, 418 (2001).

[58] Thomas Hahn, Sebastian Paßehr, and Christian Schappacher, FormCalc 9 and extensions, Proc. Sci. LL2016 (**2016**) 068 [arXiv:1604.04611].

[59] Barry R. Holstein, Graviton physics, Am. J. Phys. **74**, 1002 (2006).

[60] N. E. J. Bjerrum-Bohr, Barry R. Holstein, Ludovic Planté, and Pierre Vanhove, Graviton-photon scattering, Phys. Rev. D **91**, 064008 (2015).

[61] M. Laine and M. Meyer, Standard Model thermodynamics across the electroweak crossover, J. Cosmol. Astropart. Phys. **07** (2015) 035.

[62] R. L. Workman et al. (Particle Data Group), Review of particle physics, Prog. Theor. Exp. Phys. **2022**, 083C01 (2022).

[63] Ken'ichi Saikawa and Satoshi Shirai, Primordial gravitational waves, precisely: The role of thermodynamics in the standard model, J. Cosmol. Astropart. Phys. **05** (2018) 035.

[64] Francesco Muia, Fernando Quevedo, Andreas Schachner, and Gonzalo Villa, Testing BSM physics with gravitational waves, J. Cosmol. Astropart. Phys. **09** (2023) 006.

[65] Asher Berlin, Diego Blas, Raffaele Tito D'Agnolo, Sebastian A. R. Ellis, Roni Harnik, Yonatan Kahn, and Jan Schütte-Engel, Detecting high-frequency gravitational waves with microwave cavities, Phys. Rev. D **105**, 116011 (2022).

[66] Nancy Aggarwal et al., Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies, Living Rev. Relativity **24**, 4 (2021).

[67] G. Peter Lepage, Adaptive multidimensional integration: VEGAS enhanced, J. Comput. Phys. **439**, 110386 (2021).

[68] G. Baym, H. Monien, C. J. Pethick, and D. G. Ravenhall, Transverse interactions and transport in relativistic quark-gluon and electromagnetic plasmas, Phys. Rev. Lett. **64**, 1867 (1990).

[69] Denis Besak and Dietrich Bodeker, Thermal production of ultrarelativistic right-handed neutrinos: Complete leading-order results, J. Cosmol. Astropart. Phys. **03** (2012) 029.