# Double-graviton production from Standard Model plasma

J Ghiglieri, M Laine, J Schütte-Engel, E Speranza

► To cite this version:

J Ghiglieri, M Laine, J Schütte-Engel, E Speranza. Double-graviton production from Standard Model plasma. Journal of Cosmology and Astroparticle Physics, 2024, 04, pp.062. 10.1088/1475-7516/2024/04/062. hal-04409355

HAL Id: hal-04409355

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

Submitted on 15 Jan 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
PAPER • OPEN ACCESS

Double-graviton production from Standard Model
plasma

To cite this article: J. Ghiglieri *et al* JCAP04(2024)062

View the [article online](#) for updates and enhancements.

You may also like

- [Inflationary models with non-minimally derivative coupling](#)
Nan Yang, Qin Fei, Qing Gao et al.

- [Scalar dark energy models and scalar-tensor gravity: theoretical explanations for the accelerated expansion of the present Universe](#)
Peixiang Ji and Lijing Shao

- [Graviton Momentum: A Natural Source of Dark Energy](#)
Samuel Meng and
Double-graviton production from Standard Model
plasma

J. Ghiglieri<sup>1, ∗</sup> , <sup>a</sup> M. Laine<sup>1, †</sup> , <sup>b</sup> J. Schütte-Engel<sup>1, c,d</sup> and E. Speranza<sup>1, e</sup>

E-mail: jacopo.ghiglieri@subatech.in2p3.fr, laine@itp.unibe.ch, janschue@berkeley.edu, enrico.speranza@cern.ch

ABSTRACT: The thermal plasma filling the early universe generated a stochastic gravitational wave background that peaks in the microwave frequency range today. If the graviton production rate is expressed as a series in a fine-structure constant, α, and the temperature over the Planck mass, $T^2/m_{\text{pl}}^2$, then the lowest-order contributions come from single ($\sim \alpha T^2/m_{\text{pl}}^2$) and double ($\sim T^4/m_{\text{pl}}^4$) graviton production via $2 \to 2$ scatterings. We show that in the Standard Model, single-graviton production dominates if the maximal temperature is smaller than $4 \times 10^{18}$ GeV. This justifies previous calculations which relied solely on single-graviton production. We mention Beyond the Standard Model scenarios in which the single and double-graviton contributions could be of comparable magnitudes. Finally, we elaborate on what these results imply for the range of applicability of General Relativity as an effective theory.

KEYWORDS: particle physics — cosmology connection, primordial gravitational waves (theory), inflation

ARXIV ePRINT: 2401.08766

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

<footnote><sup>a</sup>SUBATECH, Nantes Université, IMT Atlantique, IN2P3/CNRS,
4 rue Alfred Kastler, La Chantrerie BP 20722, 44307 Nantes, France</footnote>

<footnote><sup>b</sup>AEC, Institute for Theoretical Physics, University of Bern,
Sidlerstrasse 5, CH-3012 Bern, Switzerland</footnote>

<footnote><sup>c</sup>Department of Physics, University of California,
Berkeley, CA 94720, U.S.A.</footnote>

<footnote><sup>d</sup>RIKEN iTHEMS,
Wako, Saitama 351-0198, Japan</footnote>

<footnote><sup>e</sup>Theoretical Physics Department, CERN,
CH-1211 Geneva 23, Switzerland</footnote>
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Steps of the computation</td><td>2</td></tr><tr><td>3 Numerical results</td><td>6</td></tr><tr><td>4 Conclusions and outlook</td><td>8</td></tr><tr><td>A Azimuthal and radial momentum integrals</td><td>10</td></tr></table>

## 1 Introduction

The past years have seen an explosion of interest in gravitational wave astronomy. With hopes for an upcoming observation of a primordial background through pulsar timing arrays [1–4], the construction of the LISA observatory under way, and the planning of DECIGO as well as the Einstein Telescope striding forward, theoretical efforts at identifying all conceivable cosmological sources are well motivated. In fact such efforts should cover not only the sub-kHz frequencies relevant for large interferometers, but extend up to much higher frequencies, 100 GHz and above, for which future experimental concepts are being actively developed [5].

If we restrict ourselves to the Standard Model of particle physics, perhaps minimally extended to incorporate the existence of neutrino masses, there is however only one primordial gravitational wave component that is guaranteed to be present.<fnref n="1" /> It is the gravitational waves produced by microscopic collisions [11] taking place in a Standard Model plasma [12, 13]. This background is *not* a perfect analogue of the cosmic microwave background (CMB), since gravitons were most likely never in thermal equilibrium,<fnref n="2" /> however its shape (~ Planckian) and peak frequency range ($f_0|_{\text{peak}} \sim 100$ GHz) are similar to those of the CMB. The reason for the similarity is that the graviton energies reflect the thermal energy distributions of the plasma particles from which they are produced.

A key property of the thermal gravitational wave background is that while its shape can be computed, its amplitude is not known. The reason is that gravitons do not equilibrate; they are just continuously being produced. The production rate is proportional to the strength of the gravitational interaction, ~ $1/m_{\text{pl}}^2$, where $m_{\text{pl}}$ is the Planck mass, as well as to the strength of Standard Model forces, ~ $\alpha$, where $\alpha$ is a fine-structure constant. Integrating over a Hubble time, $H^{-1} \sim m_{\text{pl}}/T^2$, and normalizing to the energy density of radiation, the fractional energy density in gravitational waves then scales as $\Omega_{\text{gw}} \sim \alpha T_{\text{max}}/m_{\text{pl}}$, where $T_{\text{max}}$ is the maximal temperature of the radiation-dominated phase in our universe.

Now, if we explore large values of $T_{\text{max}}$, it must be asked whether terms suppressed by higher powers of $T_{\text{max}}/m_{\text{pl}}$ also play a role. A concrete way to probe this question was

<footnote n="1">¹Many other sources have been analyzed, notably inflation [6], the dynamics of reheating [7, 8], and a plethora of post-reheating phenomena [9, 10], however they rely on Beyond the Standard Model (BSM) physics.</footnote>

<footnote n="2">²Exotic scenarios in which this could have been the case have also been proposed, cf., e.g., refs. [14, 15].</footnote>
identified in ref. [16]. The rate ~ $\alpha/m_{pl}^2$ mentioned above involves the production of a single graviton in association with a Standard Model particle. But in addition, there are processes in which two gravitons are produced, without any Standard Model particles. The rate of the latter processes is ~ $1/m_{pl}^4$. To complete the dimensions, this rate must come with a power of $T_{max}$. Therefore, at high $T_{max}$, the double-graviton rate overtakes the single-graviton rate. Then we also need to ask whether the notion of using General Relativity as an effective theory, reliant on the convergence of an expansion in $T_{max}/m_{pl}^2$, continues to apply.

The purpose of the current study is to promote the computation of ref. [16], carried out in scalar field theory, to the full Standard Model. This requires including fermion-antifermion pairs and gauge boson pairs as additional initial states. Summing over all processes, we provide a quantitative criterion for when single-graviton production dominates.

Our presentation is organized as follows. The main steps of the computation, as well as a difference to matrix elements squared that can be found in the literature, are discussed in section 2. Our results are illustrated numerically in section 3, and the main conclusions are formulated in section 4. Some technical details related to evaluating thermal averages of pair production cross sections have been relegated to appendix A.

## 2 Steps of the computation

We consider a moment after inflation at which the plasma filling the universe has equilibrated to an average temperature $T \ll m_{pl}$. Furthermore, the plasma contribution to the overall energy density is assumed to dominate over that from non-equilibrated fields, such as the inflaton. Then the Hubble rate, $H$, is much smaller than typical thermal momenta, $H \sim \sqrt{g_*} T^2/m_{pl} \ll \pi T$, where $g_* \sim 10^2$ is the number of relativistic degrees of freedom. Consequently, thermal wavelengths are well within the horizon, $(\pi T)^{-1} \ll H^{-1}$. We consider the production of gravitational waves (GW's) from the scatterings of particles with such wavelengths. As will become apparent, the GW spectrum reflects that of the scatterers, so that in terms of the physical momentum $k$, it is peaked at around $k|_{peak} \sim \pi T$.

Given that the physics we are interested in takes place within the horizon, the space is locally flat, and we can employ local Minkowskian coordinates for the computation. The background metric is written as $g_{\mu\nu} = \eta_{\mu\nu} + \kappa h_{\mu\nu}$, where $\eta_{\mu\nu} = \text{diag}(+---)$ is the Minkowski metric and $h_{\mu\nu}$ denotes the gravitational perturbation. The effective coupling of gravitational interactions is defined as

$$ \kappa^2 \equiv 32\pi G \equiv \frac{32\pi}{m_{pl}^2}, \qquad (2.1) $$

where $G$ is Newton's constant, and $m_{pl} = 1.22091 \times 10^{19}$ GeV is the Planck mass.

When the Einstein-Hilbert action is expanded as a series in $\kappa$, vertices are generated for the interactions of $h_{\mu\nu}$ with itself and with other fields (a pedagogic review can be found in ref. [17]). Adding a graviton leg to a process leads to a suppression of the amplitude by $\kappa \sim 1/m_{pl}$, and to a suppression of the corresponding amplitude squared by $\kappa^2$. But we can also add a normal Standard Model vertex to any given process, which leads to the suppression of the rate by ~ $\alpha$, where $\alpha$ is a fine-structure constant. Therefore, the production rate can be viewed as a double series, in $\alpha$ and $1/m_{pl}^2$. In addition, we may recall that $2 \to 1$ processes
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Diagrams leading to the pair production of gravitons (double wiggly lines) from (a) scalars (dashed lines); (b) fermions (solid lines); (c) gauge fields (wavy lines). We have omitted arrows from scalars and fermions, understanding that both particle and antiparticle states are to be included.</figcaption></figure>

are kinematically forbidden for a massless final state. If we approximate all particles as massless, the leading processes are then $2 \to 2$ amplitudes. If they involve a Standard Model vertex, their rate is $\sim \alpha/m_{\text{pl}}^2$; if they only contain gravitational vertices, their rate is $\sim 1/m_{\text{pl}}^4$. To these leading terms, additional vertices of either type can be added, leading to further suppression by $\alpha$ or $1/m_{\text{pl}}^2$ (e.g. triple-graviton production at $\sim 1/m_{\text{pl}}^6$, etc).

The practical computations can be streamlined by going over to the De Donder gauge, $\partial_\alpha h_\mu^\alpha = \frac{1}{2}\partial_\mu h_\alpha^\alpha$. This also renders the quadratic part of the gravitational action invertible, so that an internal graviton propagator can be found. On the external graviton lines, which are on-shell, i.e. with four-momentum $K = (k, k)$ with $k \equiv |k|$, the sum over polarizations $\lambda$ yields the projector

$$ \sum_{\lambda=\{\times,+\}} \epsilon_{\alpha\beta}^{\lambda}(\mathbf{k}) \epsilon_{\mu\nu}^{*\lambda}(\mathbf{k}) \equiv \mathbb{L}_{\alpha\beta;\mu\nu} = \frac{1}{2} (\mathbb{K}_{\alpha\mu}^{\text{T}} \mathbb{K}_{\beta\nu}^{\text{T}} + \mathbb{K}_{\alpha\nu}^{\text{T}} \mathbb{K}_{\beta\mu}^{\text{T}} - \mathbb{K}_{\alpha\beta}^{\text{T}} \mathbb{K}_{\mu\nu}^{\text{T}}), \quad (2.2) $$

which is traceless in $\alpha\beta$ and $\mu\nu$, transverse with respect to the graviton four-momentum $K$, and projects onto two physical states in four dimensions, i.e. $\mathbb{L}_{\alpha\beta;\alpha\beta} = 2$. The tensor $\mathbb{K}^\text{T}$ can be chosen as

$$ \mathbb{K}_{\mu\nu}^{\text{T}} \equiv \eta_{\mu i} \eta_{\nu j} \left( \delta^{ij} - \frac{k^i k^j}{k^2} \right) = -\eta_{\mu\nu} + \frac{\mathcal{K}_\mu \bar{\mathcal{K}}_\nu + \mathcal{K}_\nu \bar{\mathcal{K}}_\mu}{\mathcal{K} \cdot \bar{\mathcal{K}}} , \quad (2.3) $$

which is transverse with respect to the four-momenta $K$ and $\bar{K} \equiv (k, -k)$, and also with respect to the three-momentum $k$ (we work in the plasma rest frame).

With the Feynman rules at hand, the first step is to determine the matrix elements squared for the production processes, shown in figure 1. The four-momenta of the Standard Model particles are denoted by $P_1, P_2$ and those of the gravitons by $K_1, K_2$, with $P_1 + P_2 = K_1 + K_2$.
We introduce the usual Mandelstam invariants, $s \equiv (P_1+P_2)^2$, $t \equiv (P_1-K_1)^2$, $u \equiv (P_1-K_2)^2$. The processes take place at very high temperatures, much above any known particles masses. Therefore all particles can effectively be taken to be massless, whereby

$$s+t+u=0. \qquad (2.4)$$

For massless fermions, sums over spins yield the standard expressions $\sum_{\tau=\pm} u(\tau, \mathbf{p})\bar{u}(\tau, \mathbf{p}) = \sum_{\tau=\pm} v(\tau, \mathbf{p})\bar{v}(\tau, \mathbf{p}) = \mathcal{P}$. For massless gauge bosons, a polarization sum yields

$$\sum_{\tau=\pm} \epsilon_{\mu}^{\tau}(\mathbf{p}) \epsilon_{\nu}^{\tau*}(\mathbf{p}) \equiv \mathcal{P}_{\mu\nu}^{T}, \qquad (2.5)$$

where $\mathcal{P}_{\mu\nu}^{T}$ is defined like in eq. (2.3), but with four-momentum $\mathcal{P}$ rather than $K$. It projects onto two physical states in four dimensions, i.e. $|\mathcal{P}_{\mu}^{T}| = 2$.

Inserting the Feynman rules to the amplitudes in figure 1, with the cubic graviton vertex being particularly cumbersome [17]; constructing the amplitudes squared; and contracting with polarization or spin sums, as specified above, leads to rather lengthy expressions. We have found it practical to manipulate them with FORM [18], FEYNARTS [19], FORMCALC [20], and/or FEYNCALC [21]. In the course of these computations, we have verified the gauge independence of the contribution originating from massless gauge bosons, notably that the last term of eq. (2.3) (with $K \to \mathcal{P}$), containing $\bar{\mathcal{P}}$, does not contribute.

The matrix elements squared are best tabulated in a form where we sum not only over all spins and polarization states, both of initial and final-state particles, but also over particles and antiparticles. For a complex scalar field $\phi$, the processes from figure 1(a) then yield [16]

$$\sum_{\text{all}} |\mathcal{M}_{\phi\phi^*\to hh}|^2 + \sum_{\text{all}} |\mathcal{M}_{\phi^*\phi\to hh}|^2 = \frac{\kappa^4 t^2 u^2}{4 s^2} = \frac{\kappa^4}{4} \Phi_{\text{scalar}}. \qquad (2.6)$$

In the fermionic case, we take a Dirac fermion of definite chirality (i.e. a Weyl fermion) as a building block, with two physical degrees of freedom. However, as we sum over fermions and antifermions, we effectively have four degrees of freedom for each fermionic external leg, and therefore each of them can effectively be represented as a usual Dirac fermion $\psi$. The processes in figure 1(b) then lead to

$$\sum_{\text{all}} |\mathcal{M}_{\psi\bar{\psi}\to hh}|^2 = \frac{\kappa^4}{4} \left( -\frac{t^2 u^2}{s^2} + \frac{t u}{2} \right) = \frac{\kappa^4}{4} \Phi_{\text{fermion}}. \qquad (2.7)$$

Finally, for a U(1) gauge boson $g$, the processes in figure 1(c) produce

$$\sum_{\text{all}} |\mathcal{M}_{gg\to hh}|^2 = \frac{\kappa^4}{4} \left( \frac{t^2 u^2}{s^2} - 2t u + \frac{s^2}{2} \right) = \frac{\kappa^4}{4} \Phi_{\text{gauge}}. \qquad (2.8)$$

We note that, up to overall normalization, eqs. (2.6) and (2.7) can be extracted as the massless limits of differential cross sections given in ref. [22]. However, this is not the case for eq. (2.8), as the assignment of a mass to vector fields is ambiguous and breaks gauge invariance. If masses are given to gauge fields, then the physical origin of the longitudinal polarization states should be specified (e.g. through the Higgs mechanism), in order to obtain a gauge independent expression.
Given the matrix elements squared, we can determine the thermally averaged production rate of gravitational radiation. As a tool for this, we introduce the polarization-averaged phase-space distribution of gravitons, $f_{gw}$. The energy density carried by gravitational waves of physical momentum $k$ can now be expressed as

$$de_{gw} = 2k f_{gw} \frac{d^3 k}{(2\pi)^3}, \qquad (2.9)$$

where the factor 2 counts the physical polarization states. In a locally Minkowskian time, making use of the rotational invariance of the thermal ensemble, the production rate can subsequently be expressed as

$$\frac{de_{gw}}{dt d\ln k} = \frac{k^4 \dot{f}_{gw}}{\pi^2}. \qquad (2.10)$$

The time variation of the phase-space distribution can be extracted from a Boltzmann equation. Assuming the particle content of the Standard Model (with three generations of charged leptons $\ell$, neutrinos $\nu$, and up and down-type quarks $u, d$, each with left and right chiralities L and R), this takes at $O(\kappa^4)$ the form

$$\dot{f}_{gw} \stackrel{O(\kappa^4)}{\supset} \frac{1}{8k} \int d\Omega_{2\to 2} \frac{\kappa^4}{4} \left\{ f_B(p_1) f_B(p_2) \left[ \underbrace{2}_{SU_L(2)} \Phi_{\text{scalar}} \right] \right. \\ + f_F(p_1) f_F(p_2) \left[ 3 \left( \underbrace{3}_{\ell_L, \ell_R, \nu_L} + \underbrace{1}_{\nu_R} + \underbrace{4N_c}_{u_L, u_R, d_L, d_R} \right) \Phi_{\text{fermion}} \right] \\ \left. + f_B(p_1) f_B(p_2) \left[ \underbrace{1}_{U_Y(1)} + \underbrace{3}_{SU_L(2)} + \underbrace{N_c^2 - 1}_{SU_c(N_c)} \right] \Phi_{\text{gauge}} \right\}, \qquad (2.11)$$

where $p_i \equiv |\mathbf{p}_i|$, $N_c = 3$, and the Bose and Fermi distributions have been introduced as $f_B(\epsilon) \equiv 1/(e^{\epsilon/T} - 1)$ and $f_F(\epsilon) \equiv 1/(e^{\epsilon/T} + 1)$, respectively. The prefactor $1/(8k)$ is a combination of the standard $1/(2k)$ associated with phase space, a factor $1/2$ for cancelling the sum over graviton polarizations that we had included in $\sum |M|^2$, as well as $1/2$ for cancelling the overcounting of identical particles or the redundant sum over particles and antiparticles introduced above. Furthermore, we note that eq. (2.11) includes only the gain terms, no loss terms, because the graviton phase-space density is much below the equilibrium value, $f_{gw}(k) \ll f_B(k)$. For the same reason, there are no Bose enhancement factors for the final-state gravitons in the gain terms.

Let us remark that right-handed neutrinos have been included as degrees of freedom in eq. (2.11). However, they only interact via Yukawa couplings, whose magnitudes are unknown. If they are small, right-handed neutrinos might not equilibrate fast enough to be part of the thermal ensemble, and should be omitted as initial states. We will illustrate the numerical difference between including and not including right-handed neutrinos in figure 2. It is perhaps appropriate to mention that the equilibration rates of many other particles have also been discussed in the literature, however this typically concerns particle asymmetry changing reactions relevant for chemical equilibration (cf., e.g., ref. [23]), whereas kinetic equilibrium should be efficiently maintained by gauge interactions.
The phase-space average in eq. (2.11) can be conveniently implemented by taking the s-channel momentum transfer as the integration variable. For this, we write $P_i = (p_i, \mathbf{p}_i)$ and $K_i = (k_i, \mathbf{k}_i)$ with $K \equiv K_2$. Then,

$$ \begin{aligned} \int d\Omega_{2\leftrightarrow2} &\equiv \int \frac{d^3\mathbf{p}_1 d^3\mathbf{p}_2 d^3\mathbf{k}_1 d^4\mathcal{Q}}{(2p_1)(2p_2)(2k_1)(2\pi)^9} (2\pi)^4 \delta^{(4)}(\mathcal{P}_1 + \mathcal{P}_2 - \mathcal{Q}) \delta^{(4)}(\mathcal{Q} - \mathcal{K}_1 - \mathcal{K}) \\ &= \int \frac{d^3\mathbf{p}_2 dq_0 d^3\mathbf{q}}{|\mathbf{q} - \mathbf{p}_2| p_2 |\mathbf{q} - \mathbf{k}|} \frac{\delta(|\mathbf{q} - \mathbf{p}_2| + p_2 - q_0) \delta(q_0 - |\mathbf{q} - \mathbf{k}| - k)}{8(2\pi)^5}, \end{aligned} \quad (2.12) $$

where integrals were carried out over $\mathbf{p}_1$ and $\mathbf{k}_1$. The constraints fix two angles as

$$ \cos \theta_{\mathbf{q},\mathbf{p}_2} = \frac{q^2 + q_0(2p_2 - q_0)}{2qp_2}, \quad \cos \theta_{\mathbf{q},\mathbf{k}} = \frac{q^2 + q_0(2k - q_0)}{2qk}. \quad (2.13) $$

The remaining angle can be expressed in terms of an azimuthal variable, $\varphi$, as

$$ \cos \theta_{\mathbf{k},\mathbf{p}_2} = \cos \theta_{\mathbf{q},\mathbf{k}} \cos \theta_{\mathbf{q},\mathbf{p}_2} + \sin \theta_{\mathbf{q},\mathbf{k}} \sin \theta_{\mathbf{q},\mathbf{p}_2} \cos \varphi. \quad (2.14) $$

Carrying out the angular integrals in order to remove the Dirac-$\delta$'s, we then obtain

$$ \int d\Omega_{2\leftrightarrow2} = \frac{1}{(4\pi)^3 k} \int_k^\infty dq_0 \int_{|2k-q_0|}^{q_0} dq \int_{q_-}^{q_+} dp_2 \int_{-\pi}^\pi \frac{d\varphi}{2\pi}, \quad q_\pm \equiv \frac{q_0 \pm q}{2}. \quad (2.15) $$

Subsequently, the integrals over $p_2$ and $\varphi$ can be performed analytically, as detailed in appendix A. We represent the result as

$$ \int_{q_{-}}^{q_{+}} dp_{2} \int_{-\pi}^{\pi} \frac{d\varphi}{2\pi} f_{\text{B}}(q_{0}-p_{2}) f_{\text{B}}(p_{2}) \Phi_{\text{scalar}} \equiv f_{\text{B}}(q_{0}) (q_{0}^{2}-q^{2})^{2} \Theta_{\text{scalar}}, \quad (2.16) $$

and correspondingly for $\Theta_{\text{fermion}}$ and $\Theta_{\text{gauge}}$. The expressions originating from eqs. (2.6), (2.7) and (2.8) are given in eqs. (A.7), (A.9) and (A.11), respectively.

The physical production rate can now be assembled from eqs. (2.10), (2.11), (2.15), and (2.16). It entails a 2-dimensional integral,

$$ \frac{\mathrm{d}e_{\mathrm{gw}}}{\mathrm{d}t \mathrm{d} \ln k} \supset \frac{\kappa^4 k^2}{2(4\pi)^5} \int_k^\infty dq_0 \int_{|2k-q_0|}^{q_0} dq f_B(q_0) (q_0^2 - q^2)^2 \\ \times \{2 \Theta_{\mathrm{scalar}} + 12 (1 + N_c) \Theta_{\mathrm{fermion}} + (3 + N_c^2) \Theta_{\mathrm{gauge}}\}, \quad (2.17) $$

which is rapidly convergent and readily evaluated numerically.<fnref n="3" />

## 3 Numerical results

A numerical illustration of eq. (2.17) is given in figure 2. We have plotted the rate in units of $T^7/m_{\text{pl}}^2$, whereby it scales with the temperature as $T^2/m_{\text{pl}}^2$. The reason for this choice is that the result can then be conveniently compared with the rate originating from single-graviton

<footnote n="3">³In this respect the rate at $O(\kappa^4)$ differs from that at $O(\alpha\kappa^2)$: the latter is logarithmically IR divergent in naive perturbation theory, and requires Hard Thermal Loop resummation [24, 25], in order to obtain a finite result [13]. After the resummation, it contains a logarithmically enhanced term, of $O(\alpha\kappa^2 \ln(1/\alpha))$ [12].</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Left: the single (h) [13] and double (hh) graviton production rates from a Standard Model plasma. The slow evolution of the single-graviton rate with the temperature is due to a logarithmic running of the couplings, whereas the double-graviton rate scales with the temperature as $T^2/m_{pl}^2$ in the units chosen. At $T = 10^{18}$ GeV, the double-graviton rate is still more than an order of magnitude smaller than the single-graviton rate, but it would become dominant at $T > 4 \times 10^{18}$ GeV. Right: different contributions to the double-graviton rate from eq. (2.17) at $T = 10^{18}$ GeV. In the fermionic part and the sum, the narrow band indicates the effect of omitting right-handed neutrinos.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The time-integrated contribution of single and double-graviton production to the differential gravitational energy density today, $h^2d\Omega_{gw,0}/d\ln f_0$, as a function of the current-day frequency, $f_0$, and various maximal temperatures, $T_{max}$. We have assumed radiation-dominated expansion at $T < T_{max}$ [33], and denoted by $h$ the reduced Hubble rate. The results are compared with a curve of the same shape (grey line), which after integration over $f_0$ yields $h^2\Omega_{gw,0} \equiv 5.62 \times 10^{-6}\Delta N_{eff}$ (cf., e.g., ref. [9]). Here we have set $\Delta N_{eff} = 0.2$, which is the maximal observationally allowed contribution to the effective number of massless neutrino species [34].</figcaption></figure>
production [13], which is practically constant in these units. For both results we have assumed the Standard Model matter content. Additional particle species could be included in a fairly straightforward manner (cf., e.g., refs. [26–31]), but we do not expect them to change the overall pattern, even though they could have a quantitative influence.

The main observation from figure 2 is that in the domain of the peak frequency, the double-graviton production rate is below the single-graviton production rate at $T < 4 \times 10^{18}$ GeV. In order to put this in context, we recall that the inflationary constraint $H_{\text{max}} < 10^{-5}m_{\text{pl}}$ [32], combined with the assumption of instantaneous reheating, suggests $T_{\text{max}} < 10^{16}$ GeV. Therefore, within standard inflationary scenarios, the thermal part of the gravitational wave production rate can always be estimated from the single-graviton contribution.

If we ignore the standard inflationary paradigm, another constraint on gravitational waves originates from $N_{\text{eff}}$, parametrizing the energy density residing in light degrees of freedom at the time of primordial nucleosynthesis ($\Delta e_{\text{gw,bbn}} \equiv \Delta N_{\text{eff}} (7/8)(4/11)^{4/3}e_{\gamma,\text{bbn}}$). In order to obtain the gravitational contribution, we integrate the evolution of the phase space distribution, appearing in eq. (2.10), over the cosmological history, along a comoving trajectory [12]. The time integral can be converted into an integral over temperature, from the current $T_0$ up to a maximal one, $T_{\text{max}}$, yielding

$$ \frac{\mathrm{d}e_{\mathrm{gw},0}}{\mathrm{d}\ln k(t_0)} = - \int_{T_0}^{T_{\text{max}}} \mathrm{d}T \frac{\mathrm{d}t}{\mathrm{d}T} \frac{a^4(T)}{a^4(T_0)} \left[ \frac{\mathrm{d}e_{\mathrm{gw}}}{\mathrm{d}t \mathrm{d}\ln k} \right]_{k=k(t_0)} \frac{a(T_0)}{a(T)} \quad (3.1) $$

where $t_0$ is the current time, $a$ is the cosmological scale factor, and the relationship $t \leftrightarrow T$ is to be determined from the Friedmann equations. Subsequently, normalizing with the current critical energy density, $e_{\text{crit}} = h^2 e_{\gamma,0}/2.473 \times 10^{-5}/(T_0/2.7255 \text{ K})^4$, where $e_{\gamma,0}$ is the current photon energy density and $h$ is the reduced Hubble rate, and converting $k(t_0)$ into the current frequency, $f_0$, yields the gravitational energy fraction, $h^2 d\Omega_{\text{gw},0}/d\ln f_0$.<fnref n="4" /> The result is shown in figure 3. Once we integrate over $f_0$, obtaining $h^2 \Omega_{\text{gw},0}$, the total gravitational energy density contributes to $N_{\text{eff}}$ as $\Delta N_{\text{eff}} \approx h^2 \Omega_{\text{gw},0}/5.62 \times 10^{-6}$ (cf., e.g., ref. [9]).

Given that our result for $h^2 \Omega_{\text{gw},0}$ depends on $T_{\text{max}}$, we now get a relationship between $T_{\text{max}}$ and $\Delta N_{\text{eff}}$. This can be compared with existing theoretical and experimental knowledge about $\Delta N_{\text{eff}}$. In particular, the current estimated uncertainty of the Standard Model prediction, $\Delta N_{\text{eff}} \le 0.001$, corresponds to $T_{\text{max}} \le 2 \times 10^{17}$ GeV [13]. At $T_{\text{max}} = 4 \times 10^{18}$ GeV, when the double-graviton rate overtakes the single-graviton rate, $\Delta N_{\text{eff}} \approx 0.02$.<fnref n="5" /> This is still below the current experimental limit, $\Delta N_{\text{eff}} < 0.2$ [34]. For illustration, a spectrum that would yield $\Delta N_{\text{eff}} = 0.2$ is also shown in figure 3.

# 4 Conclusions and outlook

In this work we have computed the double-graviton production rate from a thermal Standard Model plasma present during an early radiation-dominated epoch. Previous calculations

<footnote n="4">⁴In the literature, $\Omega_{\text{gw},0}$ sometimes denotes the differential spectrum, but here we need a separate notation for the differential $(d\Omega_{\text{gw},0}/d\ln f_0)$ and integrated spectra $(\Omega_{\text{gw},0})$.</footnote>

<footnote n="5">⁵As mentioned in section 1, the single-graviton production rate involves $\sim \alpha/m_{\text{pl}}^2$, but the age of the universe scales with $\sim m_{\text{pl}}$, whereby the integrated contribution to $N_{\text{eff}}$ behaves as $\sim \alpha T_{\text{max}}/m_{\text{pl}}$.</footnote>
had only considered scalar field contributions to the double-graviton rate [16]. For momenta $k \sim \pi T$, the double-graviton rate scales as $de_{gw}/[dt \, d\ln k] \sim T^9/m_{pl}^4$. This can be compared with the single-graviton production rate [13], which is of order $de_{gw}/[dt \, d\ln k] \sim \alpha T^7/m_{pl}^2$, where $\alpha$ is a fine-structure constant.

Our result for the double-graviton production rate is given in eq. (2.17) and illustrated numerically in figure 2. A comparison with the single-graviton production rate shows that the $O(T^9/m_{pl}^4)$ contribution overtakes the $O(\alpha T^7/m_{pl}^2)$ one at $T \sim 4 \times 10^{18}$ GeV. This temperature is much higher than expected to be reached after standard inflationary scenarios. Therefore, single-graviton production is normally the dominant source. If one considers BSM scenarios associated with a non-renormalizable coupling, $\sim 1/\Lambda^2$ (cf., e.g., ref. [27]), then we should effectively replace $\alpha \to T^2/\Lambda^2$ in the single-graviton rate. In these cases, if $\Lambda \sim m_{pl}$, the single and double-graviton contributions could be of comparable magnitudes.

In order to obtain our results, we needed to evaluate the matrix elements squared corresponding to the Feynman diagrams in figure 1, producing the results given in eqs. (2.6)–(2.8). Similar matrix elements squared have been considered in previous literature, either out of theoretical interest (cf., e.g., ref. [22]), or in view of the scattering of gravitational waves on matter (cf., e.g., ref. [35]). As we are concerned with very high temperatures, the Standard Model particles can be considered as massless for the purposes of our analysis. We noted that the literature results had assigned masses to gauge fields in a way which did not permit to recover the massless limit of eq. (2.8).

Let us elaborate briefly on the more general conceptual implications of our results. As explained below eq. (2.1), the gravitational wave production rate within General Relativity is a double series, in $\alpha$ and $T^2/m_{pl}^2$. Normally, to estimate the convergence of the series in $T^2/m_{pl}^2$, we would compare terms with a same power of $\alpha$. However, the term of $O(\alpha^0 T^2/m_{pl}^2)$ is absent, because it is kinematically forbidden, and the term of $O(\alpha^1 T^4/m_{pl}^4)$ is unknown, as it is of next-to-leading order, either in $\alpha$ with respect to the known contribution of $O(\alpha^0 T^4/m_{pl}^4)$, or in $T^2/m_{pl}^2$ with respect to the known contribution of $O(\alpha^1 T^2/m_{pl}^2)$. But we can still compare the known terms. Because of the absence of $\alpha$, the term of $O(\alpha^0 T^4/m_{pl}^4)$ overtakes the $O(\alpha^1 T^2/m_{pl}^2)$ contribution sooner than the proper probe of $O(\alpha^1 T^4/m_{pl}^4)$ would. Therefore, the criterion $T < 4 \times 10^{18}$ GeV can be seen as a conservative (or sufficient) condition for the highest temperature at which General Relativity together with the Standard Model can be employed as a self-consistent effective theory.

Returning to the phenomenological side, we end by remarking that reaching a high temperature after inflation does not guarantee a substantial thermal gravitational wave background. The reason is that if the inflaton equilibration rate is small compared with the Hubble rate, as may be expected for an extremely weakly coupled field, then inflation is followed by a period of matter domination (cf., e.g., ref. [36]). Similarly to a phase transition produced signal (cf., e.g., ref. [37]), the thermal background would be diluted during such a stage.

## Acknowledgments

J.G. was supported by the Agence Nationale de la Recherche (France), under grant ANR-22-CE31-0018 (AUTOTHERM), J.S.-E. by the National Science Foundation under cooperative
agreement 202027 and by the RIKEN-Berkeley Center, and E.S. by the European Union's Horizon Europe Research and Innovation Program, under the Marie Skłodowska-Curie grant agreement no. 101109747.

## A Azimuthal and radial momentum integrals

We show here how two of the integrals in eq. (2.15) can be carried out analytically, in order to determine the factors $\Theta$ defined in eq. (2.16) that enter our final result in eq. (2.17).

Making use of the identities

$$f_B(\epsilon_2)f_B(q_0 - \epsilon_2) = f_B(q_0)[1 + f_B(\epsilon_2) + f_B(q_0 - \epsilon_2)], \quad (A.1)$$

$$f_F(\epsilon_2)f_F(q_0 - \epsilon_2) = f_B(q_0)[1 - f_F(\epsilon_2) - f_F(q_0 - \epsilon_2)], \quad (A.2)$$

the definition in eq. (2.16) as well as the matrix element squared in eq. (2.6) imply that

$$\Theta_{\text{scalar}} = \int_{q_{-}}^{q_{+}} dp_{2} \int_{-\pi}^{\pi} \frac{d\varphi}{2\pi} \left( \frac{t^2 u^2}{s^4} \right) [1 + f_{B}(p_{2}) + f_{B}(q_{0} - p_{2})]. \quad (A.3)$$

Here $u = -s - t$, so that the integrand is a polynomial in $t/s$. Inserting four-momenta, we get $t/s = 2(\mathbf{k} \cdot \mathbf{p}_2 - k p_2)/s$, which after the use of eqs. (2.13) and (2.14) can be expressed as

$$\frac{t}{s} = \frac{(2k - q_0)(2p_2 - q_0) - q^2 + \cos \varphi \sqrt{q^2 - (2k - q_0)^2} \sqrt{q^2 - (2p_2 - q_0)^2}}{2q^2}. \quad (A.4)$$

The azimuthal average in eq. (A.3) reduces to averages of $\cos^n\varphi$, trivially carried out. Left over are integrals over powers of $p_2$, weighted by a Bose or Fermi distribution. These can be expressed in terms of polylogarithms, which we define as

$$l_{1b}(\epsilon) \equiv \ln(1 - e^{-\epsilon/T}), \qquad l_{1f}(\epsilon) \equiv \ln(1 + e^{-\epsilon/T}), \quad (A.5)$$

$$l_{\text{nb}}(\epsilon) \equiv \text{Li}_n(e^{-\epsilon/T}), \qquad l_{\text{nf}}(\epsilon) \equiv \text{Li}_n(-e^{-\epsilon/T}), \quad n \ge 2. \quad (A.6)$$

Thereby

$$\begin{aligned} \Theta_{\text{scalar}} = & \frac{q}{30} + \frac{T[q^2 - (2k - q^0)^2]^2}{8q^4} [l_{1b}(q_+) - l_{1b}(q_-)] \\ & - \frac{T^2[q^2 - (2k - q^0)^2][q^2 - 5(2k - q^0)^2]}{2q^5} [l_{2b}(q_+) + l_{2b}(q_-)] \\ & - \frac{T^3[5q^4 - 42q^2(2k - q^0)^2 + 45(2k - q^0)^4]}{2q^6} [l_{3b}(q_+) - l_{3b}(q_-)] \\ & - \frac{3T^4[3q^4 - 30q^2(2k - q^0)^2 + 35(2k - q^0)^4]}{q^7} [l_{4b}(q_+) + l_{4b}(q_-)] \\ & - \frac{6T^5[3q^4 - 30q^2(2k - q^0)^2 + 35(2k - q^0)^4]}{q^8} [l_{5b}(q_+) - l_{5b}(q_-)], \quad (A.7) \end{aligned}$$
$$ \Theta_{\text{fermion}} = \int_{q_{-}}^{q_{+}} dp_{2} \int_{-\pi}^{\pi} \frac{d\varphi}{2\pi} \left( -\frac{t^2 u^2}{s^4} + \frac{tu}{2s^2} \right) [1 - f_{\text{F}}(p_2) - f_{\text{F}}(q_0 - p_2)] \quad (\text{A.8}) $$

$$ \begin{aligned} = & \frac{q}{20} + \frac{T[q^4 - (2k - q^0)^4]}{8q^4} [l_{1f}(q_+) - l_{1f}(q_-)] \\ & + \frac{T^2[-3q^2 + 5(2k - q^0)^2](2k - q^0)^2}{2q^5} [l_{2f}(q_+) + l_{2f}(q_-)] \\ & + \frac{3T^3[q^4 - 12q^2(2k - q^0)^2 + 15(2k - q^0)^4]}{2q^6} [l_{3f}(q_+) - l_{3f}(q_-)] \\ & + \frac{3T^4[3q^4 - 30q^2(2k - q^0)^2 + 35(2k - q^0)^4]}{q^7} [l_{4f}(q_+) + l_{4f}(q_-)] \\ & + \frac{6T^5[3q^4 - 30q^2(2k - q^0)^2 + 35(2k - q^0)^4]}{q^8} [l_{5f}(q_+) - l_{5f}(q_-)], \end{aligned} \quad (\text{A.9}) $$

$$ \Theta_{\text{gauge}} = \int_{q_{-}}^{q_{+}} dp_{2} \int_{-\pi}^{\pi} \frac{d\varphi}{2\pi} \left( \frac{t^2 u^2}{s^4} - \frac{2tu}{s^2} + \frac{1}{2} \right) [1 + f_{\text{B}}(p_2) + f_{\text{B}}(q_0 - p_2)] \quad (\text{A.10}) $$

$$ \begin{aligned} = & \frac{q}{5} + \frac{T[q^4 + 6q^2(2k - q^0)^2 + (2k - q^0)^4]}{8q^4} [l_{1b}(q_+) - l_{1b}(q_-)] \\ & + \frac{T^2[3q^4 - 6q^2(2k - q^0)^2 - 5(2k - q^0)^4]}{2q^5} [l_{2b}(q_+) + l_{2b}(q_-)] \\ & + \frac{3T^3[q^4 + 6q^2(2k - q^0)^2 - 15(2k - q^0)^4]}{2q^6} [l_{3b}(q_+) - l_{3b}(q_-)] \\ & - \frac{3T^4[3q^4 - 30q^2(2k - q^0)^2 + 35(2k - q^0)^4]}{q^7} [l_{4b}(q_+) + l_{4b}(q_-)] \\ & - \frac{6T^5[3q^4 - 30q^2(2k - q^0)^2 + 35(2k - q^0)^4]}{q^8} [l_{5b}(q_+) - l_{5b}(q_-)]. \end{aligned} \quad (\text{A.11}) $$

## References

[1] NANOGRAV collaboration, *The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background*, Astrophys. J. Lett. 951 (2023) L8 [arXiv:2306.16213] [INSPIRE].

[2] EPTA and INPTA collaborations, *The second data release from the European Pulsar Timing Array — III. Search for gravitational wave signals*, Astron. Astrophys. 678 (2023) A50 [arXiv:2306.16214] [INSPIRE].

[3] D.J. Reardon et al., *Search for an isotropic gravitational-wave background with the Parkes Pulsar Timing Array*, Astrophys. J. Lett. 951 (2023) L6 [arXiv:2306.16215] [INSPIRE].

[4] H. Xu et al., *Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I*, Res. Astron. Astrophys. 23 (2023) 075024 [arXiv:2306.16216] [INSPIRE].

[5] N. Aggarwal et al., *Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies*, Living Rev. Rel. 24 (2021) 4 [arXiv:2011.12414] [INSPIRE].

[6] A.A. Starobinsky, *Spectrum of relict gravitational radiation and the early state of the universe*, JETP Lett. 30 (1979) 682 [INSPIRE].

[7] S.Y. Khlebnikov and I.I. Tkachev, *Relic gravitational waves produced after preheating*, Phys. Rev. D 56 (1997) 653 [hep-ph/9701423] [INSPIRE].
[8] K.D. Lozanov, *Lectures on reheating after inflation*, arXiv:1907.04402 [inSPIRE].

[9] C. Caprini and D.G. Figueroa, *Cosmological backgrounds of gravitational waves*, Class. Quant. Grav. **35** (2018) 163001 [arXiv:1801.04268] [inSPIRE].

[10] R. Roshan and G. White, *Using gravitational waves to see the first second of the universe*, arXiv:2401.04388 [inSPIRE].

[11] S. Weinberg, *Gravitation and Cosmology*, John Wiley & Sons, New York, NY, U.S.A. (1971).

[12] J. Ghiglieri and M. Laine, *Gravitational wave background from Standard Model physics: qualitative features*, JCAP **07** (2015) 022 [arXiv:1504.02569] [inSPIRE].

[13] J. Ghiglieri, G. Jackson, M. Laine and Y. Zhu, *Gravitational wave background from Standard Model physics: complete leading order*, JHEP **07** (2020) 092 [arXiv:2004.11392] [inSPIRE].

[14] B. Allen, *Stochastic gravity-wave background in inflationary-universe models*, Phys. Rev. D **37** (1988) 2078 [inSPIRE].

[15] S. Vagnozzi and A. Loeb, *The Challenge of Ruling Out Inflation via the Primordial Graviton Background*, Astrophys. J. Lett. **939** (2022) L22 [arXiv:2208.14088] [inSPIRE].

[16] J. Ghiglieri, J. Schütte-Engel and E. Speranza, *Freezing-in gravitational waves*, Phys. Rev. D **109** (2024) 023538 [arXiv:2211.16513] [inSPIRE].

[17] S.Y. Choi, J.S. Shim and H.S. Song, *Factorization and polarization in linearized gravity*, Phys. Rev. D **51** (1995) 2751 [hep-th/9411092] [inSPIRE].

[18] J. Kuipers, T. Ueda, J.A.M. Vermaseren and J. Vollinga, *FORM version 4.0*, Comput. Phys. Commun. **184** (2013) 1453 [arXiv:1203.6543] [inSPIRE].

[19] T. Hahn, *Generating Feynman diagrams and amplitudes with FeynArts 3*, Comput. Phys. Commun. **140** (2001) 418 [hep-ph/0012260] [inSPIRE].

[20] T. Hahn, S. Paßehr and C. Schappacher, *FormCalc 9 and Extensions*, PoS LL2016 (2016) 068 [arXiv:1604.04611] [inSPIRE].

[21] V. Shtabovenko, R. Mertig and F. Orellana, *FeynCalc 9.3: New features and improvements*, Comput. Phys. Commun. **256** (2020) 107478 [arXiv:2001.04407] [inSPIRE].

[22] N.E.J. Bjerrum-Bohr, B.R. Holstein, L. Planté and P. Vanhove, *Graviton-photon scattering*, Phys. Rev. D **91** (2015) 064008 [arXiv:1410.4148] [inSPIRE].

[23] D. Bödeker and D. Schröder, *Equilibration of right-handed electrons*, JCAP **05** (2019) 010 [arXiv:1902.07220] [inSPIRE].

[24] E. Braaten and R.D. Pisarski, *Soft amplitudes in hot gauge theories: A general analysis*, Nucl. Phys. B **337** (1990) 569 [inSPIRE].

[25] J.C. Taylor and S.M.H. Wong, *The effective action of hard thermal loops in QCD*, Nucl. Phys. B **346** (1990) 115 [inSPIRE].

[26] A. Ringwald, J. Schütte-Engel and C. Tamarit, *Gravitational waves as a big bang thermometer*, JCAP **03** (2021) 054 [arXiv:2011.04731] [inSPIRE].

[27] P. Klose, M. Laine and S. Procacci, *Gravitational wave background from non-Abelian reheating after axion-like inflation*, JCAP **05** (2022) 021 [arXiv:2201.02317] [inSPIRE].

[28] L. Castells-Tiestos and J. Casalderrey-Solana, *Thermal emission of gravitational waves from weak to strong coupling*, JHEP **10** (2022) 049 [arXiv:2202.05241] [inSPIRE].
[29] F. Muia, F. Quevedo, A. Schachner and G. Villa, *Testing BSM physics with gravitational waves*, JCAP 09 (2023) 006 [arXiv:2303.01548] [INSPIRE].

[30] N. Bernal, S. Cléry, Y. Mambrini and Y. Xu, *Probing reheating with graviton bremsstrahlung*, JCAP 01 (2024) 065 [arXiv:2311.12694] [INSPIRE].

[31] M. Drewes, Y. Georis, J. Klaric and P. Klose, *Upper Bound on Thermal Gravitational Wave Backgrounds from Hidden Sectors*, arXiv:2312.13855 [INSPIRE].

[32] PLANCK collaboration, *Planck 2018 results. X. Constraints on inflation*, Astron. Astrophys. 641 (2020) A10 [arXiv:1807.06211] [INSPIRE].

[33] M. Laine, M. Meyer and Y. Schröder, *Data for the Standard Model equation of state*, http://www.laine.itp.unibe.ch/eos15/.

[34] PLANCK collaboration, *Planck 2018 results. VI. Cosmological parameters*, Astron. Astrophys. 641 (2020) A6 [Erratum ibid. 652 (2021) C4] [arXiv:1807.06209] [INSPIRE].

[35] T. Daniel, M. Howard and M. König, *An SZ-like effect on cosmological gravitational wave backgrounds*, JCAP 12 (2023) 041 [arXiv:2308.00111] [INSPIRE].

[36] D.J.H. Chung, E.W. Kolb and A. Riotto, *Production of massive particles during reheating*, Phys. Rev. D 60 (1999) 063504 [hep-ph/9809453] [INSPIRE].

[37] H. Kolesova and M. Laine, *Update on gravitational wave signals from post-inflationary phase transitions*, Phys. Lett. B 851 (2024) 138553 [arXiv:2311.03718] [INSPIRE].