Gravitational wave background from Standard Model physics: complete leading order

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 ∼ π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 Neff parameter. The current theoretical uncertainty ∆Neff≤ 10−3 corresponds to Tmax≤ 2 × 1017 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.


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 2 Pl and therefore tiny for temperatures much below the Planck scale. On the other hand, the age of the universe (inverse Hubble rate) is ∼ 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.

JHEP07(2020)092
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] f GW (t, k) = Γ(k) n B (k) − f GW (t, k) + O 1 m 4 Pl , (1.1) where k ≡ |k| and n B (k) ≡ 1/(e k/T − 1) is the Bose distribution. The differential energy density is given by de GW = 2k f GW d 3 k (2π) 3 . Adopting a logarithmic scale, the production rate of gravitational energy density can thus be expressed as (1.2) In the following we are interested in estimating the rate Γ(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 ∼ π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 n B (k), so that the right-hand side of eq. (1.1) evaluates to Γ(k)n B (k). However, the same coefficient Γ(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.

Setup
Assuming that a system is spatially homogeneous and stationary on the time scales observed, and aligning the z-axis with the momentum (k = k e z ), the production rate of the energy density carried by gravitational waves can be related to the Wightman correlator G < 12;12 ≡ X e ik(t−z) T 12 (0) T 12 (X ) , X ≡ (t, x) . (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 µ − u i u µ )(. . .) µ .
In equilibrium, G < 12;12 is related to the imaginary part of the retarded correlator as G < 12;12 = 2n B (k) Im G R 12;12 . In the following we compute a Euclidean correlator G E 12;12 as a function of a Euclidean four-momentum K = (k n , k), from which G R 12;12 is obtained by JHEP07(2020)092 an analytic continuation, G R 12;12 = G E 12;12 | k n →−i[k+i0 + ] . Here k n = 2πnT , with n ∈ Z, is a bosonic Matsubara frequency. The rate Γ(k) from eq. (1.1) is then given by [2] Γ(k) = 16π Im G R where D denotes the dimension of space-time, X ≡ (τ, x), and τ ∈ (0, 1 T ). Here we have defined the projector (L µν;αβ L αβ;γδ = L µν;γδ ) which is symmetric (L µν;αβ = L νµ;αβ = L αβ;µν ) and projects onto transverse (K µ L µν;αβ = k i δ iµ L µν;αβ = 0) and traceless (δ µν L µν;αβ = 0) modes. We also denote As T µν we take the Standard Model energy-momentum tensor, which we write in Euclidean metric. Given that L µν;αβ projects out trace parts, it is enough to include non-trace ones, where the a i label the generators of the various gauge groups; φ is the Higgs doublet; q L , L are the left-handed quark and lepton doublets, respectively; and u R , d R , ν R , e R are the corresponding right-handed components. The covariant derivative has the form where g 1 , g 2 , g 3 are gauge couplings, a L is the left-handed projector and the hypercharge assignments [8]. We note that because of their vanishing gauge charge assignments and the omission of their Yukawa couplings, the fields ν R do not contribute to 2 ↔ 2 scatterings and have thus no effect on our final results (traditionally, ν R are often omitted from the outset). 1 A simple way to verify the factor in the denominator is to consider momentum averages in the transverse plane. By rotational symmetry, q i q j q k q l = A (δ ij δ kl +δ ik δ jl +δ il δ jk ). Therefore a representative of T 12 T 12 evaluates to q 2 1 q 2 2 = A, whereas L ij;kl q i q j q k q l = A D(D − 3).

JHEP07(2020)092
In order to avoid inverse polynomials of D in section 2.2, the result for G E 12;12 is expressed as where n S = 1 is the number of Higgs doublets, n G ≡ 3 is the number of fermion generations, , and O(g 4 ) refers generically to any 3-loop contribution. 2 Here s, f, g refer to effects from scalars, fermions, and gauge bosons, respectively; Φ a is a 1-loop diagram with a particle of type a; Φ a(b) is a 2-loop diagram where a particle of type a couples to T µν and a particle of type b appears in a loop; and Φ 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.

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] where {P } denotes a fermionic Matsubara four-momentum. The indices 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 ↔ 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]. 2 The Higgs self-coupling and top Yukawa coupling appear in a Euclidean Lagrangian as After the use of symmetries related to substitutions of integration variables, and noting that terms with odd numbers of γ 5 -matrices do not contribute at this order, the results read  The computation was carried out in a general covariant gauge, and we have checked that the gauge parameter drops out exactly. The result for Φ g(g) can be crosschecked against ref. [9].

Extracting 2 ↔ 2 cuts at light cone
As discussed below eq. (2.1), from each Φ we need to extract the cut Im Φ| k n →−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 We denote the phase space of 2 ↔ 2 scatterings by . Distribution functions are denoted by so that n + = n B and n − = −n F are the Bose and Fermi distributions, respectively. Distribution functions appear in the combination Mandelstam variables are defined as usual, s ≡ (

JHEP07(2020)092
With this notation, the 2 ↔ 2 cuts for the structures in eq. (2.25) read where σ a , σ c and σ e label the statistics of the 1 st , 3 rd and 5 th 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 → 4 for simplicity, as there are no ultraviolet divergences in these cuts. Denoting by C an operation which produces an integrand for eq. (2.29), viz.
and making use of symmetries such asĨ fgh 1b101 =Ī hgf 1b101 (obtained by the substitution P → Q − P ), the non-zero contributions for the combinations appearing in eq. (2.8) read 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 · p 1 − kp 1 ), and recalling that P T can make use of energy-momentum conservation to verify that With this identity, combined with renamings p 1 ↔ p 2 as well as a repeated use of s+t+u = 0, all projectors 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 N τ 1 ;σ 1 σ 2 : 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, Im G R T , can be replaced by the full vector correlator, Im G R V = Im G R T + Im G R L , because a Ward identity guarantees the vanishing of Im G R L for zero virtuality. We are not aware of a similar operator relation between the tensor channel correlator in JHEP07(2020)092 eq. (2.3) and one without any P T 's, even if intriguing relations between photon and graviton production amplitudes are known to exist (cf. section 2.4).

Connection to Boltzmann equations
The 2 ↔ 2 cuts of section 2.3 can also be obtained from kinetic theory and Boltzmann equations. As a starting point, we may, for k ∼ πT , write the leading-order contribution to eq. (1.1) aṡ (2.39) where we have neglected f GW (t, k) on the right-hand side. The sum runs over all abc ∈ SM (Standard Model) particle and antiparticle degrees of freedom and thus over all ab → cG processes, with G denoting the graviton. |M ab cG (p 1 , p 2 ; k 1 , 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 ∼ π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 = b this factor is compensated for by their being counted twice in the sum over abc. The thermal distributions f i correspond to n B and n F for bosons and fermions, respectively, with [1 ± f c (k 1 )] implying [1 + n B (k 1 )] in the former case and [1 − n 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.
where the SM energy-momentum tensor T µν SM 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 assigments were needed). This package and its companion FormCalc [14] were then used to generate, evaluate and square all amplitudes, summing over the relevant degeneracies. 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 JHEP07(2020)092 proceeded as follows. FeynArts assigns to external tensor bosons a polarization tensor λ µν (k) which is written, using a common factorization formula (cf., e.g., refs. [15,16]), as with λ µ (k) the transverse polarization vector of a massless gauge boson. Upon taking k = k e z and the circular polarization vectors λ µ (k) = 1/ it is easy to verify that the polarization sum satisfies 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 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 To verify the agreement, relabellings p 1 ↔ p 2 (and t ↔ u) as well as use of the identity N τ 1 ; 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 γ 5 matrices in Dirac traces, since they vanish under the dΩ 2→2 integration. Specifically, these terms appear in the f g → f G 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 2 3 -part of eq. (2.43) corresponds to the matrix elements squared for the gluonic scattering gg → gG, yielding Recently, there has been much work on factorizing graviton amplitudes into photon amplitudes multiplied by kinematic factors, say f γ → f G versus f γ → f γ (cf., e.g.,
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 ∼ α 2 s T , where elementary particle states need to be replaced by hydrodynamic modes [2].

Phase space integrals
The next step is to carry out the phase space integral dΩ 2→2 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]. 4 We discuss separately the treatment of t and s-channel cases (u-channel can always be transformed into t-channel).

t-channel
Consider the phase space integral The idea is to insert 1 = d 4 Q δ (4) (P 1 − K 1 − Q) in the integral. Then the energymomentum conservation constraint inside dΩ 2→2 can be written as δ (4) (Q + P 2 − K). We can now integrate over p 2 and k 1 by using the spatial parts of the Dirac δ's, leaving q 0 , q and p 1 as the integration variables. The temporal Dirac δ's fix two angles as whereas kinematic variables become The azimuthal average of powers of k · p 1 can be computed by parametrizing q = (0, 0, q) , k = k (sin χ, 0, cos χ) , p 1 = p 1 (sin θ cos ϕ, sin θ sin ϕ, cos θ) , (2.53)

JHEP07(2020)092
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 thereby factorizing the p 1 -dependence. Denoting the integration range of p 1 can be established as (q + , ∞). The integration measure contains no powers of p 1 , whereas azimuthal averages yield powers up to p 2 1 . The integral reads where All in all this results in The integral in eq. (2.60) is logarithmically IR divergent at small q 0 , q. For the different statistics the divergent parts read (2.62) JHEP07(2020)092

s-channel
The s-channel phase space integral is defined as This time we insert 1 = d 4 Q δ (4) (P 1 + P 2 − Q) in the integral, whereby the energymomentum conservation constraint inside dΩ 2→2 can be written as δ (4) (Q − K 1 − K). We integrate over p 1 and k 1 by using the spatial parts of the Dirac δ's, leaving q 0 , q and p 2 as the integration variables. The temporal Dirac δ's fix two angles as whereas kinematic variables become The azimuthal average of powers of k · p 2 can be computed like in eqs. (2.52)-(2.53), exchanging p 1 ↔ p 2 . The phase space distributions from eq. (2.28) are now cast in the form 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 q + q − dp 2 β 0 + β 1 p 2 + β 2 p 2 2 1 + n σ 1 (q 0 − p 2 ) + n σ 2 (p 2 ) All in all, this gives

JHEP07(2020)092
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 ± = k + O(q)).

IR divergence
Let us collect together the IR divergence affecting the 2 ↔ 2 computation. Comparing eqs.
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 ↔ 2 contribution: (2.76)

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. 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 (2.77) 5 Originally this was shown in the context of photon production in QCD [22][23][24][25].

JHEP07(2020)092
Computing the diagram associated with Φ g in figure 1 with HTL-resummed propagators, the result reads 6 (2.78) where ∆ HTL is the gauge propagator, with P T being the projector defined in eq. (2.4), ξ a gauge parameter, and The tensor Θ parametrizes the cubic graviton-gauge vertex, 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 ∼ πT ) and only one is "soft" (i.e. an internal t-channel rung). 7 Adding to this that Θ projects out the longitudinal part of the propagator to which it is attached, permits us to replace ∆ HTL σλ (K − Q) → 2δ σλ /(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 (2.82) Furthermore, we may focus on the contribution that is largest in the IR domain q, q 0 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): At this point we write the Euclidean propagators in a spectral representation,

JHEP07(2020)092
carry out the Matsubara sum over q n , and take the cut, where qk ≡ |q − k|. Focussing on the soft contribution from the domain q, q 0 k, only the latter channel gets kinematically realized. Carrying out the angular integral, this contribution can be expressed as Inserting now the full structure of eq. (2.83) into eq. (2.86), we get The angular constraint implies that 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 .

(2.89)
Here the Debye mass m E reads, in the case of the different gauge groups, (2.92)

JHEP07(2020)092
In this way we find Here a function Λ has been introduced, with the property lim q 0 ,q→0 Λ = 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, we reproduce the IR divergence from eq. (2.76) in the domain where Λ = 1.
The second way is that we evaluate the HTL contribution as such. This could be computed numerically after inserting the full spectral functions ρ T,E into eq. (2.87), but through an opportune choice of the weighting function Λ 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 ≈ q 2 0 + q 2 ⊥ , and use then q ⊥ and q 0 as integration variables. Second, for q 0 T , the Bose distribution n B (q 0 ) ≈ 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 Λ = Λ , where We also note that the difference ρ T (q 0 , q 2 0 + q 2 ⊥ ) − ρ E (q 0 , q 2 0 + q 2 ⊥ ) decreases rapidly at large |q 0 |, whereby the integration range over q 0 can be extended to positive infinity. Therefore 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), (2.97) Given that for q 0 , q m E the full and expanded HTL spectral functions agree up to terms of O(g 4 ), the influence of Λ drops out in this difference, however the same choice needs to be made in both terms (we chose Λ = Λ ). 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. 2) at a few representative temperatures, normalized to T 3 /m 2 Pl . The interaction rate decreases in these units with temperature, because the most important running couplings become smaller. Right: the combination m 2 Pl k 3 Γ(k) n B (k)/T 6 that plays a role for the production rate of the energy density carried by gravitational radiation.
In figure 3, Γ(k) is plotted both as m 2 Pl Γ(k)/T 3 and in the combination appearing in the energy density production rate, m 2 Pl k 3 Γ(k) n B (k)/T 6 , at T ≈ 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 2 3 and h 2 t . We remark that Γ(k) has a (barely visible) negative dip for k/T → 0. In this region many of our approximations, taken under the assumption k ∼ πT , fail. Most importantly, HTL resummation with one hard and one soft gauge boson in Φ g , as described in section 2.6, only works correctly for k m E . 8 This is neither new nor specific to graviton production: previous calculations of gravitino [29][30][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 8 For k m E , we could actually replace the argument of the logarithm in eq. (2.96) with just 4k 2 /m 2 E , as the difference between these is parametrically of O(g 4 ). For k m E /2, however, ln(1 + 4k 2 /m 2 E ) is small and positive, whereas ln(4k 2 /m 2 E ) is large and negative. That said, our result is formally incomplete for k < ∼ m E , as is practically any available thermal production rate as of today, including that of photons from QCD.

JHEP07(2020)092
agrees with ours for k ∼ πT up to terms of O(g 4 ), it differs for k ∼ 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 higherorder 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.

Cosmological implications
As a final step we embed the production rate in an expanding cosmological background and compute where the final temperature can be chosen as T 0 ∼ 0.01 MeV and e γ ≡ π 2 T 4 0 /15 is the energy density carried by photons. The constraints originating from N eff are analogous in spirit to the constraints on e GW considered in refs. [36,37] (see also [4]), and recently N eff itself was invoked in ref. [38]. The uncertainties of the Standard Model prediction of N eff continue to be discussed in the literature (cf., e.g., refs. [39][40][41] and references therein), being around ∆N eff ∼ 10 −3 , whereas the current experimental accuracy is ∆N eff ∼ 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, ∆N eff ∼ 10 −3 , to set an interesting sensitivity goal for considerations concerning the gravitational background.
Denoting by H ≡ 8πe SM /(3m 2 Pl ) the Hubble rate, by s SM the Standard Model entropy density, and by c 2 s the speed of sound squared, the energy density at T 0 can be obtained as [2] e GW (T 0 ) where the production rate R is related to the damping coefficient Γ from eq. (1.1) through 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 ∼ 160 GeV, for its determination. The entropy dilution that takes place at low temperatures is accounted for by the factor s 4/3 SM (T 0 ) in eq. (4.2). We have adopted a prescription for s SM which permits for its use even at T < 2 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 eff , obtained from eq. (4.1), is shown in figure 4(right). Once the experimental accuracy reaches the level ∆N eff ≈ 10 −3 , maximal temperatures above 2 × 10 17 GeV can be constrained. 9 The numerical values are attached to this publication as supplementary material. JHEP07(2020)092  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 eff (cf. eq. (4.1)), as a function of the highest temperature of the radiation epoch. Once the experimental determination of N eff reaches the current theoretical precision, ∆N eff ∼ 10 −3 , reheating temperatures above T max ≈ 2 × 10 17 GeV can be constrained.

Conclusions and outlook
The main purpose of this paper has been to refine the estimate T max < ∼ 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 eff can reach the current theoretical accuracy, ∆N eff ∼ 10 −3 , and no deviations from the Standard Model prediction are found, the refined estimate reads T max ≤ 2 × 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].
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 πT ), axions JHEP07(2020)092 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. Φ f HTL , and show that no resummation is needed at leading order.
Computing the diagram associated with Φ f in figure 1 within the HTL theory, the result reads 10 (A.1) where G HTL is the HTL-resummed fermion propagator, and the tensor Υ parametrizes the cubic graviton-fermion vertex, Like in the gluonic case, we can replace one of the propagators by a free one (Π W,P → 0 in eq. (A.2)) and account for the associated symmetry by a factor 2. Taking the Dirac trace, this leads to Writing now 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 q · k.

JHEP07(2020)092
Next, we invoke a spectral representation like in eq. (2.84), carry out the Matsubara sum over q n , and take the cut, where˜ qk ≡ |q + k|. Focussing on the soft contribution from the domain q, q 0 k, only the first channel gives a contribution. It is convenient to substitute q 0 → −q 0 and make use of the antisymmetry ρ(−q 0 , q) = −ρ(q 0 , q). Carrying out the angular integral, this yields We note from the angular constraint in eq. (A.7) that for the most IR sensitive contribution we can replace q · k → −kq 0 . Combining this with eqs. (A.4) and (A.5) leads us to focus on where P T q 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 Here m A is a so-called asymptotic thermal mass [54], which for quarks reads 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 This is integrable (i.e. IR finite) at q, q 0 k, and therefore does not appear in eq. (2.76).
JHEP07(2020)092 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.
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 The integral is dominated by q ⊥ ∼ 2k, yielding a contribution of O(g 2 T 4 ) for k ∼ π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 ↔ 2 + n processes The processes we have considered in the main text, illustrated in figure 2, correspond to 2 ↔ 2 scatterings. It may be asked if 1 + n ↔ 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 ≥ 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 ↔ 2 scatterings, and have to be summed to all orders ( ∞ n=0 ), through a procedure known as Landau-Pomeranchuk-Migdal (LPM) resummation [55][56][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

JHEP07(2020)092
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 → 0 + , because eq. (B.9) turns out to be applicable for M 2 → 0 − as well, with the contribution originating from q 0 < 0 and q 0 > k in that case. When the coefficients α 0 , α 1 , α 2 are inserted into the prefactor according to eqs. (B.2)-(B.4), the functions κ i in eq. (B.9) become κ s (q 0 ) = 1 2q 0 (k − q 0 ) , κ f (q 0 ) = q 2 0 + (k − q 0 ) 2 4q 2 0 (k − q 0 ) 2 , κ g (q 0 ) = (B.10) Up to overall conventions, κ s and κ f agree with the prefactors cited for scalars and fermions in ref. [56]. The factor κ 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 ↔ 2+n contributions. For this, we can set the virtuality to be parametrically M 2 ∼ 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 ∼ k ∼ πT . Then eq. (B.9) implies that q 2 ⊥ = M 2 q 0 (k − q 0 )/k 2 and, up to logarithms in the case of κ g , Im Φ i ∼ M 4 ∼ g 4 T 4 . This is suppressed by 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.