# Gravitational wave background from Standard Model physics: Complete leading order

J. Ghiglieri, G. Jackson, M. Laine, Y. Zhu

► To cite this version:

J. Ghiglieri, G. Jackson, M. Laine, Y. Zhu. Gravitational wave background from Standard Model physics: Complete leading order. Journal of High Energy Physics, 2020, 07, pp.092. 10.1007/JHEP07(2020)092. hal-02572522

HAL Id: hal-02572522

https://hal.science/hal-02572522

Submitted on 3 Aug 2021

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
RECEIVED: May 5, 2020
ACCEPTED: June 21, 2020
PUBLISHED: July 14, 2020

# Gravitational wave background from Standard Model physics: complete leading order

J. Ghiglieri,ª G. Jackson,ᵇ M. Laineᵇ and Y. Zhuᶜ,<fnref n="1" />

ªSUBATECH, Université de Nantes, IMT Atlantique, IN2P3/CNRS,
4 rue Alfred Kastler, La Chantrerie BP 20722, 44307 Nantes, France

ᵇAEC, Institute for Theoretical Physics, University of Bern,
Sidlerstrasse 5, CH-3012 Bern, Switzerland

ᶜWetzlar, Germany

E-mail: jacopo.ghiglieri@subatech.in2p3.fr, jackson@itp.unibe.ch,
laine@itp.unibe.ch, zhuyiopp@gmail.com

**ABSTRACT:** We compute the production rate of the energy density carried by gravitational waves emitted by a Standard Model plasma in thermal equilibrium, consistently to leading order in coupling constants for momenta $k \sim \pi T$. Summing up the contributions from the full history of the universe, the highest temperature of the radiation epoch can be constrained by the so-called $N_{\text{eff}}$ parameter. The current theoretical uncertainty $\Delta N_{\text{eff}} \le 10^{-3}$ corresponds to $T_{\text{max}} \le 2 \times 10^{17}$ GeV. In the course of the computation, we show how a subpart of the production rate can be determined with the help of standard packages, even if subsequently an IR subtraction and thermal resummation need to be implemented.

**KEYWORDS:** Thermal Field Theory, Cosmology of Theories beyond the SM, Resummation, Quark-Gluon Plasma

ARXIV ePRINT: 2004.11392

<footnote n="1">¹Previous address: Physik-Department, TU München, James-Franck-Strasse 1, 85748 Garching, Germany.</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Steps of the computation</td><td>2</td></tr><tr><td>2.1</td><td>Setup</td><td>2</td></tr><tr><td>2.2</td><td>Retarded energy-momentum correlator</td><td>4</td></tr><tr><td>2.3</td><td>Extracting 2 ↔ 2 cuts at light cone</td><td>7</td></tr><tr><td>2.4</td><td>Connection to Boltzmann equations</td><td>10</td></tr><tr><td>2.5</td><td>Phase space integrals</td><td>12</td></tr><tr><td>2.5.1</td><td>t-channel</td><td>12</td></tr><tr><td>2.5.2</td><td>s-channel</td><td>14</td></tr><tr><td>2.5.3</td><td>IR divergence</td><td>15</td></tr><tr><td>2.6</td><td>Hard Thermal Loop resummation</td><td>15</td></tr><tr><td>3</td><td>Numerical results</td><td>19</td></tr><tr><td>4</td><td>Cosmological implications</td><td>20</td></tr><tr><td>5</td><td>Conclusions and outlook</td><td>21</td></tr><tr><td>A</td><td>Soft t-channel fermion exchange</td><td>22</td></tr><tr><td>B</td><td>Magnitude of 1 + n ↔ 2 + n processes</td><td>24</td></tr></table>

## 1 Introduction

A neutral plasma with charged constituents, such as the early universe before recombination, emits and absorbs photons, because scatterings between the microscopic constituents amount to changing electromagnetic currents. Similarly, a homogeneous plasma can emit and absorb gravitational waves, because scatterings also imply changing energy and momentum currents (cf., e.g., ref. [1]). The emission/absorption rate is suppressed by $1/m_{Pl}^2$ and therefore tiny for temperatures much below the Planck scale. On the other hand, the age of the universe (inverse Hubble rate) is $\sim m_{Pl}$, so that the total energy density emitted into gravitational radiation is only suppressed by $1/m_{Pl}$. This may motivate a precise computation of the production rate and its integration over the history of the universe [2].

In addition to the emission from an equilibrium plasma, there are numerous potential non-equilibrium sources for gravitational radiation. These range from tensor modes produced during inflation [3] to a multitude of post-inflationary sources (for a review see, e.g., ref. [4]). However, all of these rely on yet-to-be-established models, unlike the Standard Model background that we are interested in.
Restricting for a moment to locally Minkowskian spacetime, the rate of change of the polarization-averaged phase space distribution of gravitons ($f_{GW}$) has the form [5]

$$ \dot{f}_{GW}(t, \mathbf{k}) = \Gamma(k) [n_B(k) - f_{GW}(t, \mathbf{k})] + \mathcal{O}\left(\frac{1}{m_{Pl}^4}\right), \qquad (1.1) $$

where $k \equiv |\mathbf{k}|$ and $n_B(k) \equiv 1/(e^{k/T} - 1)$ is the Bose distribution. The differential energy density is given by $de_{GW} = 2k f_{GW} \frac{d^3\mathbf{k}}{(2\pi)^3}$. Adopting a logarithmic scale, the production rate of gravitational energy density can thus be expressed as

$$ \frac{\mathrm{d}e_{\mathrm{GW}}}{\mathrm{d}t \mathrm{d} \ln k} = \frac{k^4 \dot{f}_{\mathrm{GW}}}{\pi^2}. \qquad (1.2) $$

In the following we are interested in estimating the rate $\Gamma(k)$ defined by eq. (1.1) in the frequency range in which $de_{GW}$ peaks. This range is given by the typical thermal scale $k \sim \pi T$ [2], corresponding after red shift to the same microwave range at which most CMB photons lie. In this frequency range, the gravitational wave abundance is expected to be much below equilibrium, $f_{GW} \ll n_B(k)$, so that the right-hand side of eq. (1.1) evaluates to $\Gamma(k)n_B(k)$. However, the same coefficient $\Gamma(k)$ also governs other phenomena, for instance the damping of a gravitational wave as it passes through a thermal plasma, if produced by some astrophysical source before (cf., e.g., refs. [6, 7] for recent works).

We start by describing in some detail the technical steps of the computation, which we have implemented in two complementary ways, viz. by taking the cut of a retarded 2-point correlator of the energy-momentum tensor (sections 2.1–2.3), and by considering Boltzmann equations for graviton production (section 2.4). After phase space integration (section 2.5) and thermal resummation (section 2.6), the result is evaluated numerically (section 3) and embedded in a cosmological environment (section 4). Conclusions and an outlook are offered in section 5. Two appendices explain why two classes of contributions, frequently considered in the literature, are of subleading order for the present observable.

# 2 Steps of the computation

## 2.1 Setup

Assuming that a system is spatially homogeneous and stationary on the time scales observed, and aligning the z-axis with the momentum ($\mathbf{k} = k \mathbf{e}_z$), the production rate of the energy density carried by gravitational waves can be related to the Wightman correlator

$$ G^{<}_{12;12} \equiv \int_{\mathcal{X}} e^{ik(t-z)} \langle T_{12}(0) T_{12}(\mathcal{X}) \rangle, \quad \mathcal{X} \equiv (t, \mathbf{x}). \qquad (2.1) $$

Here we work in the medium rest frame, with its four-velocity taking the form $u = (1, 0)$, in order to permit for a simple identification of the energy density. For a general frame, spatial indices $(...)_i$ should be replaced with $(g_i^\mu - u_i u^\mu)(...)_\mu$.

In equilibrium, $G^{<}_{12;12}$ is related to the imaginary part of the retarded correlator as $G^{<}_{12;12} = 2n_B(k) \mathrm{Im} G^{\mathrm{R}}_{12;12}$. In the following we compute a Euclidean correlator $G^{\mathrm{E}}_{12;12}$ as a function of a Euclidean four-momentum $K = (k_n, \mathbf{k})$, from which $G^{\mathrm{R}}_{12;12}$ is obtained by
an analytic continuation, $G_{12;12}^R = G_{12;12}^E|_{k_n \to -i[k_i + i0^+]}$. Here $k_n = 2\pi nT$, with $n \in \mathbb{Z}$, is a bosonic Matsubara frequency. The rate $\Gamma(k)$ from eq. (1.1) is then given by [2]

$$ \Gamma(k) = \frac{16\pi \operatorname{Im} G_{12;12}^R}{k m_{\mathrm{Pl}}^2}, \qquad (2.2) $$

where $m_{\mathrm{Pl}} = 1.22091 \times 10^{19}$ GeV is the Planck mass.

We write the correlator in a covariant form as<fnref n="1" />

$$ G_{12;12}^{\mathrm{E}} = \frac{L_{\mu\nu;\alpha\beta} G_{\mu\nu;\alpha\beta}^{\mathrm{E}}}{D(D-3)}, \quad G_{\mu\nu;\alpha\beta}^{\mathrm{E}} \equiv \int_X e^{iK \cdot X} \langle T_{\mu\nu}(X) T_{\alpha\beta}(0) \rangle, \qquad (2.3) $$

where $D$ denotes the dimension of space-time, $X \equiv (\tau, \mathbf{x})$, and $\tau \in (0, \frac{1}{T})$. Here we have defined the projector $(L_{\mu\nu;\alpha\beta} L_{\alpha\beta;\gamma\delta} = L_{\mu\nu;\gamma\delta})$

$$ L_{\mu\nu;\alpha\beta} \equiv \frac{\mathbb{P}_{\mu\alpha}^{\mathrm{T}} \mathbb{P}_{\nu\beta}^{\mathrm{T}} + \mathbb{P}_{\mu\beta}^{\mathrm{T}} \mathbb{P}_{\nu\alpha}^{\mathrm{T}}}{2} - \frac{\mathbb{P}_{\mu\nu}^{\mathrm{T}} \mathbb{P}_{\alpha\beta}^{\mathrm{T}}}{D-2}, \quad \mathbb{P}_{\mu\nu}^{\mathrm{T}} \equiv \delta_{\mu i} \delta_{\nu j} \left( \delta_{ij} - \frac{k_i k_j}{k^2} \right), \qquad (2.4) $$

which is symmetric $(L_{\mu\nu;\alpha\beta} = L_{\nu\mu;\alpha\beta} = L_{\alpha\beta;\mu\nu})$ and projects onto transverse $(K_\mu L_{\mu\nu;\alpha\beta} = k_i \delta_{i\mu} L_{\mu\nu;\alpha\beta} = 0)$ and traceless $(\delta_{\mu\nu} L_{\mu\nu;\alpha\beta} = 0)$ modes. We also denote

$$ \mathbb{P}_{\mathbf{p}}^{\mathrm{T}} \equiv \mathbb{P}_{\mu\nu}^{\mathrm{T}} P_{\mu} P_{\nu} = p^2 - \frac{(\mathbf{p} \cdot \mathbf{k})^2}{k^2}. \qquad (2.5) $$

As $T_{\mu\nu}$ we take the Standard Model energy-momentum tensor, which we write in Euclidean metric. Given that $L_{\mu\nu;\alpha\beta}$ projects out trace parts, it is enough to include non-trace ones,

$$ T_{\mu\nu} \supset F_{\mu\alpha}^{a_i} F_{\nu\alpha}^{a_i} + (D_{\mu}\phi)^{\dagger}(D_{\nu}\phi) + (D_{\nu}\phi)^{\dagger}(D_{\mu}\phi) \\ + \frac{1}{4} \left[ \bar{q}_{\mathrm{L}}(\gamma_{\mu}\overleftrightarrow{D}_{\nu} + \gamma_{\nu}\overleftrightarrow{D}_{\mu})q_{\mathrm{L}} + \bar{u}_{\mathrm{R}}(\gamma_{\mu}\overleftrightarrow{D}_{\nu} + \gamma_{\nu}\overleftrightarrow{D}_{\mu})u_{\mathrm{R}} + \bar{d}_{\mathrm{R}}(\gamma_{\mu}\overleftrightarrow{D}_{\nu} + \gamma_{\nu}\overleftrightarrow{D}_{\mu})d_{\mathrm{R}} \right. \\ \left. + \bar{\ell}_{\mathrm{L}}(\gamma_{\mu}\overleftrightarrow{D}_{\nu} + \gamma_{\nu}\overleftrightarrow{D}_{\mu})\ell_{\mathrm{L}} + \bar{\nu}_{\mathrm{R}}(\gamma_{\mu}\overleftrightarrow{D}_{\nu} + \gamma_{\nu}\overleftrightarrow{D}_{\mu})\nu_{\mathrm{R}} + \bar{e}_{\mathrm{R}}(\gamma_{\mu}\overleftrightarrow{D}_{\nu} + \gamma_{\nu}\overleftrightarrow{D}_{\mu})e_{\mathrm{R}} \right], \qquad (2.6) $$

where the $a_i$ label the generators of the various gauge groups; $\phi$ is the Higgs doublet; $q_L, \ell_L$ are the left-handed quark and lepton doublets, respectively; and $u_R, d_R, \nu_R, e_R$ are the corresponding right-handed components. The covariant derivative has the form

$$ D_{\mu} = \partial_{\mu} - ig_{1}Y A_{\mu} - ig_{2}T^{a2} A_{\mu}^{a2} a_{L} - ig_{3}T^{a3} A_{\mu}^{a3}, \qquad (2.7) $$

where $g_1, g_2, g_3$ are gauge couplings, $a_L$ is the left-handed projector and the hypercharge assignments are $Y = -\frac{1}{2}, -\frac{1}{2N_c}, -\frac{N_c+1}{2N_c}, \frac{N_c-1}{2N_c}, \frac{1}{2}, 0, 1$ for $\phi, q_L, u_R, d_R, \ell_L, \nu_R, e_R$, respectively [8]. We note that because of their vanishing gauge charge assignments and the omission of their Yukawa couplings, the fields $\nu_R$ do not contribute to $2 \leftrightarrow 2$ scatterings and have thus no effect on our final results (traditionally, $\nu_R$ are often omitted from the outset).

<footnote n="1">¹A simple way to verify the factor in the denominator is to consider momentum averages in the transverse plane. By rotational symmetry, $\langle q_i q_j q_k q_l \rangle = A (\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})$. Therefore a representative of $\langle T_{12}T_{12} \rangle$ evaluates to $\langle q_1^2 q_2^2 \rangle = A$, whereas $L_{ij;kl} \langle q_i q_j q_k q_l \rangle = AD(D-3)$.</footnote>
In order to avoid inverse polynomials of $D$ in section 2.2, the result for $G_{12;12}^E$ is expressed as

$$ G_{12;12}^E = \frac{2}{D(D-2)(D-3)} \left\{ \begin{aligned} & + n_S \Phi_s + 2n_G (1 + N_c) \Phi_f + (2 + N_c C_F) \Phi_g + n_S \lambda \Phi_{s(s)} + (3g_2^2 + N_c^2 C_F g_3^2) \Phi_{g(g)} \\ & + n_S |h_t|^2 N_c [\Phi_{s(f)} + \Phi_{f(s)} + \Phi_{s|f}] + n_S (g_1^2 + 3g_2^2) [\Phi_{s(g)} + \Phi_{g(s)} + \Phi_{s|g}] \\ & + n_G \left[ \frac{(N_c + 1)(N_c + 2)g_1^2}{4N_c} + \frac{3(N_c + 1)g_2^2}{4} + 2N_c C_F g_3^2 \right] [\Phi_{f(g)} + \Phi_{g(f)} + \Phi_{f|g}] \\ & + \mathcal{O}(g^4) \end{aligned} \right\}, \quad (2.8) $$

where $n_S = 1$ is the number of Higgs doublets, $n_G = 3$ is the number of fermion generations, $C_F \equiv (N_c^2 - 1)/(2N_c)$, and $\mathcal{O}(g^4)$ refers generically to any 3-loop contribution.<fnref n="2" /> Here $s, f, g$ refer to effects from scalars, fermions, and gauge bosons, respectively; $\Phi_a$ is a 1-loop diagram with a particle of type a; $\Phi_{a(b)}$ is a 2-loop diagram where a particle of type a couples to $T_{\mu\nu}$ and a particle of type b appears in a loop; and $\Phi_{a|b}$ is a 2-loop diagram involving a cross correlation between the energy-momentum tensors of particles of types a and b (in terms of matrix elements this corresponds to an interference term). The corresponding Feynman diagrams are shown in figure 1.

## 2.2 Retarded energy-momentum correlator

As the gravitational wave production rate is dominated by very high temperatures, we treat all particles as massless for the moment (the role of thermal masses is discussed in section 2.6 and in appendices A and B). Then the results for the correlators can be expressed in terms of the “master” sum-integrals [9]

$$ J_{ab}^{c} \equiv \int_{P}^{c} j_{ab}^{c}, \qquad \tilde{J}_{ab}^{c} \equiv \int_{\{P\}}^{c} j_{ab}^{c}, \qquad j_{ab}^{c} \equiv \frac{[\mathbb{P}_{\mathbf{p}}^{\mathrm{T}}]^{c} [K^{2}]^{x}}{[P^{2}]^{a} [(K-P)^{2}]^{b}}, \quad (2.9) $$

$$ I_{abcde}^{fgh} \equiv \int_{PQ}^{c} i_{abcde}^{fgh}, \qquad \tilde{I}_{abcde}^{fgh} \equiv \int_{P\{Q\}}^{c} i_{abcde}^{fgh}, \qquad \tilde{I}_{abcde}^{fgh} \equiv \int_{\{P\}Q}^{c} i_{abcde}^{fgh}, \qquad \tilde{I}_{abcde}^{fgh} \equiv \int_{PQ}^{c} i_{abcde}^{fgh}, $$

$$ i_{abcde}^{fgh} \equiv \frac{[\mathbb{P}_{\mathbf{p}}^{\mathrm{T}}]^{f} [\mathbb{P}_{\mathbf{q}}^{\mathrm{T}}]^{g} [\mathbb{P}_{\mathbf{q}-\mathbf{p}}^{\mathrm{T}}]^{h} [K^{2}]^{y}}{[P^{2}]^{a} [Q^{2}]^{b} [(Q-P)^{2}]^{c} [(K-P)^{2}]^{d} [(K-Q)^{2}]^{e}}, \quad (2.10) $$

where $\{P\}$ denotes a fermionic Matsubara four-momentum. The indices $x \equiv a + b - c$ and $y \equiv a + b + c + d + e - f - g - h - 2$ guarantee the overall dimensionality GeV$^4$. In the fermionic cases the representation is not unique; for the class of masters discussed in section 2.3, which have a cut corresponding to a $2 \leftrightarrow 2$ scattering, we have ordered the indices such that $a, c, e$ are non-negative.

The reduction of the energy-momentum tensor correlator to the basis of eqs. (2.9) and (2.10) has been carried out with a self-designed algorithm implemented in FORM [10].

<footnote n="2">²The Higgs self-coupling and top Yukawa coupling appear in a Euclidean Lagrangian as $L_E \supset \lambda(\phi^\dagger\phi)^2 + \bar{q}_L h_t t_R \phi + \bar{\phi}^\dagger \bar{t}_R h_t^* q_L$, whereas other Yukawa couplings are omitted.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The 1 and 2-loop graphs contributing to eq. (2.8). Each subset is gauge independent. Dashed lines denote scalars; solid lines fermions; wiggly lines gauge fields; dotted lines ghosts; blobs the operator $T_{\mu\nu}$. Graphs obtained by symmetrizations have been omitted.</figcaption></figure>

After the use of symmetries related to substitutions of integration variables, and noting that terms with odd numbers of $\gamma_5$-matrices do not contribute at this order, the results read

$$ \Phi_s = 4(D-3)J_{11}^2, \qquad (2.11) $$

$$ \Phi_f = -4(D-3)\tilde{J}_{11}^2 + \frac{D(D-3)}{2}(2\tilde{J}_{10}^1 - \tilde{J}_{11}^1), \qquad (2.12) $$

$$ \Phi_g = 2(D-3) \left[ (D-2)J_{11}^2 + D(J_{11}^1 - J_{10}^1) + \frac{D(D-2)}{8} (J_{11}^0 - 2J_{10}^0 + 4J_{00}^0) \right], \qquad (2.13) $$

$$ \Phi_{s(s)} = -48(D-3)I_{21010}^{200}, \qquad (2.14) $$

$$
\begin{aligned}
\Phi_{g(g)} = \frac{D(D-2)(D-3)}{2} & \left[ -I_{11111}^{000} - I_{12101}^{010} - I_{21100}^{100} + 3I_{10101}^{100} - 12I_{12001}^{010} \right. \\
& \left. + 2(I_{121-21}^{010} - I_{11100}^{000} - I_{11001}^{010} - I_{11000}^{000} - I_{20010}^{100}) \right. \\
& \left. + 4(I_{121-11}^{100} + I_{111-11}^{100} + I_{11101}^{100} + I_{11101}^{000} - I_{21100}^{010} - I_{11001}^{100} - I_{12001}^{100}) \right] \\
& + 2D(D-3) \left[ 4I_{11101}^{010} - 2I_{11111}^{100} - I_{11111}^{001} \right] + 2D(D-6) \left[ 2I_{11111}^{101} + I_{11111}^{110} \right] \\
& - (3D^2 - 16D + 12) \left[ 2I_{11111}^{200} + I_{11111}^{002} \right] - \frac{D(D-3)(3D-10)}{2} I_{11100}^{100} \\
& + D(D-2) \left[ 4(I_{12101}^{110} - I_{21100}^{110}) + 2(I_{12101}^{101} - I_{12101}^{200} + I_{21100}^{020} - I_{21100}^{011}) \right. \\
& \left. + I_{21100}^{200} - I_{12101}^{020} \right] + \frac{(D-2)^2(D-3)}{2} \left[ D(I_{12000}^{000} - I_{12000}^{000}) - 8I_{12001}^{020} \right],
\end{aligned}
\qquad (2.15)
$$
$$ \Phi_{s(f)} = 8(D-3)[2\hat{I}_{12001}^{020} - \hat{I}_{11101}^{020}], \quad (2.16) $$

$$ \begin{aligned} \Phi_{f(s)} = \frac{D(D-3)}{2} [ & 4(\bar{I}_{111-11}^{010} + \bar{I}_{10101}^{100} - \bar{I}_{10110}^{100} - \bar{I}_{11100}^{100} - \bar{I}_{20100}^{100} + \bar{I}_{20110}^{100} + \bar{I}_{21000}^{100} - \bar{I}_{21010}^{100}) \\ & -2(\bar{I}_{10101}^{001} + \bar{I}_{11011}^{100} + \bar{I}_{11011}^{001}) + 16(D-3)[\bar{I}_{20110}^{200} - \bar{I}_{21010}^{200}] - \frac{3D-8}{2}\bar{I}_{11111}^{002} \\ & +2(D-2)[4(\bar{I}_{11011}^{101} - \bar{I}_{11101}^{101} - \bar{I}_{11101}^{110}) + 2(\bar{I}_{11101}^{002} + \bar{I}_{11101}^{020} + \bar{I}_{11101}^{200} - \bar{I}_{11011}^{200}) \\ & + \bar{I}_{11111}^{101} - \bar{I}_{11011}^{002}] + (D-4)[8\bar{I}_{11101}^{011} - 4\bar{I}_{11011}^{110} + \bar{I}_{11111}^{200} - \bar{I}_{11111}^{110}], \quad (2.17) \end{aligned} $$

$$ \begin{aligned} \Phi_{s|f} = {}& 2(D-2)[4(\tilde{I}_{11101}^{011} + \tilde{I}_{11101}^{101}) - 2(\tilde{I}_{11111}^{011} + \tilde{I}_{11101}^{200} + \tilde{I}_{11101}^{020} + \tilde{I}_{11101}^{002}) + \tilde{I}_{11111}^{020} \\ & + \tilde{I}_{11111}^{002}] + 8(D-3)\tilde{I}_{11101}^{020} - 4[\tilde{I}_{11111}^{110} + \tilde{I}_{11111}^{101}] - 2(D-4)[\tilde{I}_{11111}^{200} + 4\tilde{I}_{11101}^{110}], \quad (2.18) \end{aligned} $$

$$ \begin{aligned} \Phi_{s(g)} = \frac{D-2}{2} [ & 4I_{11111}^{101} - 2I_{11111}^{200} - I_{11111}^{002} - (D-4)I_{11111}^{110} \\ & + \frac{D-3}{2}[3DI_{1010}^{010} - 8I_{11101}^{020} - 4(D-1)I_{21010}^{200}], \quad (2.19) \end{aligned} $$

$$ \begin{aligned} \Phi_{g(s)} = \frac{D(D-3)}{4} [ & 4(I_{121-11}^{100} - I_{111-11}^{100} + I_{11010}^{010} - I_{21100}^{010} - I_{21010}^{010}) \\ & + 2(I_{121-21}^{010} + I_{11101}^{010} - I_{20010}^{100}) - I_{21100}^{100} - I_{12101}^{010} - I_{10101}^{010} + 6I_{11010}^{100} \\ & + 7I_{11100}^{100} - 12I_{21010}^{100}] + 2(D-2)[I_{11101}^{200} - I_{11101}^{101}] - 4I_{11101}^{110} - (D-4)I_{11101}^{020} \\ & + \frac{D}{2}[4(I_{12101}^{110} - I_{21100}^{110}) + 2(I_{12101}^{101} - I_{12101}^{200} + I_{21100}^{020} - I_{21100}^{011}) + I_{21100}^{200} - I_{12101}^{020}] \\ & + \frac{(D-2)(D-3)}{4}[D(2I_{11010}^{000} + I_{21000}^{000} - I_{21010}^{000} - 4I_{11000}^{000}) - 8I_{21010}^{200}], \quad (2.20) \end{aligned} $$

$$ \begin{aligned} \Phi_{s|g} = {}& (D-2)[2(I_{11111}^{110} + I_{11101}^{101} - I_{11111}^{200} - I_{11101}^{200}) - I_{11111}^{002}] + 4I_{11111}^{101} \\ & + \frac{D(D-3)}{2}[4(I_{111-11}^{100} - I_{11010}^{100} - I_{11010}^{010}) + 2I_{11101}^{100} - I_{11101}^{010} - I_{10101}^{010} - 5I_{11100}^{100}] \\ & + 4I_{11101}^{110} + (5D-16)I_{11101}^{020} + \frac{D(D-2)(D-3)}{2}[2I_{11000}^{000} - I_{11010}^{000}], \quad (2.21) \end{aligned} $$

$$ \begin{aligned} \Phi_{f(g)} = \frac{D(D-2)(D-3)}{2} [ & 2(\bar{I}_{111-11}^{010} + \bar{I}_{101-11}^{000} + \bar{I}_{02101}^{010} - \bar{I}_{12001}^{010} - \bar{I}_{01101}^{010} + \bar{I}_{12000}^{010} - \bar{I}_{02100}^{010}) \\ & - \bar{I}_{10101}^{000} - \bar{I}_{01010}^{000}] + \frac{D(D-3)}{2} [\bar{I}_{11111}^{100} + \bar{I}_{11111}^{010} - \bar{I}_{11111}^{001} + 2\bar{I}_{11100}^{001} - 2D\bar{I}_{11100}^{100}] \\ & + \frac{(D-4)(D+2)}{4} [\bar{I}_{11111}^{200} + \bar{I}_{11111}^{020}] - \frac{3D^2-18D+32}{4}\bar{I}_{11111}^{002} - \frac{D^2-18D+40}{2}\bar{I}_{11111}^{110} \\ & + (D-2)^2 [2(\bar{I}_{11101}^{200} + \bar{I}_{11101}^{002} + \bar{I}_{11011}^{101} + \bar{I}_{11011}^{011}) - \bar{I}_{11011}^{200} - \bar{I}_{11011}^{020} - \bar{I}_{11011}^{002} - 4\bar{I}_{11101}^{101}] \\ & + \frac{D(D-3)(D-10)}{2} [2\bar{I}_{10101}^{100} - \bar{I}_{10101}^{001}] - \frac{D(D-3)(D-6)}{4} [\bar{I}_{11011}^{100} + \bar{I}_{11011}^{010} - \bar{I}_{11011}^{001}] \\ & + 2(D^2+4D-20)\bar{I}_{11101}^{020} - 2(D-2)(D-4)\bar{I}_{11011}^{110} + \frac{D^2-8D+20}{2}[\bar{I}_{11111}^{101} + \bar{I}_{11111}^{011}] \\ & + 4(D^2-10D+20)\bar{I}_{11101}^{011} - 4(D-4)^2\bar{I}_{11101}^{110} + 8(D-2)(D-3)[\bar{I}_{02101}^{020} - \bar{I}_{12001}^{020}], \quad (2.22) \end{aligned} $$

$$ \begin{aligned} \Phi_{g(f)} = \frac{D(D-2)(D-3)}{2} [ & \tilde{I}_{11100}^{000} - \hat{I}_{11101}^{000} - \tilde{I}_{101-11}^{000} \\ & + 2(\tilde{I}_{10101}^{000} + \tilde{I}_{21010}^{000} - \tilde{I}_{21000}^{000}) - 3\tilde{I}_{01100}^{000} - 4\tilde{I}_{11010}^{000} + 8\tilde{I}_{11000}^{000}] \end{aligned} $$
$$
\begin{align}
&+D(D-3)[2(\hat{I}_{111-11}^{100}-\hat{I}_{111-11}^{001}+\hat{I}_{121-11}^{001}-\hat{I}_{121-21}^{010}-\hat{I}_{121-11}^{100}-\hat{I}_{12101}^{100} \nonumber \\
&-\tilde{I}_{10101}^{001}-\tilde{I}_{21100}^{001}+\tilde{I}_{21000}^{001}-\tilde{I}_{21000}^{100}-\tilde{I}_{21000}^{010})+\hat{I}_{12101}^{010}+\tilde{I}_{21100}^{100}+3(\tilde{I}_{10101}^{100}+\tilde{I}_{11100}^{100}) \nonumber \\
&+4(\hat{I}_{11101}^{100}-\tilde{I}_{11010}^{100}+\tilde{I}_{21100}^{010})-6(\hat{I}_{11101}^{010}+\tilde{I}_{11100}^{010})+10\tilde{I}_{21010}^{100}] \nonumber \\
&+2D[4(\hat{I}_{21100}^{110}-\hat{I}_{12101}^{110})+2(\hat{I}_{12101}^{200}-\hat{I}_{12101}^{101}+\tilde{I}_{21100}^{011}-\tilde{I}_{21100}^{020})+\hat{I}_{12101}^{020}-\tilde{I}_{21100}^{200}] \nonumber \\
&+16\hat{I}_{11101}^{110}+8(D-2)[\hat{I}_{11101}^{101}-\hat{I}_{11101}^{200}] \nonumber \\
&+8(D-2)(D-3)\tilde{I}_{21010}^{200}-4(D^2-6D+10)\hat{I}_{11101}^{020}, \tag{2.23}
\end{align}
$$

$$
\begin{align}
\Phi_{f|g} ={}& D(D-2)(D-3)[\tilde{I}_{101-11}^{000} - \tilde{I}_{10101}^{000} + \tilde{I}_{01100}^{000} + 2\tilde{I}_{11010}^{000} - 4\tilde{I}_{11000}^{000}] \nonumber \\
& + D(D-3)[\tilde{I}_{11111}^{010} + \tilde{I}_{11111}^{001} - \tilde{I}_{11111}^{100} - 2\tilde{I}_{11100}^{100} + 3\tilde{I}_{10101}^{010} + 5\tilde{I}_{10101}^{001} - 6\tilde{I}_{10101}^{100} \nonumber \\
& + 4(\tilde{I}_{111-11}^{001} - \hat{I}_{111-11}^{100} + \hat{I}_{11101}^{010} - \tilde{I}_{11101}^{010} + \tilde{I}_{11101}^{100}) + 8(\tilde{I}_{11100}^{010} - \hat{I}_{11101}^{100})] \nonumber \\
& + 2(D^2 - 10D + 20)[2(\tilde{I}_{11011}^{110} - \tilde{I}_{11101}^{110}) - \tilde{I}_{11111}^{011}] + (D^2 - 2D - 4)[\tilde{I}_{11111}^{002} + \tilde{I}_{11111}^{020}] \nonumber \\
& + 8(D-2)[\hat{I}_{11101}^{200} - \hat{I}_{11101}^{101}] + 4(D^2 - 6D + 10)\hat{I}_{11101}^{020} - 16\hat{I}_{11101}^{110} \nonumber \\
& + 2(D-2)^2[2(\tilde{I}_{11101}^{101} - \tilde{I}_{11011}^{101}) + \tilde{I}_{11011}^{002} - \tilde{I}_{11101}^{002} + \tilde{I}_{11011}^{200} - \tilde{I}_{11101}^{200}] \nonumber \\
& + 2(D^2 - 12D + 28)\tilde{I}_{11011}^{020} - 2(D^2 + 4D - 20)\tilde{I}_{11101}^{020} \nonumber \\
& + (D-4)^2[4(\tilde{I}_{11101}^{011} - \tilde{I}_{11011}^{011}) - \tilde{I}_{11111}^{200}] - 2(3D-10)[\tilde{I}_{11111}^{101} + \tilde{I}_{11111}^{110}]. \tag{2.24}
\end{align}
$$

The computation was carried out in a general covariant gauge, and we have checked that the gauge parameter drops out exactly. The result for $\Phi_{g(g)}$ can be crosschecked against ref. [9].

## 2.3 Extracting 2 ↔ 2 cuts at light cone

As discussed below eq. (2.1), from each $\Phi$ we need to extract the cut $\mathrm{Im}\,\Phi|_{k_n \to -i[k+i0^+]}$. For the moment we only consider the cuts corresponding to 2 ↔ 2 scatterings, which originate from the masters $I$, with the discussion of 1 ↔ 2 reactions postponed to appendix B. As we restrict ourselves to the light cone, structures which have a positive power $y$ in eq. (2.10) yield no contribution. This implies that the only structures playing a role are of the types

$$ I_{101-11}^{000}, \quad I_{10101}^{100}, \quad I_{111-11}^{100}, \quad I_{121-21}^{100}, \quad I_{11101}^{200}. \qquad (2.25) $$

We denote the phase space of 2 ↔ 2 scatterings by

$$ \int d\Omega_{2\to2} \equiv \int \frac{d^3\mathbf{p}_1}{(2\pi)^3 2p_1} \int \frac{d^3\mathbf{p}_2}{(2\pi)^3 2p_2} \int \frac{d^3\mathbf{k}_1}{(2\pi)^3 2k_1} (2\pi)^4 \delta^{(4)}(\mathcal{P}_1 + \mathcal{P}_2 - \mathcal{K}_1 - \mathcal{K}_2), \quad (2.26) $$

where $\mathcal{P}_i \equiv (p_i, \mathbf{p}_i)$ with $p_i \equiv |\mathbf{p}_i|$, and $\mathcal{K}_2 \equiv \mathcal{K} \equiv (k, \mathbf{k})$. Distribution functions are denoted by

$$ n_{\sigma}(\epsilon) \equiv \frac{\sigma}{e^{\epsilon/T} - \sigma}, \quad \sigma = \pm, \qquad (2.27) $$

so that $n_+ = n_B$ and $n_- = -n_F$ are the Bose and Fermi distributions, respectively. Distribution functions appear in the combination

$$ \mathcal{N}_{\tau_1; \sigma_1 \sigma_2} \equiv n_{\tau_1}(k_1) [1 + n_{\sigma_1}(p_1)] [1 + n_{\sigma_2}(p_2)] - n_{\sigma_1}(p_1) n_{\sigma_2}(p_2) [1 + n_{\tau_1}(k_1)]. \quad (2.28) $$

Mandelstam variables are defined as usual, $s \equiv (\mathcal{P}_1 + \mathcal{P}_2)^2$, $t \equiv (\mathcal{P}_1 - \mathcal{K}_1)^2$, $u \equiv (\mathcal{P}_2 - \mathcal{K}_1)^2$.
With this notation, the 2 ↔ 2 cuts for the structures in eq. (2.25) read

$$
\mathrm{Im}\{I_{1b1d1}^{fgh}\}_{k_n \to -i[k+i0^+]} = \frac{1}{2} \int d\Omega_{2 \to 2} \left\{ \frac{[\mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}}]^f [\mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}}]^g [\mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}}]^h \mathcal{N}_{\sigma_a; \sigma_e \sigma_c}}{[-u]^b [-s]^d} \right. \\
\left. + \frac{[\mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}}]^f [\mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}}]^g [\mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}}]^h \mathcal{N}_{\sigma_c; \sigma_a \sigma_e}}{[-t]^b [-u]^d} \right. \\
\left. + \frac{[\mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}}]^f [\mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}}]^g [\mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}}]^h \mathcal{N}_{\sigma_e; \sigma_c \sigma_a}}{[-s]^b [-t]^d} \right\}, \quad (2.29)
$$

where $\sigma_a, \sigma_c$ and $\sigma_e$ label the statistics of the 1st, 3rd and 5th subscript of $I$, respectively.
The diagram illustrates the cuts, with crosses on the propagators $b$ and $d$ of which at least
one comes with a zero or negative power.

We can now collect together the cuts from eqs. (2.14)–(2.24). In so doing we also set
$D \to 4$ for simplicity, as there are no ultraviolet divergences in these cuts. Denoting by $\mathcal{C}$
an operation which produces an integrand for eq. (2.29), viz.

$$
\lim_{D \to 4} \mathrm{Im}\{\Phi\} |_{k_n \to -i[k+i0^+]}^{2 \leftrightarrow 2} \equiv \int d\Omega_{2 \to 2} \mathcal{C}\Phi, \quad (2.30)
$$

and making use of symmetries such as $\tilde{I}_{1b101}^{fgh} = \bar{I}_{1b101}^{hgf}$ (obtained by the substitution $P \to Q - P$), the non-zero contributions for the combinations appearing in eq. (2.8) read

$$
\mathcal{C}\Phi_{g(g)} = 4\mathcal{C}[\Phi_{s(g)} + \Phi_{g(s)} + \Phi_{s|g}] \quad (2.31)
$$

$$
= 4C [2I_{121-21}^{010} + 4I_{111-11}^{100} + 3I_{10101}^{100}]
$$

$$
= 2 \mathcal{N}_{+;+} \left\{ \mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}} \left( 3 + \frac{4u}{t} + \frac{2s^2}{u^2} \right) + \mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}} \left( 3 + \frac{4t}{s} + \frac{2u^2}{t^2} \right) \right. \\
\left. + \mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}} \left( 3 + \frac{4s}{u} + \frac{2t^2}{s^2} \right) \right\}, \quad (2.32)
$$

$$
\mathcal{C}[\Phi_{s(f)} + \Phi_{f(s)} + \Phi_{s|f}] = 4\mathcal{C}[2(\bar{I}_{111-11}^{010} + \bar{I}_{10101}^{100}) - \bar{I}_{10101}^{001}] \quad (2.31)
$$

$$
= 2 \mathcal{N}_{-;-} \left\{ \frac{2s \mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}}}{u} + 2 \mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}} - \mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}} \right\}
$$

$$
+ 2 \mathcal{N}_{-;-} \left\{ \frac{2t \mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}}}{s} + 2 \mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}} - \mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}} \right\}
$$

$$
+ 2 \mathcal{N}_{+-} \left\{ \frac{2u \mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}}}{t} + 2 \mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}} - \mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}} \right\}, \quad (2.33)
$$

$$
\mathcal{C}[\Phi_{f(g)} + \Phi_{g(f)} + \Phi_{f|g}] = 4\mathcal{C}[2(\tilde{I}_{111-11}^{001} - \tilde{I}_{111-11}^{100} - \tilde{I}_{121-21}^{010} + \tilde{I}_{111-11}^{010} + \tilde{I}_{101-11}^{000}) + \tilde{I}_{101-11}^{000}] \quad (2.34)
$$

$$
= 4 \mathcal{N}_{-;-} \left\{ \frac{s \mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}}}{u} + \frac{u [\mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}} - \mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}}]}{t} - \frac{u^2 \mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}}}{t^2} \right\}
$$

$$
+ 4 \mathcal{N}_{-;-} \left\{ \frac{t \mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}}}{s} + \frac{s [\mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}} - \mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}}]}{u} - \frac{s^2 \mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}}}{u^2} \right\}
$$

$$
+ 4 \mathcal{N}_{+-;--} \left\{ \frac{u \mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}}}{t} + \frac{t [\mathbb{P}_{\mathbf{p}_1}^{\mathrm{T}} - \mathbb{P}_{\mathbf{p}_2}^{\mathrm{T}}]}{s} - \frac{t^2 \mathbb{P}_{\mathbf{k}_1}^{\mathrm{T}}}{s^2} \right\}. \quad (2.34)
$$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** t-channel $2 \leftrightarrow 2$ scatterings contributing to gravitational wave production (further processes are obtained with u and s-channel reflections). The notation is as in figure 1, with the double line indicating a graviton. Up to numerical prefactors, the amplitudes squared originating from these processes, after summing over the physical polarization states of the gravitons and Standard Model particles, correspond to the cuts shown in eqs. (2.36)–(2.38) (cf. section 2.4).</figcaption></figure>

At the light cone, there is a further identity that has not been employed yet and that permits for a remarkable simplification of eqs. (2.32)–(2.34). Noting that for massless particles $u = 2(k \cdot p_1 - k p_1)$, and recalling that $\mathbb{P}_{\mathbf{p}_1}^T = (k p_1 - k \cdot \mathbf{p}_1)(k p_1 + k \cdot \mathbf{p}_1)/k^2$, we can make use of energy-momentum conservation to verify that

$$ \frac{\mathbb{P}_{\mathbf{k}_1}^T}{s} + \frac{\mathbb{P}_{\mathbf{p}_2}^T}{t} + \frac{\mathbb{P}_{\mathbf{p}_1}^T}{u} = -1. \qquad (2.35) $$

With this identity, combined with renamings $\mathbf{p}_1 \leftrightarrow \mathbf{p}_2$ as well as a repeated use of $s+t+u = 0$, all projectors $\mathbb{P}^T$ can be eliminated, and the cuts in eqs. (2.31)–(2.34) can be written in a form where the breaking of Lorentz invariance through the medium manifests itself only through the distribution functions $\mathcal{N}_{\tau_1; \sigma_1 \sigma_2}$:

$$ C\Phi_{g(g)} = 4C\left[\Phi_{s(g)} + \Phi_{g(s)} + \Phi_{s|g}\right] = 2\mathcal{N}_{+;+} \left\{-2\left(\frac{s^2+u^2}{t} + \frac{t^2}{s}\right)\right\}, \qquad (2.36) $$

$$ C\left[\Phi_{s(f)} + \Phi_{f(s)} + \Phi_{s|f}\right] = 2\mathcal{N}_{-;-} \{2t\} + 2\mathcal{N}_{+;-} \{s\}, \qquad (2.37) $$

$$ C\left[\Phi_{f(g)} + \Phi_{g(f)} + \Phi_{f|g}\right] = 4\mathcal{N}_{-;-} \left\{\frac{s^2+u^2}{t}\right\} + 4\mathcal{N}_{+;-} \left\{\frac{t^2}{s}\right\}. \qquad (2.38) $$

We note that eq. (2.36) could be written in a more symmetric form, but for later convenience we prefer to use the same structures as in eq. (2.38). Eqs. (2.36)–(2.38) correspond to amplitudes squared for processes illustrated in figure 2 (cf. section 2.4).

The drastic simplification that we have observed when going on the light-cone has a known precedent: it also takes place for photon production from a thermal medium. Furthermore, in that case it is well understood. The transverse correlator to which physical photons couple, $\mathrm{Im} G_T^R$, can be replaced by the full vector correlator, $\mathrm{Im} G_V^R = \mathrm{Im} G_T^R + \mathrm{Im} G_L^R$, because a Ward identity guarantees the vanishing of $\mathrm{Im} G_L^R$ for zero virtuality. We are not aware of a similar operator relation between the tensor channel correlator in
eq. (2.3) and one without any $\mathbb{P}^T$'s, even if intriguing relations between photon and graviton production amplitudes are known to exist (cf. section 2.4).

## 2.4 Connection to Boltzmann equations

The $2 \leftrightarrow 2$ cuts of section 2.3 can also be obtained from kinetic theory and Boltzmann equations. As a starting point, we may, for $k \sim \pi T$, write the leading-order contribution to eq. (1.1) as

$$ \dot{f}_{\text{GW}}(t, \mathbf{k}) = \Gamma(k) n_{\text{B}}(k) = \frac{1}{8k} \int d\Omega_{2 \to 2} \sum_{abc} |\mathcal{M}_{cG}^{ab}(\mathbf{p}_1, \mathbf{p}_2; \mathbf{k}_1, \mathbf{k})|^2 f_a(p_1) f_b(p_2) [1 \pm f_c(k_1)], \quad (2.39) $$

where we have neglected $f_{\text{GW}}(t, \mathbf{k})$ on the right-hand side. The sum runs over all $abc \in \text{SM}$ (Standard Model) particle and antiparticle degrees of freedom and thus over all $ab \to cG$ processes, with $G$ denoting the graviton. $|\mathcal{M}_{cG}^{ab}(\mathbf{p}_1, \mathbf{p}_2; \mathbf{k}_1, \mathbf{k})|^2$ is the corresponding matrix element squared, summed over all degeneracies of each species. For the SM in the symmetric phase, these are spin, polarization, colour, weak isospin and generation. For $k \sim \pi T$ the contribution of thermal masses is suppressed, so the external states can be considered massless (thermal masses are only needed for the IR-divergent part of the squared amplitudes, cf. section 2.6). The prefactor $1/8k$ is a combination of $1/2k$ from the phase space measure, $1/2$ for the graviton polarization degeneracy, and $1/2$ for the symmetry factor for identical initial state particles; in the cases where $a \neq b$ this factor is compensated for by their being counted twice in the sum over $abc$. The thermal distributions $f_i$ correspond to $n_{\text{B}}$ and $n_{\text{F}}$ for bosons and fermions, respectively, with $[1 \pm f_c(k_1)]$ implying $[1 + n_{\text{B}}(k_1)]$ in the former case and $[1 - n_{\text{F}}(k_1)]$ in the latter.

The main challenge is the determination of the matrix elements squared, which requires the derivation of Feynman rules for all graviton-SM couplings and the computation of the tree-level amplitudes. Given the large number of vertices and processes, and the associated opportunities for error, we have adopted automated techniques, originally developed for collider physics. We first used FEYNRULES [11], which can derive Feynman rules from a given Lagrangian. We applied it to the Lagrangian describing the symmetric-phase SM coupled to gravitons, i.e.

$$ \mathcal{L}_{\text{SM+G}} = \mathcal{L}_{\text{SM}} - \frac{\sqrt{32\pi}}{2m_{\text{Pl}}} h_{\mu\nu} T_{\text{SM}}^{\mu\nu}, \quad (2.40) $$

where the SM energy-momentum tensor $T_{\text{SM}}^{\mu\nu}$ contains also the trace part. The kinetic term for gravitons can be omitted, as they are external states in our computation.

Using the appropriate interface [12], FEYNRULES can generate a *model file* for FEYN-ARTS [13] (unfortunately, sometimes manual fixes of the generation and SU(2) index assignments were needed). This package and its companion FORMCALC [14] were then used to generate, evaluate and square all amplitudes, summing over the relevant degeneracies.<fnref n="3" /> The handling of spin, vector boson polarization and colour is available in FORMCALC, whereas SU(2) algebra and tensor boson polarization had to be implemented. For the latter, we

<footnote n="3">³We have also looked into several other packages, however have not identified a procedure that would be simpler than the one described here.</footnote>
proceeded as follows. FEYNARTS assigns to external tensor bosons a polarization tensor $\epsilon_{\mu\nu}^{\lambda}(\mathbf{k})$ which is written, using a common factorization formula (cf., e.g., refs. [15, 16]), as

$$ \epsilon_{\mu\nu}^{\lambda}(\mathbf{k}) \equiv \epsilon_{\mu}^{\lambda}(\mathbf{k}) \epsilon_{\nu}^{\lambda}(\mathbf{k}), \qquad (2.41) $$

with $\epsilon_{\mu}^{\lambda}(\mathbf{k})$ the transverse polarization vector of a massless gauge boson. Upon taking $\mathbf{k} = k \mathbf{e}_z$ and the circular polarization vectors $\epsilon_{\mu}^{\lambda}(\mathbf{k}) = 1/\sqrt{2}(0, -1^{\lambda}, -i, 0)$, $\lambda = 1, 2$, it is easy to verify that the polarization sum satisfies

$$ \sum_{\lambda} \epsilon_{\mu\nu}^{\lambda}(\mathbf{k}) \epsilon_{\alpha\beta}^{*}(\mathbf{k}) = L_{\mu\nu;\alpha\beta}, \qquad (2.42) $$

with $L$ as defined in eq. (2.4). We implemented this form of the tensor polarization sum as a MATHEMATICA routine interfaced with the MATHEMATICA output of FEYNARTS/FORMCALC. The resulting matrix elements have an apparent dependence on the projectors $\mathbb{P}^T$, which again disappears by applying eq. (2.35).

Upon generating and evaluating all processes and plugging the results in eq. (2.39), we find

$$ \Gamma(k) n_B(k) = \frac{1}{8k} \frac{32\pi}{m_{Pl}^2} \int d\Omega_{2\to2} \left\{ \right. \\ + n_B(p_1) n_B(p_2) [1+n_B(k_1)] (g_1^2+15g_2^2+48g_3^2) \left( \frac{st}{u} + \frac{su}{t} + \frac{tu}{s} \right) \qquad (2.43) $$

$$ -n_F(p_1) n_B(p_2) [1-n_F(k_1)] \left[ 6|h_t|^2 t + (10g_1^2+18g_2^2+48g_3^2) \frac{s^2+u^2}{t} \right] \qquad (2.44) $$

$$ -n_B(p_1) n_F(p_2) [1-n_F(k_1)] \left[ 6|h_t|^2 u + (10g_1^2+18g_2^2+48g_3^2) \frac{s^2+t^2}{u} \right] \qquad (2.45) $$

$$ +n_F(p_1) n_F(p_2) [1+n_B(k_1)] \left[ 6|h_t|^2 s + (10g_1^2+18g_2^2+48g_3^2) \frac{t^2+u^2}{s} \right] \}. \qquad (2.46) $$

This expression agrees with the one obtained by plugging eqs. (2.36)-(2.38) into eqs. (2.30), (2.8) and (2.2). To verify the agreement, relabellings $\mathbf{p}_1 \leftrightarrow \mathbf{p}_2$ (and $t \leftrightarrow u$) as well as use of the identity $N_{\tau_1;\sigma_1\sigma_2} = n_{\sigma_1}(p_1) n_{\sigma_2}(p_2) [1 + n_{\tau_1}(k_1)] n_{\tau_1\sigma_1\sigma_2}^{-1}(p_1 + p_2 - k_1)$ are needed.

In obtaining the fermionic parts of the total rate, i.e. eqs. (2.44)-(2.46), we have not written out terms which arise from an odd number of $\gamma_5$ matrices in Dirac traces, since they vanish under the $\int d\Omega_{2\to2}$ integration. Specifically, these terms appear in the $fg \to fG$ processes and their crossings, with $f$ a fermion and $g$ a gauge boson.

We also note that the automated procedure fixes the gauge group factors, multiplicities and charge assignments to those specific for the SM; the coefficients multiplying the coupling constants are not obtained in terms of $N_c$, $n_G$ and $n_S$. Focussing on sub-processes, it is easy to reinstate group theory factors. For instance, the $g_3^2$-part of eq. (2.43) corresponds to the matrix elements squared for the gluonic scattering $gg \to gG$, yielding

$$ |\mathcal{M}_{gG}^{gg}(\mathbf{p}_1, \mathbf{p}_2; \mathbf{k}_1, \mathbf{k})|^2 = \frac{32\pi}{m_{Pl}^2} 2(N_c^2 - 1) N_c g_3^2 \left( \frac{st}{u} + \frac{su}{t} + \frac{tu}{s} \right). \qquad (2.47) $$

Recently, there has been much work on factorizing graviton amplitudes into photon amplitudes multiplied by kinematic factors, say $f\gamma \to fG$ versus $f\gamma \to f\gamma$ (cf., e.g.,
refs. [15, 16] and references therein). It is not clear to us, however, whether all the terms in eqs. (2.43)-(2.46) could be related to photon production or scattering rates.

We conclude this section by stressing that kinetic theory and its automated implementation are not sufficient for determining the leading-order gravitational wave production rate. Indeed, as discussed in sections 2.5.3 and 2.6, phase space integrals over matrix elements squared lead to IR divergences, related to soft gauge-boson exchange. The divergences need to be subtracted and subsequently Hard Thermal Loop resummed. An even more dramatic departure from the simple scattering picture is needed at smaller momenta, $k \sim \alpha_s^2 T$, where elementary particle states need to be replaced by hydrodynamic modes [2].

## 2.5 Phase space integrals

The next step is to carry out the phase space integral $\int d\Omega_{2\to2}$ for the cuts in eqs. (2.36)-(2.38) or the matrix elements squared in eqs. (2.43)-(2.46). For this task it is helpful to employ the parametrization introduced in ref. [17].<fnref n="4" /> We discuss separately the treatment of $t$ and $s$-channel cases ($u$-channel can always be transformed into $t$-channel).

### 2.5.1 $t$-channel

Consider the phase space integral

$$ \Gamma_{\tau_1; \sigma_1 \sigma_2}^t \equiv \int d\Omega_{2\to2} \mathcal{N}_{\tau_1; \sigma_1 \sigma_2} \left\{ a_1 \frac{s^2 + u^2}{t} + a_2 t \right\}. \qquad (2.48) $$

The idea is to insert $1 = \int d^4Q \delta^{(4)}(P_1 - K_1 - Q)$ in the integral. Then the energy-momentum conservation constraint inside $d\Omega_{2\to2}$ can be written as $\delta^{(4)}(Q + P_2 - K)$. We can now integrate over $p_2$ and $k_1$ by using the spatial parts of the Dirac $\delta$'s, leaving $q_0, q$ and $p_1$ as the integration variables. The temporal Dirac $\delta$'s fix two angles as

$$ q \cdot k = \frac{q^2 - q_0^2 + 2kq_0}{2}, \quad q \cdot p_1 = \frac{q^2 - q_0^2 + 2p_1q_0}{2}, \qquad (2.49) $$

whereas kinematic variables become

$$ t = q_0^2 - q^2, \quad u = 2(k \cdot p_1 - k p_1), \quad s = -t - u. \qquad (2.50) $$

The azimuthal average of powers of $k \cdot p_1$ can be computed by parametrizing

$$ q = (0, 0, q), \quad k = k(\sin \chi, 0, \cos \chi), \quad p_1 = p_1(\sin \theta \cos \varphi, \sin \theta \sin \varphi, \cos \theta), \qquad (2.51) $$

and integrating over $\varphi$. Denoting $\langle \dots \rangle \equiv \frac{1}{2\pi} \int_0^{2\pi} d\varphi (\dots)$, this yields

$$ \langle k \cdot p_1 \rangle = \frac{(q \cdot k) (q \cdot p_1)}{q^2}, \qquad (2.52) $$

$$ \langle (k \cdot p_1)^2 \rangle = \frac{1}{2} \left[ k^2 p_1^2 - \frac{p_1^2 (q \cdot k)^2}{q^2} - \frac{k^2 (q \cdot p_1)^2}{q^2} + \frac{3 (q \cdot k)^2 (q \cdot p_1)^2}{q^4} \right]. \qquad (2.53) $$

<footnote n="4">⁴If one is considering spectral functions off the light cone, more complicated structures $\sim P^T K^4/(ut)$ appear, which require a refined parametrization if a two-dimensional integral representation is desired [18, 19].</footnote>
The scalar products appearing here can be eliminated through eq. (2.49). Finally, the phase space distributions from eq. (2.28) can be cast in the form

$$ \begin{aligned} \mathcal{N}_{\tau_1; \sigma_1 \sigma_2} &= [1 + n_{\tau_1 \sigma_1}(p_1 - k_1) + n_{\sigma_2}(p_2)] [n_{\tau_1}(k_1) - n_{\sigma_1}(p_1)] \\ &= [1 + n_{\tau_1 \sigma_1}(q_0) + n_{\sigma_2}(k - q_0)] [n_{\tau_1}(p_1 - q_0) - n_{\sigma_1}(p_1)], \end{aligned} \quad (2.54) $$

thereby factorizing the $p_1$-dependence.

Denoting

$$ q_{\pm} \equiv \frac{q_0 \pm q}{2}, \quad (2.55) $$

the integration range of $p_1$ can be established as $(q_+, \infty)$. The integration measure contains no powers of $p_1$, whereas azimuthal averages yield powers up to $p_1^2$. The integral reads

$$ \begin{aligned} & \int_{q_+}^{\infty} \mathrm{d}p_1 (\beta_0 + \beta_1 p_1 + \beta_2 p_1^2) [n_{\tau_1}(p_1 - q_0) - n_{\sigma_1}(p_1)] \\ & = (\beta_0 + \beta_1 q_+ + \beta_2 q_+^2) L_1 + (\beta_1 + 2\beta_2 q_+) L_2 + (2\beta_2) L_3, \end{aligned} \quad (2.56) $$

where

$$ L_1 \equiv T[\ln(1 - \sigma_1 e^{-q_+/T}) - \ln(1 - \tau_1 e^{q_-/T})], \quad (2.57) $$

$$ L_2 \equiv T^2[\mathrm{Li}_2(\tau_1 e^{q_-/T}) - \mathrm{Li}_2(\sigma_1 e^{-q_+/T})], \quad (2.58) $$

$$ L_3 \equiv T^3[\mathrm{Li}_3(\tau_1 e^{q_-/T}) - \mathrm{Li}_3(\sigma_1 e^{-q_+/T})]. \quad (2.59) $$

All in all this results in

$$ \Gamma_{\tau_1; \sigma_1 \sigma_2}^{t} = \frac{1}{(4\pi)^3 k} \int_{-\infty}^{k} \mathrm{d}q_0 \int_{|q_0|}^{2k-q_0} \mathrm{d}q [1 + n_{\tau_1 \sigma_1}(q_0) + n_{\sigma_2}(k-q_0)] (q^2 - q_0^2) \\ \times \left\{ \frac{a_1[q^2 - 3(q_0 - 2k)^2][12L_3 + 6qL_2 + q^2L_1]}{6q^4} - \left(a_2 + \frac{2a_1}{3}\right)L_1 \right\}. \quad (2.60) $$

The integral in eq. (2.60) is logarithmically IR divergent at small $q_0$, $q$. For the different statistics the divergent parts read

$$ \begin{aligned} \Gamma_{+; ++|_{\mathrm{IR}}}^{t} &\equiv -2\Gamma_{-; -+|_{\mathrm{IR}}}^{t} \equiv \frac{1}{(4\pi)^3 k} \int_{-\infty}^{k} \mathrm{d}q_0 \int_{|q_0|}^{2k-q_0} \mathrm{d}q [1 + n_{\mathrm{B}}(q_0) + n_{\mathrm{B}}(k-q_0)] \\ &\quad \times \left\{ -\frac{4a_1 q_0 (q^2 - q_0^2) k^2 \pi^2 T^2}{q^4} \right\}, \end{aligned} \quad (2.61) $$

$$ \begin{aligned} \Gamma_{-; -+|_{\mathrm{IR}}}^{t} &\equiv -\Gamma_{+; --|_{\mathrm{IR}}}^{t} \equiv \frac{1}{(4\pi)^3 k} \int_{-\infty}^{k} \mathrm{d}q_0 \int_{|q_0|}^{2k-q_0} \mathrm{d}q [1 - n_{\mathrm{F}}(q_0) - n_{\mathrm{F}}(k-q_0)] \\ &\quad \times \left\{ \frac{42a_1(q^2 - q_0^2)k^2\zeta(3)T^3}{q^4} \right\}. \end{aligned} \quad (2.62) $$
### 2.5.2 s-channel

The s-channel phase space integral is defined as

$$ \Gamma_{\tau_1; \sigma_1 \sigma_2}^{s} \equiv \int d\Omega_{2 \to 2} \mathcal{N}_{\tau_1; \sigma_1 \sigma_2} \left\{ b_1 \frac{t^2}{s} + b_2 s \right\}. \qquad (2.63) $$

This time we insert $1 = \int d^4Q \delta^{(4)}(P_1 + P_2 - Q)$ in the integral, whereby the energy-momentum conservation constraint inside $d\Omega_{2 \to 2}$ can be written as $\delta^{(4)}(Q - K_1 - K)$. We integrate over $\mathbf{p}_1$ and $\mathbf{k}_1$ by using the spatial parts of the Dirac $\delta$'s, leaving $q_0, \mathbf{q}$ and $\mathbf{p}_2$ as the integration variables. The temporal Dirac $\delta$'s fix two angles as

$$ \mathbf{q} \cdot \mathbf{k} = \frac{q^2 - q_0^2 + 2kq_0}{2}, \quad \mathbf{q} \cdot \mathbf{p}_2 = \frac{q^2 - q_0^2 + 2p_2q_0}{2}, \qquad (2.64) $$

whereas kinematic variables become

$$ s = q_0^2 - q^2, \quad t = 2(\mathbf{k} \cdot \mathbf{p}_2 - k p_2), \quad u = -s - t. \qquad (2.65) $$

The azimuthal average of powers of $\mathbf{k} \cdot \mathbf{p}_2$ can be computed like in eqs. (2.52)-(2.53), exchanging $\mathbf{p}_1 \leftrightarrow \mathbf{p}_2$. The phase space distributions from eq. (2.28) are now cast in the form

$$ \begin{aligned} \mathcal{N}_{\tau_1; \sigma_1 \sigma_2} &= [1 + n_{\sigma_1}(p_1) + n_{\sigma_2}(p_2)] [n_{\tau_1}(k_1) - n_{\sigma_1 \sigma_2}(p_1 + p_2)] \\ &= [1 + n_{\sigma_1}(q_0 - p_2) + n_{\sigma_2}(p_2)] [n_{\tau_1}(q_0 - k) - n_{\sigma_1 \sigma_2}(q_0)], \end{aligned} \qquad (2.66) $$

factorizing the dependence on $p_2$. The integration range of $p_2$ can be established as $(q_-, q_+)$, and powers up to $p_2^2$ appear, whereby the general integral reads

$$ \begin{aligned} & \int_{q_{-}}^{q_{+}} d p_{2} (\beta_{0}+\beta_{1} p_{2}+\beta_{2} p_{2}^{2}) [1+n_{\sigma_{1}}(q_{0}-p_{2})+n_{\sigma_{2}}(p_{2})] \\ & = \beta_{0} q+\frac{\beta_{1} q q_{0}}{2}+\frac{\beta_{2} q(q^{2}+3 q_{0}^{2})}{12} \\ & \quad-(\beta_{0}+\beta_{1} q_{+}+\beta_{2} q_{+}^{2}) L_{1}^{+}-(\beta_{1}+2 \beta_{2} q_{+}) L_{2}^{+}-(2 \beta_{2}) L_{3}^{+} \\ & \quad+(\beta_{0}+\beta_{1} q_{-}+\beta_{2} q_{-}^{2}) L_{1}^{-}+(\beta_{1}+2 \beta_{2} q_{-}) L_{2}^{-}+(2 \beta_{2}) L_{3}^{-}, \end{aligned} \qquad (2.67) $$

where

$$ L_1^{\pm} \equiv T \left[ \ln \left(1 - \sigma_1 e^{-q_{\mp}/T}\right) - \ln \left(1 - \sigma_2 e^{-q_{\pm}/T}\right) \right], \qquad (2.68) $$

$$ L_2^{\pm} \equiv T^2 \left[ \mathrm{Li}_2 \left( \sigma_2 e^{-q_{\pm}/T} \right) + \mathrm{Li}_2 \left( \sigma_1 e^{-q_{\pm}/T} \right) \right], \qquad (2.69) $$

$$ L_3^{\pm} \equiv T^3 \left[ \mathrm{Li}_3 \left( \sigma_2 e^{-q_{\pm}/T} \right) - \mathrm{Li}_3 \left( \sigma_1 e^{-q_{\pm}/T} \right) \right]. \qquad (2.70) $$

All in all, this gives

$$ \Gamma_{\tau_1; \sigma_1 \sigma_2}^{s} = \frac{1}{(4\pi)^3 k} \int_k^\infty dq_0 \int_{|2k-q_0|}^{q_0} dq [n_{\tau_1}(q_0-k) - n_{\sigma_1 \sigma_2}(q_0)] (q^2 - q_0^2) \\ \times \left\{ \frac{b_1[q^2 - 3(q_0 - 2k)^2][12(L_3^- - L_3^+)^2 - 6q(L_2^- + L_2^+) + q^2(L_1^- - L_1^+)]}{12q^4} \right. \\ \left. - \frac{b_1(q_0 - 2k)[2(L_2^- - L_2^+) - q(L_1 + L_1^+)]}{2q^2} - \left(\frac{b_1}{3} + b_2\right)(L_1^- - L_1^+ + q) \right\}. \qquad (2.71) $$
There is no IR divergence in the s-channel: would-be singular terms contain inverse powers of $q$, but the integration domain extends to small $q$ only around $q_0 = 2k$, where the integrand vanishes for all statistics ($q_{\pm} = k + O(q)$).

### 2.5.3 IR divergence

Let us collect together the IR divergence affecting the $2 \leftrightarrow 2$ computation. Comparing eqs. (2.36)–(2.38) with eqs. (2.48) and (2.63) we can extract the coefficients appearing in eqs. (2.61) and (2.62):

$$ C\Phi_{g(g)}: \qquad a_1|_{+;++} = -4, \qquad b_1|_{+;++} = -4, \qquad (2.72) $$

$$ C[\Phi_{s(g)} + \Phi_{g(s)} + \Phi_{s|g}]: \qquad a_1|_{+;++} = -1, \qquad b_1|_{+;++} = -1, \qquad (2.73) $$

$$ C[\Phi_{s(f)} + \Phi_{f(s)} + \Phi_{s|f}]: \qquad a_2|_{-;-} = 4, \qquad b_2|_{+-} = 2, \qquad (2.74) $$

$$ C[\Phi_{f(g)} + \Phi_{g(f)} + \Phi_{f|g}]: \qquad a_1|_{-;-} = 4, \qquad b_1|_{+-} = 4. \qquad (2.75) $$

The coefficient $a_1$ only comes with the statistical factors that were considered in eq. (2.61), so that the IR divergence shown in eq. (2.62) is absent. Adding prefactors according to eq. (2.8) yields the total IR divergence of the $2 \leftrightarrow 2$ contribution:

$$ \begin{aligned} & \lim_{D \to 4} \mathrm{Im}\{G_{12;12}^{\mathrm{IR}}\}_{2 \leftrightarrow 2} \\ &= \frac{1}{(4\pi)^3 k} \int_{-\infty}^{k} dq_0 \int_{|q_0|}^{2k-q_0} dq [1 + n_B(q_0) + n_B(k-q_0)] \frac{6q_0(q^2-q_0^2)k^2\pi^2 T^2}{q^4} \\ & \times \left\{ g_1^2 \left[ \frac{n_S}{6} + \frac{n_G(N_c+1)(N_c+2)}{12N_c} \right] + 3g_2^2 \left[ \frac{2}{3} + \frac{n_S}{6} + \frac{n_G(N_c+1)}{12} \right] \right. \\ & \left. + (N_c^2-1)g_3^2 \left( \frac{N_c}{3} + \frac{n_G}{3} \right) \right\}. \end{aligned} \qquad (2.76) $$

## 2.6 Hard Thermal Loop resummation

The logarithmic IR divergence in eq. (2.76) can be eliminated through Hard Thermal Loop resummation [20, 21]. More precisely, as shown in ref. [17] for a fermionic production rate and in ref. [2] for the present observable, the infrared divergence is shielded through the so-called Landau damping part of a resummed propagator, corresponding physically to soft t-channel exchange.<fnref n="5" /> Thermal scatterings give an effective mass to the exchanged gauge boson, whereby the logarithmic divergence turns into a finite logarithm, as we show in the remainder of this section. In principle there could be a similar contribution from soft t-channel fermion exchange, however in practice there is no divergence at leading order, as we demonstrate in appendix A. Scalar fields do not experience Landau damping, so no discussion is needed for them. In the notation of eq. (2.8), we thus need to evaluate

$$ G_{12;12}^{\mathrm{E}}|_{\mathrm{HTL}} = \frac{2}{D(D-2)(D-3)} \left\{ 2n_{\mathrm{G}}(1 + N_{\mathrm{c}})\Phi_{f}|_{\mathrm{HTL}} + (2 + N_{\mathrm{c}}C_{\mathrm{F}})\Phi_{g}|_{\mathrm{HTL}} \right\}. \qquad (2.77) $$

<footnote n="5">⁵Originally this was shown in the context of photon production in QCD [22–25].</footnote>
Computing the diagram associated with $\Phi_g$ in figure 1 with HTL-resummed propagators, the result reads$^6$

$$ \Phi_g|_{\text{HTL}} = \frac{(D-2)L_{\mu\nu;\alpha\beta}}{2} \oint_Q 4\Theta_{\alpha\beta;\rho\sigma}(Q, K-Q) \Delta_{\sigma\lambda}^{\text{HTL}}(K-Q) \Theta_{\mu\nu;\lambda\kappa}(K-Q, Q) \Delta_{\kappa\rho}^{\text{HTL}}(Q), \quad (2.78) $$

where $\Delta^{\text{HTL}}$ is the gauge propagator,

$$ \Delta_{\mu\nu}^{\text{HTL}}(K) = \frac{\mathbb{P}_{\mu\nu}^{\text{T}}}{K^2 + \Pi_{\text{T}}(K)} + \frac{\mathbb{P}_{\mu\nu}^{\text{E}}}{K^2 + \Pi_{\text{E}}(K)} + \frac{\xi K_{\mu} K_{\nu}}{K^4}, \quad (2.79) $$

with $\mathbb{P}^\text{T}$ being the projector defined in eq. (2.4), $\xi$ a gauge parameter, and

$$ \mathbb{P}_{\mu\nu}^{\text{E}} = \delta_{\mu\nu} - \frac{K_{\mu}K_{\nu}}{K^4} - \mathbb{P}_{\mu\nu}^{\text{T}}. \quad (2.80) $$

The tensor $\Theta$ parametrizes the cubic graviton-gauge vertex,

$$ \Theta_{\alpha\beta;\rho\sigma}(P, Q) \equiv (P_{\alpha}\delta_{\mu\rho} - P_{\mu}\delta_{\alpha\rho})(Q_{\beta}\delta_{\mu\sigma} - Q_{\mu}\delta_{\beta\sigma}). \quad (2.81) $$

The full HTL computation can be simplified by noting that in the diagrams of figure 2, one of the gauge bosons attaching to the graviton vertex is always “hard” (i.e. with an external momentum $q \sim \pi T$) and only one is “soft” (i.e. an internal t-channel rung).$^7$ Adding to this that $\Theta$ projects out the longitudinal part of the propagator to which it is attached, permits us to replace $\Delta_{\sigma\lambda}^{\text{HTL}}(K-Q) \to 2\delta_{\sigma\lambda}/(K-Q)^2$, where the factor 2 accounts for the two possibilities of picking the hard line. Subsequently, after carrying out the contractions, we get

$$ \Phi_g|_{\text{HTL}} \approx \oint_Q \frac{4}{(K-Q)^2} \left\{ \left( \frac{1}{Q^2 + \Pi_T} - \frac{1}{Q^2 + \Pi_E} \right) \left[ (D-3)[\mathbb{P}_{\mathbf{q}}^{\text{T}}]^2 \left(D-2-\frac{Q^2}{q^2}+\frac{Dk^2}{2q^2}\right) \right. \right. \\ \left. - \frac{D(D-3)Q^2 \mathbb{P}_{\mathbf{q}}^{\text{T}}}{2} \left( \frac{\mathbf{q} \cdot \mathbf{k}}{q^2} + \frac{Q^2}{4q^2} \right) + \frac{D(D-2)(D-3)Q^4}{8} \right] \right. \\ \left. + \frac{1}{Q^2 + \Pi_E} \left[ (D-3)(D-2)[\mathbb{P}_{\mathbf{q}}^{\text{T}}]^2 - \frac{D(D-3)Q^2 \mathbb{P}_{\mathbf{q}}^{\text{T}}}{2} + \frac{D(D-2)(D-3)Q^4}{8} \right] \right\}. \quad (2.82) $$

Furthermore, we may focus on the contribution that is largest in the IR domain $q, q_0 \ll k$. This arises from the highest power of $k$ in the numerator, i.e. the term proportional to $k^2$ on the first line of eq. (2.82):

$$ \Phi_g|_{\text{HTL}}^{\text{IR}} = \oint_Q \frac{2D(D-3)}{(K-Q)^2} \left( \frac{1}{Q^2 + \Pi_T} - \frac{1}{Q^2 + \Pi_E} \right) \frac{k^2 [\mathbb{P}_{\mathbf{q}}^{\text{T}}]^2}{q^2}. \quad (2.83) $$

At this point we write the Euclidean propagators in a spectral representation,

$$ \frac{1}{Q^2 + \Pi(Q)} = \int_{-\infty}^{\infty} \frac{dq_0}{\pi} \frac{\rho(q_0, q)}{q_0 - iq_n}, \quad \rho(q_0, q) \equiv \text{Im} \left\{ \frac{1}{Q^2 + \Pi(Q)} \right\}_{q_n \to -i[q_0 + i0^+]}, \quad (2.84) $$

<footnote>$^6$The structure is the same for all three gauge groups, so we consider one of them as a representative.</footnote>

<footnote>$^7$This is also the reason for why vertices do not need to be resummed.</footnote>
carry out the Matsubara sum over $q_n$, and take the cut,

$$ \begin{aligned} \Gamma_{\text{HTL}} &\equiv \operatorname{Im}\left\{\oint_Q \frac{1}{(K-Q)^2[Q^2 + \Pi(Q)]}\right\}_{k_n \to -i[k+i0^+]} \\ &= \int_{-\infty}^{\infty} dq_0 \int_{\mathbf{q}} \frac{\rho(q_0, q)}{2\epsilon_{qk}} \left\{\delta(q_0 - k - \epsilon_{qk}) [n_{\text{B}}(\epsilon_{qk}) - n_{\text{B}}(q_0)] \right. \\ &\qquad \left. + \delta(q_0 - k + \epsilon_{qk}) [1 + n_{\text{B}}(q_0) + n_{\text{B}}(\epsilon_{qk})] \right\}, \end{aligned} \quad (2.85) $$

where $\epsilon_{qk} \equiv |\mathbf{q} - \mathbf{k}|$. Focussing on the soft contribution from the domain $q, q_0 \ll k$, only the latter channel gets kinematically realized. Carrying out the angular integral, this contribution can be expressed as

$$ \Gamma_{\text{HTL}} \supset \frac{1}{8\pi^2 k} \int_{-\infty}^{k} dq_0 \int_{|q_0|}^{2k-q_0} dq q [1 + n_{\text{B}}(q_0) + n_{\text{B}}(k-q_0)] \rho(q_0, q) \Big|_{\mathbf{q}\cdot\mathbf{k}=\frac{q^2-q_0^2+2kq_0}{2}}. \quad (2.86) $$

Inserting now the full structure of eq. (2.83) into eq. (2.86), we get

$$ \operatorname{Im}\left\{\Phi_g^{|_{\text{HTL}}|^{\text{IR}}}\right\}_{k_n \to -i[k+i0^+]} \stackrel{D \to 4}{\rightleftharpoons} \frac{1}{8\pi^2 k} \int_{-\infty}^{k} dq_0 \int_{|q_0|}^{2k-q_0} dq q [1+n_{\text{B}}(q_0)+n_{\text{B}}(k-q_0)] \\ \times \left\{ \frac{8k^2[\mathbb{P}_{\mathbf{q}}^{\text{T}}]^2}{q^2} [\rho_{\text{T}}(q_0,q)-\rho_{\text{E}}(q_0,q)] \right\}_{\mathbf{q}\cdot\mathbf{k}=\frac{q^2-q_0^2+2kq_0}{2}}. \quad (2.87) $$

The angular constraint implies that

$$ q_{\perp}^2 \equiv \mathbb{P}_{\mathbf{q}}^{\text{T}} = (q^2 - q_0^2) \frac{(k - q_+)(k - q_-)}{k^2} \stackrel{q_{\pm} \ll k}{\approx} q^2 - q_0^2. \quad (2.88) $$

The last step is invoked in order to carry out the resummation only for the leading term in an expansion in $q_0, q$, i.e. in the regime where there is an actual IR-divergence.

We now apply eq. (2.87) combined with the insertion of eq. (2.88) in two different ways. The first is to “re-expand” the result in the form of a weak-coupling expansion. In other words, the HTL spectral functions are evaluated for large $q, q_0$, whereby they become

$$ \rho_{\text{T}} \to \frac{\pi m_{\text{E}}^2 q_0}{4q^3(q^2 - q_0^2)}, \quad \rho_{\text{E}} \to -\frac{\pi m_{\text{E}}^2 q_0}{2q^3(q^2 - q_0^2)}. \quad (2.89) $$

Here the Debye mass $m_{\text{E}}$ reads, in the case of the different gauge groups,

$$ m_{\text{E1}}^2 = g_1^2 T^2 \left[ \frac{n_S}{6} + \frac{n_G(N_c+1)(N_c+2)}{12N_c} \right], \quad (2.90) $$

$$ m_{\text{E2}}^2 = g_2^2 T^2 \left[ \frac{2}{3} + \frac{n_S}{6} + \frac{n_G(N_c+1)}{12} \right], \quad (2.91) $$

$$ m_{\text{E3}}^2 = g_3^2 T^2 \left( \frac{N_c}{3} + \frac{n_G}{3} \right). \quad (2.92) $$
In this way we find

$$ \mathrm{Im}\left\{\Phi_g\right\}_{\mathrm{HTL}}^{\mathrm{IR}}\right\}_{k_n \to -i[k+i0^+]} = \frac{1}{8\pi^2 k} \int_{-\infty}^{k} dq_0 \int_{|q_0|}^{2k-q_0} dq [1 + n_B(q_0) + n_B(k-q_0)] \Lambda(q_0, q) \times \frac{6\pi q_0 (q^2 - q_0^2) k^2 m_E^2}{q^4}. \quad (2.93) $$

Here a function $\Lambda$ has been introduced, with the property $\lim_{q_0, q \to 0} \Lambda = 1$. It can be chosen at will outside of the domain where the resummation is implemented, given that its effects cancel up to higher-order corrections (cf. the discussion below eq. (2.97)).

Adding the prefactor from eq. (2.77) and resolving the different gauge groups,

$$ \frac{2(2 + N_c C_F) m_E^2}{D(D-2)(D-3)} \xrightarrow{D \to 4} \frac{m_{E1}^2 + 3m_{E2}^2 + (N_c^2 - 1)m_{E3}^2}{8}, \quad (2.94) $$

we reproduce the IR divergence from eq. (2.76) in the domain where $\Lambda = 1$.

The second way is that we evaluate the HTL contribution as such. This could be computed numerically after inserting the full spectral functions $\rho_{T,E}$ into eq. (2.87), but through an opportune choice of the weighting function $\Lambda$ it can also be determined analytically, by making use of a sum rule [26, 27]. First, according to eq. (2.88), we can substitute $q^2 \approx q_0^2 + q_\perp^2$, and use then $q_\perp$ and $q_0$ as integration variables. Second, for $q_0 \ll T$, the Bose distribution $n_B(q_0) \approx T/q_0$ dominates over the terms $1 + n_B(k-q_0)$ that are of order unity. It is helpful to employ this simplification, which can be implemented by choosing $\Lambda = \Lambda^*$, where

$$ [1 + n_B(q_0) + n_B(k - q_0)] \Lambda^*(q_0, q) \equiv \frac{T}{q_0}. \quad (2.95) $$

We also note that the difference $\rho_T(q_0, \sqrt{q_0^2 + q_\perp^2}) - \rho_E(q_0, \sqrt{q_0^2 + q_\perp^2})$ decreases rapidly at large $|q_0|$, whereby the integration range over $q_0$ can be extended to positive infinity. Therefore

$$ \begin{aligned} \mathrm{Im}\left\{\Phi_g\right\}_{\mathrm{HTL}}^{\mathrm{IR}}\right\}_{k_n \to -i[k+i0^+]} & \stackrel{\Lambda=\Lambda^*}{\approx} \frac{kT}{\pi^2} \int_{-\infty}^{\infty} \frac{dq_0}{q_0} \int_0^{2k} dq_\perp q_\perp \frac{q_\perp^4 [\rho_T(q_0, q) - \rho_E(q_0, q)]}{q^2} \\ & \stackrel{[26, 27]}{=} \frac{kT}{\pi} \int_0^{2k} dq_\perp q_\perp^3 \left( \frac{1}{q_\perp^2} - \frac{1}{q_\perp^2 + m_E^2} \right) \\ & = \frac{kT m_E^2}{2\pi} \ln \left( 1 + \frac{4k^2}{m_E^2} \right). \end{aligned} \quad (2.96) $$

This logarithmically enhanced term corresponds to that determined in ref. [2].

The full contribution of HTL resummation can now be obtained by subtracting the term in eq. (2.93) and adding that in eq. (2.96),

$$ \Delta \mathrm{Im}\left\{\Phi_g\right\}_{\mathrm{HTL}} = \mathrm{Im}\left\{\Phi_g\right\}_{\mathrm{HTL}}^{\mathrm{IR}}\right\}_{\mathrm{full}} - \mathrm{Im}\left\{\Phi_g\right\}_{\mathrm{HTL}}^{\mathrm{IR}}\right\}_{\mathrm{expanded}}. \quad (2.97) $$

Given that for $q_0, q \gg m_E$ the full and expanded HTL spectral functions agree up to terms of $O(g^4)$, the influence of $\Lambda$ drops out in this difference, however the same choice needs to be made in both terms (we chose $\Lambda = \Lambda^*$). The gauge groups are resolved as in eq. (2.94). The subtraction term is evaluated together with eq. (2.60), rendering the latter IR finite.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left: examples of the interaction rate $\Gamma(k)$ from eq. (2.2) at a few representative temperatures, normalized to $T^3/m_{Pl}^2$. The interaction rate decreases in these units with temperature, because the most important running couplings become smaller. Right: the combination $m_{Pl}^2 k^3 \Gamma(k) n_B(k)/T^6$ that plays a role for the production rate of the energy density carried by gravitational radiation.</figcaption></figure>

### 3 Numerical results

Inserting the integrals from eqs. (2.60) and (2.71), with coefficients from eqs. (2.72)–(2.75), into eq. (2.8), and adding the resummation from eq. (2.97), we can determine the interaction rate $\Gamma(k)$ from eq. (2.2). For the running couplings and Debye masses appearing in these expressions, we use values specified in section 4 of ref. [28].

In figure 3, $\Gamma(k)$ is plotted both as $m_{Pl}^2 \Gamma(k)/T^3$ and in the combination appearing in the energy density production rate, $m_{Pl}^2 k^3 \Gamma(k) n_B(k)/T^6$, at $T \approx 10^3, 10^9, 10^{15}$ GeV. In the units chosen, the rates decrease slowly with the temperature, due to the running of $g_2^2$, $g_3^2$ and $h_t^2$.

We remark that $\Gamma(k)$ has a (barely visible) negative dip for $k/T \to 0$. In this region many of our approximations, taken under the assumption $k \sim \pi T$, fail. Most importantly, HTL resummation with one hard and one soft gauge boson in $\Phi_g$, as described in section 2.6, only works correctly for $k \gg m_E$.<fnref n="8" /> This is neither new nor specific to graviton production: previous calculations of gravitino [29–31], axion [32, 33] and axino [34] production saw the same issue. In fact, the negative dips were typically much larger (cf., e.g., figure 3 of ref. [34]). The reason for the difference can be traced back to the way in which HTL resummation was implemented in these works, following ref. [35]. Even if the method

<footnote n="8">⁸For $k \gg m_E$, we could actually replace the argument of the logarithm in eq. (2.96) with just $4k^2/m_E^2$, as the difference between these is parametrically of $O(g^4)$. For $k \ll m_E/2$, however, $\ln(1 + 4k^2/m_E^2)$ is small and positive, whereas $\ln(4k^2/m_E^2)$ is large and negative. That said, our result is formally incomplete for $k \lesssim m_E$, as is practically any available thermal production rate as of today, including that of photons from QCD.</footnote>
agrees with ours for $k \sim \pi T$ up to terms of $O(g^4)$, it differs for $k \sim m_E$, in ways related to the discussion in footnote 8. Remarkably, our implementation of HTL resummation avoids large negative dips without resorting to partial, gauge-dependent resummations of higher-order effects that were introduced in refs. [31] and [33] for gravitino and axion production, respectively. These calculations could be revisited with our method, by finding the appropriate coefficients $a_i$ and $b_i$ for eqs. (2.60) and (2.71), and taking over our implementation of HTL resummation.

## 4 Cosmological implications

As a final step we embed the production rate in an expanding cosmological background and compute

$$ \Delta N_{\text{eff}} = \frac{8}{7} \left(\frac{11}{4}\right)^{\frac{4}{3}} \frac{e_{\text{GW}}(T_0)}{e_{\gamma}(T_0)}, \qquad (4.1) $$

where the final temperature can be chosen as $T_0 \sim 0.01 \text{ MeV}$ and $e_\gamma \equiv \pi^2 T_0^4 / 15$ is the energy density carried by photons. The constraints originating from $N_{\text{eff}}$ are analogous in spirit to the constraints on $e_{\text{GW}}$ considered in refs. [36, 37] (see also [4]), and recently $N_{\text{eff}}$ itself was invoked in ref. [38]. The uncertainties of the Standard Model prediction of $N_{\text{eff}}$ continue to be discussed in the literature (cf., e.g., refs. [39–41] and references therein), being around $\Delta N_{\text{eff}} \sim 10^{-3}$, whereas the current experimental accuracy is $\Delta N_{\text{eff}} \sim 10^{-1}$ [42], which is expected to be reduced by an order of magnitude by future facilities [43]. We consider the uncertainty of the Standard Model prediction, $\Delta N_{\text{eff}} \sim 10^{-3}$, to set an interesting sensitivity goal for considerations concerning the gravitational background.

Denoting by $H \equiv \sqrt{8\pi e_{\text{SM}}/(3m_{\text{Pl}}^2)}$ the Hubble rate, by $s_{\text{SM}}$ the Standard Model entropy density, and by $c_s^2$ the speed of sound squared, the energy density at $T_0$ can be obtained as [2]

$$ \frac{e_{\text{GW}}(T_0)}{s_{\text{SM}}^{4/3}(T_0)} = \int_0^{\ln\left(\frac{T_{\text{max}}}{T_0}\right)} \frac{dx}{3c_s^2 H} \int_{\mathbf{k}} R(T, k), \quad x \equiv \ln\left(\frac{T_{\text{max}}}{T}\right), \qquad (4.2) $$

where the production rate $R$ is related to the damping coefficient $\Gamma$ from eq. (1.1) through

$$ R(T, k) \equiv 2k \Gamma(k) n_{\text{B}}(k). \qquad (4.3) $$

The integrand of eq. (4.2) is illustrated in figure 4(left) as a function of the temperature. Clearly the integral is dominated by the high-temperature end, so in practice we may restrict to temperatures above the electroweak crossover, $T \sim 160 \text{ GeV}$, for its determination.

The entropy dilution that takes place at low temperatures is accounted for by the factor $s_{\text{SM}}^{4/3}(T_0)$ in eq. (4.2). We have adopted a prescription for $s_{\text{SM}}$ which permits for its use even at $T < 2 \text{ MeV}$ when neutrinos have decoupled (cf. the web page associated with ref. [44] for the specification and for the numerical values that have been used$^9$).

Putting everything together, the contribution of the gravitational wave background to $N_{\text{eff}}$, obtained from eq. (4.1), is shown in figure 4(right). Once the experimental accuracy reaches the level $\Delta N_{\text{eff}} \approx 10^{-3}$, maximal temperatures above $2 \times 10^{17} \text{ GeV}$ can be constrained.

<footnote>$^9$The numerical values are attached to this publication as supplementary material.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Left: the integrated production rate of the energy density carried by gravitational radiation, normalized as in eq. (4.2), as a function of the temperature. Only the high-temperature end plays a significant role. Right: the contribution of the gravitational energy density to the parameter $N_{\text{eff}}$ (cf. eq. (4.1)), as a function of the highest temperature of the radiation epoch. Once the experimental determination of $N_{\text{eff}}$ reaches the current theoretical precision, $\Delta N_{\text{eff}} \sim 10^{-3}$, reheating temperatures above $T_{\text{max}} \approx 2 \times 10^{17}$ GeV can be constrained.</figcaption></figure>

## 5 Conclusions and outlook

The main purpose of this paper has been to refine the estimate $T_{\text{max}} \lesssim 10^{17}...18$ GeV that was obtained for the maximal temperature of the radiation epoch in ref. [2], by promoting the previous leading-logarithmic analysis to a full leading-order computation of the energy density carried by gravitational radiation emitted by a Standard Model plasma. If the experimental determination of the parameter $N_{\text{eff}}$ can reach the current theoretical accuracy, $\Delta N_{\text{eff}} \sim 10^{-3}$, and no deviations from the Standard Model prediction are found, the refined estimate reads $T_{\text{max}} \le 2 \times 10^{17}$ GeV. It is remarkable that this model-independent constraint is not much weaker than typical bounds on the reheating temperature that are obtained by comparing model-dependent inflationary predictions with Planck data [42].

Most of the energy density carried by thermally produced gravitational radiation peaks in the microwave frequency range today. Conceivably, this physics can be probed by tabletop experiments in the future [45–53], even if the sensitivity goal is quite formidable.

With future extensions in mind, we have displayed the technical steps of the computation in quite some detail (cf. section 2). The partly automatized procedure to determine the matrix elements squared in eqs. (2.43)–(2.46) can be straightforwardly extended to other models. The IR subtraction and thermal resummation that were described in section 2.6 must still be adjusted accordingly, however we hope that our exposition lays out these steps in a digestible fashion. Apart from graviton production in Beyond the Standard Model theories, this machinery can be applied to the production rates of other particles coupling to a heat bath via non-renormalizable operators, such as gravitinos (with $M \ll \pi T$), axions
and axinos. Indeed, as mentioned in section 3, the phase space integration and resummation prescriptions of sections 2.5, 2.6, which do not suffer from large, unphysical negative contributions at small $k/T$, can be directly applied to the known matrix elements squared in the literature [29, 32, 34].

## Acknowledgments

This work was partly supported by the Swiss National Science Foundation (SNF) under grant 200020B-188712.

## A Soft t-channel fermion exchange

We analyze in this appendix the fermion exchange part of eq. (2.77), viz. $\Phi_f|_{HTL}$, and show that no resummation is needed at leading order.

Computing the diagram associated with $\Phi_f$ in figure 1 within the HTL theory, the result reads<fnref n="10" />

$$ \Phi_f|_{HTL} = \frac{(D-2)L_{\mu\nu;\alpha\beta}}{2} \oint_{\{Q\}} \mathrm{Tr} \{ \Upsilon_{\mu\nu}(Q, K+Q) G^{\mathrm{HTL}}(K+Q) \Upsilon_{\alpha\beta}(K+Q, Q) G^{\mathrm{HTL}}(Q) \}, \quad (A.1) $$

where $G^{\mathrm{HTL}}$ is the HTL-resummed fermion propagator,

$$ G^{\mathrm{HTL}}(K) = \frac{ik_n \gamma_0}{K^2 + \Pi_W(K)} + \frac{ik_i \gamma_i}{K^2 + \Pi_P(K)}, \quad (A.2) $$

and the tensor $\Upsilon$ parametrizes the cubic graviton-fermion vertex,

$$ \Upsilon_{\alpha\beta}(P, Q) = \frac{\gamma_{\alpha}(P_{\beta} + Q_{\beta}) + \gamma_{\beta}(P_{\alpha} + Q_{\alpha})}{4}. \quad (A.3) $$

Like in the gluonic case, we can replace one of the propagators by a free one ($\Pi_{W,P} \to 0$ in eq. (A.2)) and account for the associated symmetry by a factor 2. Taking the Dirac trace, this leads to

$$ \Phi_f|_{HTL} \approx \oint_{\{Q\}} \frac{2(D-3)}{(K+Q)^2} \left\{ \frac{1}{Q^2 + \Pi_W} [-D\mathbb{P}_{\mathbf{q}}^{\mathrm{T}}(q_n^2 + q_n k_n)] + \frac{1}{Q^2 + \Pi_P} [4[\mathbb{P}_{\mathbf{q}}^{\mathrm{T}}]^2 - D\mathbb{P}_{\mathbf{q}}^{\mathrm{T}}(q^2 + \mathbf{q} \cdot \mathbf{k})] \right\}. \quad (A.4) $$

Writing now

$$ q_n^2 + q_n k_n = -(q^2 + \mathbf{q} \cdot \mathbf{k}) + \frac{(K+Q)^2 + Q^2 - K^2}{2}, \quad (A.5) $$

and noting that $K^2$ vanishes on the light cone after analytic continuation and that $(K+Q)^2$ gives no cut as it cancels the free propagator, we can identify the most IR sensitive terms as those proportional to $\mathbf{q} \cdot \mathbf{k}$.

<footnote n="10">¹⁰The structure is the same for all fermions, so we consider one Dirac-like fermion as a representative.</footnote>
Next, we invoke a spectral representation like in eq. (2.84), carry out the Matsubara sum over $q_n$, and take the cut,

$$ \begin{aligned} \tilde{\Gamma}_{\text{HTL}} &\equiv \operatorname{Im}\left\{\oint_{\{Q\}} \frac{1}{(K+Q)^2[Q^2 + \Pi(Q)]}\right\}_{k_n \to -i[k+i0^+]} \\ &= -\int_{-\infty}^{\infty} \mathrm{d}q_0 \int_{\mathbf{q}} \frac{\rho(q_0, q)}{2\tilde{\epsilon}_{qk}} \left\{\delta(q_0 + k - \tilde{\epsilon}_{qk})[n_F(q_0) - n_F(\tilde{\epsilon}_{qk})] \right. \\ &\qquad \left. + \delta(q_0 + k + \tilde{\epsilon}_{qk})[1 - n_F(q_0) - n_F(\tilde{\epsilon}_{qk})]\right\}, \end{aligned} \quad (A.6) $$

where $\tilde{\epsilon}_{qk} \equiv |\mathbf{q} + \mathbf{k}|$. Focussing on the soft contribution from the domain $q, q_0 \ll k$, only the first channel gives a contribution. It is convenient to substitute $q_0 \to -q_0$ and make use of the antisymmetry $\rho(-q_0, q) = -\rho(q_0, q)$. Carrying out the angular integral, this yields

$$ \tilde{\Gamma}_{\text{HTL}} \supset \frac{1}{8\pi^2 k} \int_{-\infty}^{k} \mathrm{d}q_0 \int_{|q_0|}^{2k-q_0} \mathrm{d}q \, q [1-n_F(q_0)-n_F(k-q_0)] \rho(q_0, q) \Big|_{\mathbf{q}\cdot\mathbf{k}=\frac{q_0^2-q^2-2kq_0}{2}}. \quad (A.7) $$

We note from the angular constraint in eq. (A.7) that for the most IR sensitive contribution we can replace $\mathbf{q} \cdot \mathbf{k} \to -kq_0$. Combining this with eqs. (A.4) and (A.5) leads us to focus on

$$ \operatorname{Im}\left\{\Phi_f\right\}_{\text{HTL}}^{\text{IR}} \Big|_{k_n \to -i[k+i0^+]} \stackrel{D \to 4}{\equiv} \frac{1}{8\pi^2} \int_{-\infty}^{k} \mathrm{d}q_0 \int_{|q_0|}^{2k-q_0} \mathrm{d}q \, q [1 - n_F(q_0) - n_F(k-q_0)] \times 8q_0 \mathbb{P}_{\mathbf{q}}^{\mathrm{T}} [\rho_{\mathrm{P}}(q_0, q) - \rho_{\mathrm{W}}(q_0, q)], \quad (A.8) $$

where $\mathbb{P}_{\mathbf{q}}^{\mathrm{T}}$ can be taken over from eq. (2.88).

Again, we evaluate eq. (A.8) in two ways. Re-expanding in a strict weak-coupling expansion, the spectral functions become

$$ \rho_{\mathrm{P}} \to -\frac{\pi m_{\mathrm{A}}^2 q_0}{4q^3(q^2 - q_0^2)}, \quad \rho_{\mathrm{W}} \to -\frac{\pi m_{\mathrm{A}}^2}{4qq_0(q^2 - q_0^2)}. \quad (A.9) $$

Here $m_{\mathrm{A}}$ is a so-called asymptotic thermal mass [54], which for quarks reads

$$ m_{q_L}^2 = \frac{(g_1^2 Y^2 + 3g_2^2/4 + g_3^2 C_F)T^2}{4}, \quad m_{u_R, d_R}^2 = \frac{(g_1^2 Y^2 + g_3^2 C_F)T^2}{4}, \quad (A.10) $$

where Y denotes the hypercharge assignment as listed below eq. (2.7). For the leptons, the SU(3) parts are absent. Inserting eq. (A.9) into eq. (A.8) yields

$$ \operatorname{Im}\left\{\Phi_f\right\}_{\text{HTL}}^{\text{IR}} \Big|_{k_n \to -i[k+i0^+]} = \frac{1}{8\pi^2} \int_{-\infty}^{k} \mathrm{d}q_0 \int_{|q_0|}^{2k-q_0} \mathrm{d}q \, [1 - n_F(q_0) - n_F(k-q_0)] \times \frac{2\pi(q^2 - q_0^2)m_{\mathrm{A}}^2}{q^2}. \quad (A.11) $$

This is integrable (i.e. IR finite) at $q, q_0 \ll k$, and therefore does not appear in eq. (2.76).
<figure><img src="image_6.png" /><figcaption>Figure 5. An example of a 2 ↔ 3 scattering contributing to gravitational wave production. The notation is as in figure 2, and the magnitude of these scatterings is estimated in appendix B.</figcaption></figure>

A complementary view on the soft fermion contribution can be obtained by evaluating eq. (A.8) like we did for the gauge contribution in eq. (2.96). Making use of a sum rule derived in ref. [17], and making a choice analogous to eq. (2.95), this gives

$$
\begin{align}
\mathrm{Im}\left\{\Phi_f^\mathrm{IR}|_{\mathrm{HTL}}\right\}_{k_n \to -i[k+i0^+]}^\mathrm{full} & \approx \frac{1}{\pi^2} \int_{-\infty}^{\infty} \mathrm{d}q_0 q_0 \int_0^{2k} \mathrm{d}q_{\perp} q_{\perp} \left[\frac{1}{2} - n_{\mathrm{F}}(k)\right] q_{\perp}^2 [\rho_{\mathrm{P}}(q_0, q) - \rho_{\mathrm{W}}(q_0, q)] \nonumber \\
& \stackrel{[17]}{=} \frac{1}{2\pi} \left[\frac{1}{2} - n_{\mathrm{F}}(k)\right] \int_0^{2k} \mathrm{d}q_{\perp} q_{\perp}^3 \frac{m_{\mathrm{A}}^2}{q_{\perp}^2 + m_{\mathrm{A}}^2}. \tag{A.12}
\end{align}
$$

The integral is dominated by $q_{\perp} \sim 2k$, yielding a contribution of $O(g^2T^4)$ for $k \sim \pi T$. This is of leading order, but just a part of the full result, not justifying any resummation.

All in all, soft fermion exchange does not need to be resummed at leading order.

## B Magnitude of $1+n \leftrightarrow 2+n$ processes

The processes we have considered in the main text, illustrated in figure 2, correspond to $2 \leftrightarrow 2$ scatterings. It may be asked if $1+n \leftrightarrow 2+n$ reactions also contribute. As Standard Model particles obtain thermal masses, whereas gravitons remain massless, there is no phase space for such a process at the Born level ($n=0$). However, if one of the particles interacts before emitting a gravitational wave ($n \ge 1$), so that it is set slightly off-shell, this argument no longer applies. An example of this type of a “bremsstrahlung” process is shown in figure 5. In the context of producing photons or massless fermions from a thermal plasma, such processes do contribute at the same order as $2 \leftrightarrow 2$ scatterings, and have to be summed to all orders ($\sum_{n=0}^{\infty}$), through a procedure known as Landau-Pomeranchuk-Migdal (LPM) resummation [55–57]. In footnote 1 of ref. [33], it has however been pointed out that such reactions are of subleading order for gravitational wave production. The purpose of this appendix is to confirm the assertion of ref. [33], which we do by employing light-cone variables similar to those normally adopted for LPM resummation.

In the notation of eq. (2.8), i.e. treating the gauge groups on equal footing for a moment, the LPM contribution reads

$$
G_{12;12}^{\mathrm{E}}|_{\mathrm{LPM}} = \frac{2}{D(D-2)(D-3)} \left\{ n_S \Phi_s|_{\mathrm{LPM}} + 2n_G(1 + N_c)\Phi_f|_{\mathrm{LPM}} + (2 + N_c C_F)\Phi_g|_{\mathrm{LPM}} \right\}. \qquad (B.1)
$$
In order to determine the three terms, we start by writing their (vanishing) Born limits in a suggestive form.<fnref n="11" /> According to eqs. (2.11)-(2.13), the cuts read

$$ \lim_{D \to 4} \mathrm{Im}\{\Phi_s\}|_{k_n \to -i[\omega+i0^+]} = 4 \mathrm{Im}\{J_{11}^2\}|_{k_n \to -i[\omega+i0^+]}, \qquad (B.2) $$

$$ \lim_{D \to 4} \mathrm{Im}\{\Phi_f\}|_{k_n \to -i[\omega+i0^+]} = -2 \mathrm{Im}\{2\tilde{J}_{11}^2 + \tilde{J}_{11}^1\}|_{k_n \to -i[\omega+i0^+]}, \qquad (B.3) $$

$$ \lim_{D \to 4} \mathrm{Im}\{\Phi_g\}|_{k_n \to -i[\omega+i0^+]} = 2 \mathrm{Im}\{2J_{11}^2 + 4J_{11}^1 + J_{11}^0\}|_{k_n \to -i[\omega+i0^+]}, \qquad (B.4) $$

where the masters $J, \tilde{J}$ were defined in eq. (2.9) and we have kept $\omega \neq k$. Let us approach the light cone from above, setting $\omega \equiv \sqrt{k^2 + M^2}$ with $M^2 \to 0^+$. Adopting results from eqs. (2.85) and (A.6) and setting $Q \to -Q$ in the latter, we can write

$$ \begin{aligned} \Gamma_{\mathrm{LPM}} &\equiv \mathrm{Im} \left\{ \oint_Q \frac{\alpha_0 [\mathbb{P}_q^\mathrm{T}]^2 + \alpha_1 \mathbb{P}_q^\mathrm{T} K^2 + \alpha_2 K^4]}{(K-Q)^2 Q^2} \right\}_{k_n \to -i[\omega+i0^+]} \\ &= \int_{-\infty}^{\infty} \mathrm{d}q_0 \int_{\mathbf{q}} \frac{\rho_{\mathrm{free}}(q_0, q)}{2\epsilon_{qk}} \{\alpha_0 [\mathbb{P}_q^\mathrm{T}]^2 - \alpha_1 \mathbb{P}_q^\mathrm{T} M^2 + \alpha_2 M^4\} \\ &\qquad \times \{\delta(q_0 - \omega - \epsilon_{qk}) [n_\sigma(\epsilon_{qk}) - n_\sigma(q_0)] + \delta(q_0 - \omega + \epsilon_{qk}) [1 + n_\sigma(q_0) + n_\sigma(\epsilon_{qk})]\}, \end{aligned} \qquad (B.5) $$

where $\sigma = \pm$ takes care of statistics according to eq. (2.27). The free spectral function reads

$$ \rho_{\mathrm{free}}(q_0, q) = \frac{\pi[\delta(q_0 - q) - \delta(q_0 + q)]}{2q}. \qquad (B.6) $$

For $M^2 > 0$ the contribution comes from the second kinematic channel in eq. (B.5) combined with the first term in eq. (B.6).

We now go over to light-cone coordinates, $\mathbf{q} = q_\parallel \mathbf{e}_k + \mathbf{q}_\perp$, so that

$$ \epsilon_{qk} = \sqrt{(k - q_\parallel)^2 + q_\perp^2}, \quad \mathbb{P}_{\mathbf{q}}^\mathrm{T} = q_\perp^2. \qquad (B.7) $$

The constraint $\delta(q_0 - q)$ is eliminated by integrating over $q_\parallel$, which sets $q_\parallel = \sqrt{q_0^2 - q_\perp^2}$ (here we anticipate the overall sign to be positive, $q_\parallel \sim q_0 \in (0, k)$, cf. eq. (B.9)). The remaining constraint $\delta(q_0 - \omega + \epsilon_{qk})$ implies

$$ M^2 = \omega^2 - k^2 = 2 \left[ q_0^2 - k \sqrt{q_0^2 - q_\perp^2} + q_0 \sqrt{k^2 + q_0^2 - 2k \sqrt{q_0^2 - q_\perp^2}} \right]. \qquad (B.8) $$

This can be expanded in $q_\perp^2/q_0^2$ and $q_\perp^2/(k-q_0)^2$, assuming again $0 < q_0 < k$ to fix signs. Keeping contributions up to $q_\perp^4$ in $\alpha_0 q_\perp^4 - \alpha_1 q_\perp^2 M^2 + \alpha_2 M^4$ and contributions up to $q_\perp^2$ inside $\delta(q_0 - \omega + \epsilon_{qk})$, we find

$$ \lim_{D \to 4} \mathrm{Im}\{\Phi_i\}|_{k_n \to -i[\omega+i0^+]} \stackrel{M^2 \approx 0}{\sim} \int_{-\infty}^{\infty} \mathrm{d}q_0 \kappa_i(q_0) [1 + n_\sigma(q_0) + n_\sigma(k-q_0)] \\ \times \int_{\mathbf{q}_\perp} q_\perp^4 \delta\left(-\frac{M^2}{2k} + \frac{q_\perp^2}{2(k-q_0)} + \frac{q_\perp^2}{2q_0}\right). \qquad (B.9) $$

<footnote n="11">¹¹More precisely, to extract the information, all terms contributing to the “slope” towards the vanishing limit need to be included, which in the current context amounts to terms $\propto K^4$.</footnote>
It is clear from here that for $M^2 > 0$ the contribution originates from $0 < q_0 < k$. However, we have removed the specifier $M^2 \to 0^+$, because eq. (B.9) turns out to be applicable for $M^2 \to 0^-$ as well, with the contribution originating from $q_0 < 0$ and $q_0 > k$ in that case.

When the coefficients $\alpha_0, \alpha_1, \alpha_2$ are inserted into the prefactor according to eqs. (B.2)–(B.4), the functions $\kappa_i$ in eq. (B.9) become

$$ \kappa_s(q_0) = \frac{1}{2q_0(k-q_0)}, \quad \kappa_f(q_0) = \frac{q_0^2 + (k-q_0)^2}{4q_0^2(k-q_0)^2}, \quad \kappa_g(q_0) = \frac{q_0^4 + (k-q_0)^4}{4q_0^3(k-q_0)^3}. \quad (B.10) $$

Up to overall conventions, $\kappa_s$ and $\kappa_f$ agree with the prefactors cited for scalars and fermions in ref. [56]. The factor $\kappa_g$ is similar to the prefactor for the gluon contribution to gluon emission that was discussed in ref. [58], however it is not exactly the same: the latter has an additional $k^4$ in the numerator, guaranteeing a symmetry between the three gluons involved.

Let us now estimate the magnitude of the $1+n \leftrightarrow 2+n$ contributions. For this, we can set the virtuality to be parametrically $M^2 \sim g^2 T^2$, as it is at this scale that thermal masses and scatterings of the type in figure 5 play a role if $q_0 \sim k \sim \pi T$. Then eq. (B.9) implies that $q_\perp^2 = M^2 q_0 (k-q_0)/k^2$ and, up to logarithms in the case of $\kappa_g$, $\mathrm{Im}\{\Phi_i\} \sim M^4 \sim g^4 T^4$. This is suppressed by $\mathcal{O}(g^2)$ compared with the effects that we are interested in.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] S. Weinberg, *Gravitation and Cosmology* John Wiley & Sons, New York (1971).

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

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

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

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

[6] G. Baym, S.P. Patil and C.J. Pethick, *Damping of gravitational waves by matter*, Phys. Rev. D 96 (2017) 084033 [arXiv:1707.05192] [INSPIRE].

[7] R. Flauger and S. Weinberg, *Absorption of gravitational waves from distant sources*, Phys. Rev. D 99 (2019) 123030 [arXiv:1906.04853] [INSPIRE].

[8] A. Abbas, *Anomalies and charge quantization in the standard model with arbitrary number of colours*, Phys. Lett. B 238 (1990) 344 [INSPIRE].

[9] Y. Schröder, M. Vepsäläinen, A. Vuorinen and Y. Zhu, *The ultraviolet limit and sum rule for the shear correlator in hot Yang-Mills theory*, JHEP 12 (2011) 035 [arXiv:1109.6548] [INSPIRE].
[10] 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].

[11] A. Alloul, N.D. Christensen, C. Degrande, C. Duhr and B. Fuks, *FeynRules 2.0 — A complete toolbox for tree-level phenomenology*, Comput. Phys. Commun. **185** (2014) 2250 [arXiv:1310.1921] [INSPIRE].

[12] N.D. Christensen et al., *A comprehensive approach to new physics simulations*, Eur. Phys. J. C **71** (2011) 1541 [arXiv:0906.2474] [INSPIRE].

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

[14] T. Hahn and M. Pérez-Victoria, *Automated one-loop calculations in four and D dimensions*, Comput. Phys. Commun. **118** (1999) 153 [hep-ph/9807565] [INSPIRE].

[15] B.R. Holstein, *Graviton physics*, Am. J. Phys. **74** (2006) 1002 [gr-qc/0607045] [INSPIRE].

[16] 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].

[17] D. Besak and D. Bödeker, *Thermal production of ultrarelativistic right-handed neutrinos: complete leading-order results*, JCAP **03** (2012) 029 [arXiv:1202.1288] [INSPIRE].

[18] M. Laine, *Thermal 2-loop master spectral function at finite momentum*, JHEP **05** (2013) 083 [arXiv:1304.0202] [INSPIRE].

[19] G. Jackson, *Two-loop thermal spectral functions with general kinematics*, Phys. Rev. D **100** (2019) 116019 [arXiv:1910.07552] [INSPIRE].

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

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

[22] J.I. Kapusta, P. Lichard and D. Seibert, *High-energy photons from quark-gluon plasma versus hot hadronic gas*, Phys. Rev. D **44** (1991) 2774 [Erratum ibid. **47** (1993) 4171] [INSPIRE].

[23] R. Baier, H. Nakkagawa, A. Niégawa and K. Redlich, *Production rate of hard thermal photons and screening of quark mass singularity*, Z. Phys. C **53** (1992) 433 [INSPIRE].

[24] T. Altherr and P.V. Ruuskanen, *Low-mass dileptons at high momenta in ultra-relativistic heavy-ion collisions*, Nucl. Phys. B **380** (1992) 377 [INSPIRE].

[25] P.B. Arnold, G.D. Moore and L.G. Yaffe, *Photon emission from quark gluon plasma: complete leading order results*, JHEP **12** (2001) 009 [hep-ph/0111107] [INSPIRE].

[26] P. Aurenche, F. Gelis and H. Zaraket, *A Simple sum rule for the thermal gluon spectral function and applications*, JHEP **05** (2002) 043 [hep-ph/0204146] [INSPIRE].

[27] S. Caron-Huot, *$O(g)$ plasma effects in jet quenching*, Phys. Rev. D **79** (2009) 065039 [arXiv:0811.1603] [INSPIRE].

[28] M. Laine, P. Schicho and Y. Schröder, *A QCD Debye mass in a broad temperature range*, Phys. Rev. D **101** (2020) 023532 [arXiv:1911.09123] [INSPIRE].

[29] M. Bolz, A. Brandenburg and W. Buchmüller, *Thermal production of gravitinos*, Nucl. Phys. B **606** (2001) 518 [Erratum ibid. **790** (2008) 336] [hep-ph/0012052] [INSPIRE].
[30] J. Pradler and F.D. Steffen, Thermal gravitino production and collider tests of leptogenesis, Phys. Rev. D 75 (2007) 023509 [hep-ph/0608344] [INSPIRE].

[31] V.S. Rychkov and A. Strumia, Thermal production of gravitinos, Phys. Rev. D 75 (2007) 075011 [hep-ph/0701104] [INSPIRE].

[32] P. Graf and F.D. Steffen, Thermal axion production in the primordial quark-gluon plasma, Phys. Rev. D 83 (2011) 075011 [arXiv:1008.4528] [INSPIRE].

[33] A. Salvio, A. Strumia and W. Xue, Thermal axion production, JCAP 01 (2014) 011 [arXiv:1310.6982] [INSPIRE].

[34] A. Brandenburg and F.D. Steffen, Axino dark matter from thermal production, JCAP 08 (2004) 008 [hep-ph/0405158] [INSPIRE].

[35] E. Braaten and T.C. Yuan, Calculation of screening in a hot plasma, Phys. Rev. Lett. 66 (1991) 2183 [INSPIRE].

[36] T.L. Smith, E. Pierpaoli and M. Kamionkowski, New Cosmic Microwave Background Constraint to Primordial Gravitational Waves, Phys. Rev. Lett. 97 (2006) 021301 [astro-ph/0603144] [INSPIRE].

[37] S. Henrot-Versillé et al., Improved constraint on the primordial gravitational-wave density using recent cosmological data and its impact on cosmic string models, Class. Quant. Grav. 32 (2015) 045003 [arXiv:1408.5299] [INSPIRE].

[38] B.X. Hu and A. Loeb, An Upper Limit on the Initial Temperature of the Radiation-Dominated Universe, arXiv:2004.02895 [INSPIRE].

[39] P.F. de Salas and S. Pastor, Relic neutrino decoupling with flavour oscillations revisited, JCAP 07 (2016) 051 [arXiv:1606.06986] [INSPIRE].

[40] J.J. Bennett, G. Buldgen, M. Drewes and Y.Y.Y. Wong, Towards a precision calculation of the effective number of neutrinos $N_{\text{eff}}$ in the Standard Model: the QED equation of state, JCAP 03 (2020) 003 [arXiv:1911.04504] [INSPIRE].

[41] M. Escudero Abenza, Precision early universe thermodynamics made simple: $N_{\text{eff}}$ and neutrino decoupling in the Standard Model and beyond, JCAP 05 (2020) 048 [arXiv:2001.04466] [INSPIRE].

[42] PLANCK collaboration, Planck 2018 results. VI. Cosmological parameters, arXiv:1807.06209 [INSPIRE].

[43] K. Abazajian et al., CMB-S4 Science Case, Reference Design and Project Plan, arXiv:1907.04473 [INSPIRE].

[44] M. Laine and M. Meyer, Standard Model thermodynamics across the electroweak crossover, JCAP 07 (2015) 035 [arXiv:1503.04935] [INSPIRE] [http://www.laine.itp.unibe.ch/eos15/].

[45] A.M. Cruise and R.M.J. Ingley, A prototype gravitational wave detector for 100 MHz, Class. Quant. Grav. 23 (2006) 6185 [INSPIRE].

[46] T. Akutsu et al., Search for a Stochastic Background of 100-MHz Gravitational Waves with Laser Interferometers, Phys. Rev. Lett. 101 (2008) 101101 [arXiv:0803.4094] [INSPIRE].

[47] F. Li, J. Baker, Robert M. L., Z. Fang, G.V. Stephenson and Z. Chen, Perturbative photon fluxes generated by high-frequency gravitational waves and their physical effects, Eur. Phys. J. C 56 (2008) 407 [arXiv:0806.1989] [INSPIRE].
[48] M.-l. Tong, Y. Zhang and F.-Y. Li, *Using a polarized maser to detect high-frequency relic gravitational waves*, Phys. Rev. D **78** (2008) 024041 [arXiv:0807.0885] [INSPIRE].

[49] C. Sabín, D.E. Bruschi, M. Ahmadi and I. Fuentes, *Phonon creation by gravitational waves*, *New J. Phys.* **16** (2014) 085003 [arXiv:1402.7009] [INSPIRE].

[50] M. Goryachev and M.E. Tobar, *Gravitational wave detection with high frequency phonon trapping acoustic cavities*, Phys. Rev. D **90** (2014) 102005 [arXiv:1410.2334] [INSPIRE].

[51] D. Hartley, T. Bravo, D. Rätzel, R. Howl and I. Fuentes, *Analogue simulation of gravitational waves in a 3+1 dimensional Bose-Einstein condensate*, Phys. Rev. D **98** (2018) 025011 [arXiv:1712.01140] [INSPIRE].

[52] R. Schützhold, *Interaction of a Bose-Einstein condensate with a gravitational wave*, Phys. Rev. D **98** (2018) 105019 [arXiv:1807.07046] [INSPIRE].

[53] M.P.G. Robbins, N. Afshordi and R.B. Mann, *Bose-Einstein condensates as gravitational wave detectors*, JCAP **07** (2019) 032 [arXiv:1811.04468] [INSPIRE].

[54] H. Weldon, *Effective fermion masses of order gT in high-temperature gauge theories with exact chiral invariance*, Phys. Rev. D **26** (1982) 2789 [INSPIRE].

[55] P. Aurenche, F. Gelis and H. Zaraket, *Landau-Pomeranchuk-Migdal effect in thermal field theory*, Phys. Rev. D **62** (2000) 096012 [hep-ph/0003326] [INSPIRE].

[56] P.B. Arnold, G.D. Moore and L.G. Yaffe, *Photon emission from ultrarelativistic plasmas*, JHEP **11** (2001) 057 [hep-ph/0109064] [INSPIRE].

[57] A. Anisimov, D. Besak and D. Bödeker, *Thermal production of relativistic Majorana neutrinos: strong enhancement by multiple soft scattering*, JCAP **03** (2011) 042 [arXiv:1012.3784] [INSPIRE].

[58] P.B. Arnold, G.D. Moore and L.G. Yaffe, *Photon and gluon emission in relativistic plasmas*, JHEP **06** (2002) 030 [hep-ph/0204343] [INSPIRE].