Double-graviton production from Standard Model plasma

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


Introduction
The past years have seen an explosion of interest in gravitational wave astronomy.With hopes for an upcoming observation of a primordial background through pulsar timing arrays [1][2][3][4], the construction of the LISA observatory under way, and the planning of DECIGO as well as the Einstein Telescope striding forward, theoretical efforts at identifying all conceivable cosmological sources are well motivated.In fact such efforts should cover not only the sub-kHz frequencies relevant for large interferometers, but extend up to much higher frequencies, 100 GHz and above, for which future experimental concepts are being actively developed [5].
If we restrict ourselves to the Standard Model of particle physics, perhaps minimally extended to incorporate the existence of neutrino masses, there is however only one primordial gravitational wave component that is guaranteed to be present. 1It is the gravitational waves produced by microscopic collisions [11] taking place in a Standard Model plasma [12,13].This background is not a perfect analogue of the cosmic microwave background (CMB), since gravitons were most likely never in thermal equilibrium, 2 however its shape (∼ Planckian) and peak frequency range (f 0 | peak ∼ 100GHz) are similar to those of the CMB.The reason for the similarity is that the graviton energies reflect the thermal energy distributions of the plasma particles from which they are produced.
A key property of the thermal gravitational wave background is that while its shape can be computed, its amplitude is not known.The reason is that gravitons do not equilibrate; they are just continuously being produced.The production rate is proportional to the strength of the gravitational interaction, ∼ 1/m 2 pl , where m pl is the Planck mass, as well as to the strength of Standard Model forces, ∼ α, where α is a fine-structure constant.Integrating over a Hubble time, H −1 ∼ m pl /T 2 , and normalizing to the energy density of radiation, the 1 Many other sources have been analyzed, notably inflation [6], the dynamics of reheating [7,8], and a plethora of post-reheating phenomena [9,10], however they rely on Beyond the Standard Model (BSM) physics. 2 Exotic scenarios in which this could have been the case have also been proposed, cf., e.g., refs.[14,15].
fractional energy density in gravitational waves then scales as Ω gw ∼ αT max /m pl , where T max is the maximal temperature of the radiation-dominated phase in our universe.Now, if we explore large values of T max , it must be asked whether terms suppressed by higher powers of T max /m pl also play a role.A concrete way to probe this question was identified in ref. [16].The rate ∼ α/m 2 pl mentioned above involves the production of a single graviton in association with a Standard Model particle.But in addition, there are processes in which two gravitons are produced, without any Standard Model particles.The rate of the latter processes is ∼ 1/m 4 pl .To complete the dimensions, this rate must come with a power of T max .Therefore, at high T max , the double-graviton rate overtakes the single-graviton rate.Then we also need to ask whether the notion of using General Relativity as an effective theory, reliant on the convergence of an expansion in T 2 max /m 2 pl , continues to apply.The purpose of the current study is to promote the computation of ref. [16], carried out in scalar field theory, to the full Standard Model.This requires including fermion-antifermion pairs and gauge boson pairs as additional initial states.Summing over all processes, we provide a quantitative criterion for when single-graviton production dominates.
Our presentation is organized as follows.The main steps of the computation, as well as a difference to matrix elements squared that can be found in the literature, are discussed in sec.2. Our results are illustrated numerically in sec.3, and the main conclusions are formulated in sec.4. Some technical details related to evaluating thermal averages of pair production cross sections have been relegated to appendix A.

Steps of the computation
We consider a moment after inflation at which the plasma filling the universe has equilibrated to an average temperature T ≪ m pl .Furthermore, the plasma contribution to the overall energy density is assumed to dominate over that from non-equilibrated fields, such as the inflaton.Then the Hubble rate, H, is much smaller than typical thermal momenta, H ∼ √ g * T 2 /m pl ≪ πT , where g * ∼ 10 2 is the number of relativistic degrees of freedom.
Consequently, thermal wavelengths are well within the horizon, (πT ) −1 ≪ H −1 .We consider the production of gravitational waves (GW's) from the scatterings of particles with such wavelengths.As will become apparent, the GW spectrum reflects that of the scatterers, so that in terms of the physical momentum k, it is peaked at around k| peak ∼ πT .
Given that the physics we are interested in takes place within the horizon, the space is locally flat, and we can employ local Minkowskian coordinates for the computation.The background metric is written as g µν = η µν +κ h µν , where η µν = diag(+−−−) is the Minkowski metric and h µν denotes the gravitational perturbation.The effective coupling of gravitational  interactions is defined as where G is Newton's constant, and m pl = 1.22091 × 10 19 GeV is the Planck mass.
When the Einstein-Hilbert action is expanded as a series in κ, vertices are generated for the interactions of h µν with itself and with other fields (a pedagogic review can be found in ref. [17]).Adding a graviton leg to a process leads to a suppression of the amplitude by κ ∼ 1/m pl , and to a suppression of the corresponding amplitude squared by κ 2 .But we can also add a normal Standard Model vertex to any given process, which leads to the suppression of the rate by ∼ α, where α is a fine-structure constant.Therefore, the production rate can be viewed as a double series, in α and 1/m 2 pl .In addition, we may recall that 2 → 1 processes are kinematically forbidden for a massless final state.If we approximate all particles as massless, the leading processes are then 2 → 2 amplitudes.If they involve a Standard Model vertex, their rate is ∼ α/m 2 pl ; if they only contain gravitational vertices, their rate is ∼ 1/m 4 pl .To these leading terms, additional vertices of either type can be added, leading to further suppression by α or 1/m 2 pl (e.g.triple-graviton production at ∼ 1/m 6 pl , etc).
The practical computations can be streamlined by going over to the De Donder gauge, This also renders the quadratic part of the gravitational action invertible, so that an internal graviton propagator can be found.On the external graviton lines, which are on-shell, i.e. with four-momentum K = (k, k) with k ≡ |k|, the sum over polarizations λ yields the projector λ={×,+} which is traceless in αβ and µν, transverse with respect to the graviton four-momentum K, and projects onto two physical states in four dimensions, i.e.Ä αβ; αβ = 2.The tensor Ã T can be chosen as which is transverse with respect to the four-momenta K and K ≡ (k, −k), and also with respect to the three-momentum k (we work in the plasma rest frame).
With the Feynman rules at hand, the first step is to determine the matrix elements squared for the production processes, shown in fig. 1.The four-momenta of the Standard Model particles are denoted by P 1 , P 2 and those of the gravitons by The processes take place at very high temperatures, much above any known particles masses.Therefore all particles can effectively be taken to be massless, whereby For massless fermions, sums over spins yield the standard expressions τ =± u(τ, p)ū(τ, p) = τ =± v(τ, p)v(τ, p) = / P .For massless gauge bosons, a polarization sum yields where È T µν is defined like in eq.(2.3), but with four-momentum P rather than K.It projects onto two physical states in four dimensions, i.e. |È Tµ µ | = 2. Inserting the Feynman rules to the amplitudes in fig. 1, with the cubic graviton vertex being particularly cumbersome [17]; constructing the amplitudes squared; and contracting with polarization or spin sums, as specified above, leads to rather lengthy expressions.We have found it practical to manipulate them with FORM [18], FeynArts [19], FormCalc [20], and/or FeynCalc [21].In the course of these computations, we have verified the gauge independence of the contribution originating from massless gauge bosons, notably that the last term of eq.(2.3) (with K → P), containing P, does not contribute.
The matrix elements squared are best tabulated in a form where we sum not only over all spins and polarization states, both of initial and final-state particles, but also over particles and antiparticles.For a complex scalar field φ, the processes from fig. 1(a) then yield [16] all In the fermionic case, we take a Dirac fermion of definite chirality (i.e. a Weyl fermion) as a building block, with two physical degrees of freedom.However, as we sum over fermions and antifermions, we effectively have four degrees of freedom for each fermionic external leg, and therefore each of them can effectively be represented as a usual Dirac fermion ψ.The processes in fig. 1 Finally, for a U(1) gauge boson g, the processes in fig.1(c) produce We note that, up to overall normalization, eqs.(2.6) and (2.7) can be extracted as the massless limits of differential cross sections given in ref. [22].However, this is not the case for eq.(2.8), as the assignment of a mass to vector fields is ambiguous and breaks gauge invariance.If masses are given to gauge fields, then the physical origin of the longitudinal polarization states should be specified (e.g. through the Higgs mechanism), in order to obtain a gauge independent expression.Given the matrix elements squared, we can determine the thermally averaged production rate of gravitational radiation.As a tool for this, we introduce the polarization-averaged phase-space distribution of gravitons, f gw .The energy density carried by gravitational waves of physical momentum k can now be expressed as where the factor 2 counts the physical polarization states.In a locally Minkowskian time, making use of the rotational invariance of the thermal ensemble, the production rate can subsequently be expressed as (2.10) The time variation of the phase-space distribution can be extracted from a Boltzmann equation.Assuming the particle content of the Standard Model (with three generations of charged leptons ℓ, neutrinos ν, and up and down-type quarks u, d, each with left and right chiralities L and R), this takes at O(κ 4 ) the form ḟgw where p i ≡ |p i |, N c = 3, and the Bose and Fermi distributions have been introduced as f B (ǫ) ≡ 1/(e ǫ/T − 1) and f F (ǫ) ≡ 1/(e ǫ/T + 1), respectively.The prefactor 1/(8k) is a combination of the standard 1/(2k) associated with phase space, a factor 1/2 for cancelling the sum over graviton polarizations that we had included in |M| 2 , as well as 1/2 for cancelling the overcounting of identical particles or the redundant sum over particles and antiparticles introduced above.Furthermore, we note that eq.(2.11) includes only the gain terms, no loss terms, because the graviton phase-space density is much below the equilibrium value, f gw (k) ≪ f B (k).For the same reason, there are no Bose enhancement factors for the final-state gravitons in the gain terms.
Let us remark that right-handed neutrinos have been included as degrees of freedom in eq.(2.11).However, they only interact via Yukawa couplings, whose magnitudes are unknown.If they are small, right-handed neutrinos might not equilibrate fast enough to be part of the thermal ensemble, and should be omitted as initial states.We will illustrate the numerical difference between including and not including right-handed neutrinos in fig. 2. It is perhaps appropriate to mention that the equilibration rates of many other particles have also been discussed in the literature, however this typically concerns particle asymmetry changing reactions relevant for chemical equilibration (cf., e.g., ref. [23]), whereas kinetic equilibrium should be efficiently maintained by gauge interactions.
The phase-space average in eq.(2.11) can be conveniently implemented by taking the schannel momentum transfer as the integration variable.For this, we write P i = (p i , p i ) and where integrals were carried out over p 1 and k 1 .The constraints fix two angles as cos θ q,p 2 = q 2 + q 0 (2p 2 − q 0 ) 2qp 2 , cos θ q,k = q 2 + q 0 (2k − q 0 ) 2qk . (2.13) The remaining angle can be expressed in terms of an azimuthal variable, ϕ, as cos θ k,p 2 = cos θ q,k cos θ q,p 2 + sin θ q,k sin θ q,p 2 cos ϕ . (2.14) Carrying out the angular integrals in order to remove the Dirac-δ's, we then obtain (2.15) Subsequently, the integrals over p 2 and ϕ can be performed analytically, as detailed in appendix A. We represent the result as which is rapidly convergent and readily evaluated numerically. 3

Numerical results
A numerical illustration of eq.(2.17) is given in fig. 2. We have plotted the rate in units of T 7 /m 2 pl , whereby it scales with the temperature as T 2 /m 2 pl .The reason for this choice is that the result can then be conveniently compared with the rate originating from single-graviton production [13], which is practically constant in these units.For both results we have assumed the Standard Model matter content.Additional particle species could be included in a fairly straightforward manner (cf., e.g., refs.[26][27][28][29][30][31]), but we do not expect them to change the overall pattern, even though they could have a quantitative influence.
The main observation from fig. 2 is that in the domain of the peak frequency, the doublegraviton production rate is below the single-graviton production rate at T < 4 × 10 18 GeV.In Left: the single (h) [13] and double (hh) graviton production rates from a Standard Model plasma.The slow evolution of the single-graviton rate with the temperature is due to a logarithmic running of the couplings, whereas the double-graviton rate scales with the temperature as T 2 /m 2 pl in the units chosen.At T = 10 18 GeV, the double-graviton rate is still more than an order of magnitude smaller than the single-graviton rate, but it would become dominant at T > 4 × 10 18 GeV.Right: different contributions to the double-graviton rate from eq. ( 2.17) at T = 10 18 GeV.In the fermionic part and the sum, the narrow band indicates the effect of omitting right-handed neutrinos.order to put this in context, we recall that the inflationary constraint H max < 10 −5 m pl [32], combined with the assumption of instantaneous reheating, suggests T max < 10 16 GeV.Therefore, within standard inflationary scenarios, the thermal part of the gravitational wave production rate can always be estimated from the single-graviton contribution.
If we ignore the standard inflationary paradigm, another constraint on gravitational waves originates from N eff , parametrizing the energy density residing in light degrees of freedom at the time of primordial nucleosynthesis (∆e gw,bbn ≡ ∆N eff (7/8)(4/11) 4/3 e γ,bbn ).In order to obtain the gravitational contribution, we integrate the evolution of the phase space distribution, appearing in eq.(2.10), over the cosmological history, along a comoving trajectory [12].The time integral can be converted into an integral over temperature, from the current T 0 up to a maximal one, T max , yielding where t 0 is the current time, a is the cosmological scale factor, and the relationship t ↔ T is to be determined from the Friedmann equations.Subsequently, normalizing with the current The time-integrated contribution of single and double-graviton production to the differential gravitational energy density today, h 2 dΩ gw,0 /d ln f 0 , as a function of the current-day frequency, f 0 , and various maximal temperatures, T max .We have assumed radiation-dominated expansion at T < T max [33], and denoted by h the reduced Hubble rate.The results are compared with a curve of the same shape (grey line), which after integration over f 0 yields h 2 Ω gw,0 ≡ 5.62 × 10 −6 ∆N eff (cf., e.g., ref. [9]).Here we have set ∆N eff = 0.2, which is the maximal observationally allowed contribution to the effective number of massless neutrino species [34].
Given that our result for h 2 Ω gw,0 depends on T max , we now get a relationship between T max and ∆N eff .This can be compared with existing theoretical and experimental knowledge about ∆N eff .In particular, the current estimated uncertainty of the Standard Model prediction, ∆N eff ≤ 0.001, corresponds to T max ≤ 2 × 10 17 GeV [13].At T max = 4 × 10 18 GeV, when the double-graviton rate overtakes the single-graviton rate, ∆N eff ≈ 0.02. 5 This is still below the current experimental limit, ∆N eff < 0.2 [34].For illustration, a spectrum that would yield ∆N eff = 0.2 is also shown in fig. 3.

Conclusions and outlook
In this work we have computed the double-graviton production rate from a thermal Standard Model plasma present during an early radiation-dominated epoch.Previous calculations had only considered scalar field contributions to the double-graviton rate [16].For momenta k ∼ πT , the double-graviton rate scales as de gw /[dt dln k] ∼ T 9 /m 4 pl .This can be compared with the single-graviton production rate [13], which is of order de gw /[dt dln k] ∼ α T 7 /m 2 pl , where α is a fine-structure constant.
Our result for the double-graviton production rate is given in eq.(2.17) and illustrated numerically in fig. 2. A comparison with the single-graviton production rate shows that the O(T 9 /m 4 pl ) contribution overtakes the O(α T 7 /m 2 pl ) one at T ∼ 4 × 10 18 GeV.This temperature is much higher than expected to be reached after standard inflationary scenarios.Therefore, single-graviton production is normally the dominant source.If one considers BSM scenarios associated with a non-renormalizable coupling, ∼ 1/Λ 2 (cf., e.g., ref. [27]), then we should effectively replace α → T 2 /Λ 2 in the single-graviton rate.In these cases, if Λ ∼ m pl , the single and double-graviton contributions could be of comparable magnitudes.
In order to obtain our results, we needed to evaluate the matrix elements squared corresponding to the Feynman diagrams in fig. 1, producing the results given in eqs.(2.6)-(2.8).Similar matrix elements squared have been considered in previous literature, either out of theoretical interest (cf., e.g., ref. [22]), or in view of the scattering of gravitational waves on matter (cf., e.g., ref. [35]).As we are concerned with very high temperatures, the Standard Model particles can be considered as massless for the purposes of our analysis.We noted that the literature results had assigned masses to gauge fields in a way which did not permit to recover the massless limit of eq.(2.8).
Let us elaborate briefly on the more general conceptual implications of our results.As explained below eq.(2.1), the gravitational wave production rate within General Relativity is a double series, in α and T 2 /m 2 pl .Normally, to estimate the convergence of the series in T 2 /m 2 pl , we would compare terms with a same power of α.However, the term of O(α 0 T 2 /m 2 pl ) is absent, because it is kinematically forbidden, and the term of O(α 1 T 4 /m 4 pl ) is unknown, as it is of next-to-leading order, either in α with respect to the known contribution of O(α 0 T 4 /m 4 pl ), or in T 2 /m 2 pl with respect to the known contribution of O(α 1 T 2 /m 2 pl ).But we can still compare the known terms.Because of the absence of α, the term of O(α 0 T 4 /m 4 pl ) overtakes the O(α 1 T 2 /m 2 pl ) contribution sooner than the proper probe of O(α 1 T 4 /m 4 pl ) would.Therefore, the criterion T < 4 × 10 18 GeV can be seen as a conservative (or sufficient) condition for the highest temperature at which General Relativity together with the Standard Model can be employed as a self-consistent effective theory.
Returning to the phenomenological side, we end by remarking that reaching a high temperature after inflation does not guarantee a substantial thermal gravitational wave background.

Figure 1 :
Figure 1: Diagrams leading to the pair production of gravitons (double wiggly lines) from (a) scalars (dashed lines); (b) fermions (solid lines); (c) gauge fields (wavy lines).We have omitted arrows from scalars and fermions, understanding that both particle and antiparticle states are to be included.
Figure2: Left: the single (h)[13] and double (hh) graviton production rates from a Standard Model

Figure 3 :
Figure 3: The time-integrated contribution of single and double-graviton production to the differential