Freezing-In Gravitational Waves

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


I. INTRODUCTION
The first detection of gravitational waves (GWs) from back hole and neutron star mergers [1,2] opened up a new window to explore our universe.While the GWs that have been detected so far were emitted in the late-time universe, GWs can also be produced in the early universe.These GWs are stochastic in nature and their detection would yield unprecedented information about early universe cosmology as well as high energy particle physics.To give a few examples, GWs in the early universe can be produced from inflation [3][4][5][6], preheating [7,8], inflaton annihilation into gravitons [9][10][11], first-order phase transitions [12,13], cosmic defects such as cosmic strings [14,15], noisy turbulent motion [16][17][18][19][20] and equilibrated gravitons [21,22].For a review on early universe GW sources see Ref. [23].The full GW spectrum for a specific particle physics model that can describe the entire cosmological history was worked out in Ref. [24].
In this paper we consider GWs that were produced from the thermal plasma [25][26][27][28] in the early universe.Every plasma, even in thermal equilibrium, produces GWs due to microscopic particle collisions and macroscopic hydrodynamic fluctuations, cf.Ref. [26].In the former case the GW momenta are on the order of the temperature, k ∼ T and in the latter case they are much smaller, k ≪ T .Here, we focus on GWs produced by microscopic particle collisions, since it enables us to probe elementary particle physics theories at high energies.Furthermore, the GW contribution from microscopic particle collisions to the final spectrum is larger compared to the contribution from hydrodynamic fluctuations [26].
Our main assumption is that after the hot Big Bang a thermal plasma of particles in thermal equilibrium at a maximum temperature T max was present.In addition we assume that at this time no GWs are present.In an expanding universe GWs from the thermal plasma are continuously produced as the temperature decreases.The spectrum of the produced GWs peaks at a frequency on the order of the temperature at the time of production.If the redshift-temperature relation is linear, the GW spectra that are produced at different temperatures add up such that the observed GW spectrum today is enhanced.The spectrum of the produced GWs today peaks in the microwave regime and is hence dubbed the cosmic gravitational microwave background (CGMB).
In principle, the maximum temperature, T max , of the thermal plasma can be as high as the Planck mass m p ≈ 1.2 × 10 19 GeV [29].However in slow roll inflationary cosmology it cannot be much higher than 10 −3 m p .This bound follows by first inferring the energy scale of inflation from the amplitude of scalar perturbations and the tensor-toscalar ratio and then assuming an instantaneous and a maximally efficient reheating [30] to a radiation dominated universe, cf.Ref. [28] and Ref. [31] for a review.Note that T max > 10 −3 m p can be achieved in non-inflationary scenarios.One particular example are bouncing cosmology scenarios which can lead to T max which goes up to the Planck scale: T max < m p , cf.Refs.[32][33][34].The maximum temperature of the thermal plasma is also bounded from below.The most conservative estimates set a lower limit around a few MeV [35][36][37][38][39], shortly before Big Bang nucleosynthesis took place.However, most scenarios require temperatures reaching well above the electroweak scale such that, e.g., sphalerons can be active in leptogenesis scenarios [40].
In previous works the CGMB spectrum has been calculated within the Standard Model (SM) [26,27] and for Beyond Standard Model (BSM) theories [28,[41][42][43].Those works considered only single graviton production processes.In this case the resulting GW energy density per logarithmic momentum interval, Ω gw , is proportional to g 2 × (T max /m p ), where g is the internal coupling in the thermal bath.Here, we extend previous works by also including GW production processes with two gravitons in the final state.These give a contribution to Ω gw that is proportional to (T max /m p ) 3 .
We refer to this GW production channel as graviton pair production.Depending on the values of T max /m p and g, the graviton pair production channel can be the dominating contribution to the CGMB spectrum.In analogy to dark matter production from the thermal plasma we dub the GW production from the thermal plasma GW freeze-in production.
We also identify at which order quantum gravity and back-reaction effects would appear in the CGMB spectrum.Observing these effects in the CGMB spectrum would therefore probe the quantization of gravity and reveal fundamental information about particle physics, if the GW production occurs at high energy scales that cannot be probed with particle colliders on Earth.
Throughout this paper we work with a complex scalar field with quartic coupling λ that is the internal coupling in the thermal bath, i.e., for our model the previously mentioned generic coupling g is the quartic coupling λ.In previous works [26][27][28] such a coupling has not been considered even though the SM has such a coupling in the Higgs sector.That is because Refs.[26,27] worked under the assumption that in the SM the three gauge couplings and the top Yukawa are of order of the square root of the Higgs self-coupling.Note that in BSM theories this is not necessarily the case.This paper is organized as follows: in Sec.II we introduce our model, which is a complex scalar field coupled to gravity.This is then followed by Sec.III, where we introduce the full evolution equations for the two distribution functions f ϕ and f h which describe the scalars and gravitons, respectively.Furthermore, we perturbatively expand the distribution functions around their initial states.This enables us to find a solution of the coupled nonlinear integral-differential equations for the distribution functions.In Sec.IV we calculate the Matrix elements squared for the graviton production processes.We then compute the GW spectrum in Sec.V. Finally, conclusions are given in Sec.VI.Throughout this paper we use natural units with ℏ = c = k B = 1, where k B is the Boltzmann constant.

II. SCALAR MODEL
The action for a complex scalar field on curved space-time is where ∇ µ is a covariant derivative, g µν the metric tensor, g = Det [g µν ] and U is the potential.The flat space-time metric is defined as η µν ≡ diag(−1, 1, 1, 1).Note that for a scalar field the covariant derivative reduces to a partial derivative: ∇ µ ϕ = ∂ µ ϕ.The considered complex scalar is not charged under a local transformation and we consider a quartic potential U = − λ 4 |ϕ| 4 .We assume that the scalar field is massless, which is justified if the considered temperature in the thermal plasma is larger than the mass of the scalar field.
On top of the action in Eq. (1) we need the Einstein-Hilbert action where R is the Ricci scalar and G ≡ 1/m 2 p is the gravitational constant.In the following we expand the metric around flat-space-time: g µν = η µν + h µν , with h µν ≪ 1.A detailed expansion has been worked out before in Ref. [44] and yields Note that the h-fields in Eq. ( 3) have been rescaled with a factor κ ≡ √ 32πG and have now mass dimension one.Furthermore, we have adopted the so-called transverse-traceless (TT) gauge, which includes the De Donder gauge: ∂ α h α µ = 1 2 ∂ µ h together with the requirement that the trace h = h µ µ is zero.In the first and second line of Eq. (3) we wrote down the zeroth and first order terms coming from L EH and L ϕ .In the second line we only write down the second order term coming from L ϕ , since the second order term from L EH will not be needed for our calculations.The lowest order Feynman vertices for our theory are shown in Fig. 1.
FIG. 1: Lowest order vertices that arise from the expansion of the scalar field and Einstein-Hilbert Lagrangian.Scalars are represented by dashed lines and gravitons by double lines.

III. EVOLUTION EQUATIONS FOR THE DISTRIBUTION FUNCTIONS
We describe the thermal plasma of ϕ-particles and the produced gravitons with two distribution functions defined as V is the considered volume, and N k ϕ and N k h are the numbers of ϕ-states and gravitons with momentum k = |k| in the interval d 3 k.Since we shall expand f ϕ around an isotropic equilibrium state, f h is understood to be the polarization-averaged distribution function.Also we do not introduce a distribution function for ϕ † since our model and initial conditions are CP symmetric and therefore it would always be equal to f ϕ .
In the regime where the momentum k is on the order of the hard scale which in equilibrium corresponds to the temperature, i.e. k ∼ T , kinetic theory is expected to be a good approximation for our system.The evolution equations for the ϕ and graviton distribution functions can thus be written in the following Boltzmann-like form where the G and L terms describe the gain and loss terms of particle states.A generic expression for the graviton production term G h is given by all processes r with at least one final state graviton where the index r labels all possible processes.We call the momenta of the incoming ϕ and graviton states respectively.The momenta of the outgoing ϕ's and gravitons are p 1 , • • • , p i and In our notation the m incoming and i outgoing ϕ-states can be ϕ or ϕ † states.In Eq. ( 7) the symmetry factor S r has to be included if two or more indistinguishable particles appear in the initial or final state.We will make the symmetry factor explicit in section Sec.IV where we calculate the graviton rate.The sum in Eq. ( 7) runs over combinations of all processes with at least one graviton with momentum k in the final state.The pre-factor 1/(4k) is a combination of 1/(2k) from the phase space measure and 1/2 from the graviton polarization degeneracy.We need the factor 2 from the polarization degeneracy since the matrix element squared is summed over polarizations and the distribution function is defined to be averaged over both polarizations.The loss term is analogous to the gain term in Eq. (7) with the difference that one sums over all processes with at least one graviton in the initial state.The Boltzmann-like Eqs.
(5), ( 6) and (7) come with important caveats on their validity beyond leading order, which we shall discuss later.The integral that appears in G h is the phase space integral that has to be performed over all momenta, except k: where we use capital letters to denote four-vectors.For further use we introduce the shorthand notation: We have written Eqs.( 5), ( 6), ( 7) and ( 8) in a rather generic form which includes all possible processes.In our specific model of a massless complex scalar field, 1 ↔ 2 processes are only allowed in the collinear limit, i.e., when the three-momenta of all three particles are exactly parallel.However, in this case the thermal and vacuum masses of the scalars have to be taken into account, which leads to the fact that 1 ↔ 2 graviton production processes are not even allowed in the collinear limit.The first kinematically allowed processes for graviton production are 2 → 2 and 2 ↔ 3 processes, which have a lowest order matrix element squared of O κ 4 and O λ 2 κ 2 , respectively.We are interested in tracking the evolution of f h under the assumption that it starts from an initially vanishing value in a bath of equilibrated scalars. 1 In our stated freeze-in scenario, we further assume that throughout the entire evolution f h ≪ 1 and |n B − f ϕ | ≪ 1, where n B is the Bose-Einstein distribution n B (k) ≡ 1/(e k/T − 1).We can thus expand the distribution functions up to fourth order in λ and κ: where the superscript stands for the order of λ and κ that is considered, i.e.O(λ i κ j ) = (i, j).Note that in our expansion of the distribution functions we have also implicitly expanded the matrix element squared.As κ is dimensionful, the expansion in κ ∼ 1/m p has to be understood as an expansion in T κ, the corresponding dimensionless quantity.The zeroth order term of f ϕ is set to be n B and the zeroth order term of the f h distribution function is zero.We have suppressed the time arguments in Eqs. ( 10) and (11).Note that in the expansion we treat λ and T κ on equal footing and for the distribution functions we have only written out the non-vanishing terms up to fourth order.Terms of O λ n κ 2 with n > 2 are also expected, which -depending on the temperature and the value for λ -can be larger than O κ 4 .However, as we discuss later in this section, these terms cannot in general be included in a straightforward manner into the Boltzmann-like ansatz, cf.Eqs. ( 5), ( 6) and (7).
In the following, we discuss based on three examples why the terms shown in Eqs.(10) and ( 11) are the only nonzero contributions.The evolution Eqs. ( 5) and ( 6) yield ḟ (2,0) h = 0 since, in order to produce or annihilate a graviton, one has to go at least to second order in κ.Furthermore, ḟ (0,2) (2,0) . After a redefinition of variables in the phase space integral one can show that both terms cancel each other.Similar arguments hold for the other terms in ḟ (2,0) ϕ such that overall ḟ (2,0) ϕ = 0. Note that if ḟ (i,j) = 0 then f (i,j) = 0 for all times for (i, j) ̸ = (0, 0) which follows from the initial conditions: f h (t = 0) = 0 and f ϕ (t = 0) = n B .In complete analogy to the discussed examples one can show with kinematic and detailed balance arguments that the other terms that are not shown in Eqs. ( 10) and ( 11) also vanish.
Next we discuss the non-zero fourth order terms for the graviton distribution function.The ḟ (2,2) h rate is non-zero and it is given by where the dots stand for other processes of the same order, e.g., ϕϕ † ϕ → ϕh.All possible processes are written down in Sec.IV.In Fig. 2, we show the Feynman diagrams for the single graviton production process ϕϕ → ϕϕh at order FIG. 2: Lowest order Feynman diagrams for the ϕϕ → ϕϕh process.The matrix element squared is of the order O λ 2 κ 2 .In addition to the diagrams shown, there are two more diagrams that can be obtained by crossing the final state graviton.
FIG. 3: Lowest order Feynman diagrams for the process ϕ † ϕ → hh.The matrix element squared is of the order O κ 4 .There is one additional diagram that is not shown since it can be obtained by crossing the final state gravitons in the second diagram.
The contributions to ḟh at order O κ 4 are sourced by processes that have two gravitons in the final state.These processes can be relevant in high-temperature early universe scenarios since the dimensionless expansion parameter is T /m p which can be relatively large if the temperature is close to the Planck mass.The explicit form of ḟ (0,4) Processes with a graviton in the initial state do not contribute, since the initial state graviton always comes with a factor f h which makes the whole term of higher order.The same holds for final state 1 + f h amplification factors.
The corresponding Feynman diagrams that contribute to two graviton production are shown in Fig. 3.In Sec.IV we calculate these diagrams.
In the following we discuss how one could extend our calculation of f h to higher orders.In particular we point out the limitations and challenges that one would face.The evolution Eq. ( 5) for the graviton distribution function can definitively be used without problems at lowest order in perturbation theory, i.e., in our case these are the 2 → 2 and 2 ↔ 3 processes that are of the order O κ 4 and O λ 2 κ 2 respectively.These lowest order terms have real corrections, and virtual corrections.Virtual corrections are loop correction, while real corrections come from tree-level processes with extra initial-or final-state particles.If they are finite, the latter are easily incorporated into the Boltzmann equation formalism, i.e., Eq. (7).Incorporating the former, on the other hand, is not straightforward as the matrix elements squared contain not only standard vacuum fluctuations but also statistical fluctuations which, in turn, depend on the distribution functions themselves [45].Furthermore, while real and renormalized virtual corrections might separately be finite in a standalone scalar theory, in more complex systems such as the scalar theory coupled to gravity or gauge theories they are in general not finite, with infrared (IR) divergences canceling between the two, as in the case of the Kinoshita-Lee-Nauenberg theorem [46,47].In conclusion it is a challenging task to incorporate higher order effects with the Boltzmann-like approach since it is only possible to incorporate the finite higher order effects.
FIG. 4: The interference of the left and center diagrams yields a term O λ 2 κ 4 in the matrix element squared.The right diagram denotes a real correction that is on the same order when it's squared.
Virtual gravitons arise already at order O κ 4 in the graviton distribution function, cf.Fig. 3.Quantum gravity effects start to play a role at order O λ 2 κ 4 , since at this order diagrams with graviton loops exist.The three diagrams that we show in Fig. 4 are of order O λ 1 κ 1 , O λ 1 κ 3 and O λ 1 κ 2 respectively.The interference term of the first two diagrams is of order O λ 2 κ 4 and is the first virtual correction involving loops of gravitons.Conversely, the square of the third diagram is O λ 2 κ 4 and is part of the real corrections at that order.This further exemplifies the challenge in going beyond leading order: the real corrections can be dealt with in a Boltzmann form in a rather straightforward way, while the virtual corrections cannot, as their matrix element squared will depend in non-trivial ways on the statistical factors -see Ref. [48] for a recent work on this problem in a non-gravitational setting.
In an alternative way, one could systematically study quantum effects in kinetic theory from first principles using the Wigner-function formalism.By performing an expansion in ℏ of the Wigner function, one can in principle derive quantum corrections to the classical Boltzmann equation, see, e.g., Refs [49][50][51][52][53].The development of such quantum kinetic theory is left for future studies.
Along the lines that we discussed before, we want to mention that it is possible to calculate ḟ (all,2) h in a full quantum picture.As shown in Ref. [27] for the specific case of GWs and more generally in Ref. [54] for any state that is feebly coupled to a thermal bath, the Boltzmann-equation based approach that we use here agrees, at O λ 2 κ 2 , with the thermal-field-theoretical approach of production and equilibration rates.Namely, for GWs one has ḟ (all,2) [26] where the production/equilibration rate Γ(k) is proportional to κ 2 times the imaginary part of the retarded two-point function of the T 12 component of the energy-momentum tensor of the equilibrium particles, i.e., in our case, the scalars.This formalism defines the single-graviton production rate to all orders in λ.Within this formalism, higher orders in λ naturally incorporate both real and virtual corrections, without the issues that would plague direct attempts in the Boltzmann-like approach.However note that, in principle, we do not want to go to higher orders in λ but to higher orders in κ to identify quantum gravity effects.The discussed thermal-field-theoretical approach is not suited for this and new strategies have to be developed for a full quantum treatment of the graviton production rate.
Finally let us discuss back-reaction effects, which can be incorporated into the Boltzmann-like formalism.If we stay at order O κ 4 and go to non-zero order in λ we can identify back-reaction effects.These appear at lowest order at O(κ 4 λ 4 ).The graviton production rate at this order contains the following back-reaction terms: where the dots above stand for other terms that we have omitted here.We call the terms in Eq. ( 14) back-reaction terms since the small corrections on top of the Bose Einstein distribution, f (2,2) ϕ and f (0,4) ϕ , appear in the phase space integral.Back-reaction effects also appear in the ḟϕ rate.
We further note that the RHS in Eqs. ( 12) and ( 13

IV. MATRIX ELEMENTS AND PHASE SPACE INTEGRALS
In this section we calculate the matrix elements squared for graviton production at order O λ 2 κ 2 and O κ 4 .Let us start with the O λ 2 κ 2 component.As argued previously, it arises from 2 → 3 and 3 → 2 processes.The corresponding expressions for the distribution functions are: where at order O λ 2 κ 2 there is no f h on the right-hand side.The sums run over all abcd scalar and antiscalar degrees of freedom and thus over all ab → cdh and abc → dh processes, with h denoting the graviton.The quantities are the corresponding matrix elements squared summed over the graviton polarizations.For k ∼ T the contribution of the thermal mass m ϕ T = λ/12T is suppressed, so the external states can be considered massless.The prefactor 1/(16k) is a combination of 1/(2k) from the phase space measure, 1/2 for the graviton polarization degeneracy, and 1/(2!) 2 for the symmetry factors for identical initial and final state particles.In the cases where a ̸ = b or c ̸ = d the sum over abcd counts the process two times and compensates for this factor.Similarly, 1/(24k) is a combination of 1/(2k) from the phase space measure, 1/2 for the graviton polarization degeneracy, and 1/(3!) for the symmetry factors for identical initial state particles.
The phase spaces can be read off from Eq. ( 8).For the matrix element squared we used the automated pipeline introduced in Ref. [27].We first used FeynRules [55] to derive Feynman rules for the Lagrangian in Eq. ( 3).Using the appropriate interface [56], FeynRules generates a model file for FeynArts [57].This package and its companion FormCalc [58] were then used to generate, evaluate and square all amplitudes.Tensor boson polarization sums had to be implemented following the method discussed in Ref. [27].For ϕϕ † → ϕϕ † h, the results is Equation ( 18) arises from diagrams such as the ones in Fig. 2. The four structures in the denominator, e.g.P ′ 1 • K, correspond to the propagator of the virtual, intermediate scalar connecting the ϕϕh vertex with the ϕ 4 one.It is reassuring to see that this matrix element squared is finite even when one of the scalar products in the denominators vanishes.For example, the term P ′ 1 • K can vanish either for p ′ 1 → 0 or when p ′ 1 is parallel to k.In the former case the powers of p ′ 1 at the numerator immediately remove the divergence, and similarly the phase space is free of endpoint divergences for p ′ 1 → 0, even in the presence of Bose enhancement (n B (p ′ 1 → 0) ≈ T /p ′ 1 ).In the collinear case, p ′ 1 ∥ k, one can show that the sum of the divergent terms is finite.
Let us now discuss the terms that are obtained by crossing.By crossing an initial state ϕ (ϕ † ) to the final state one obtains a ϕ † (ϕ).Hence one finds Furthermore, M ϕϕ † ϕϕ † h 2 is symmetric under permutations within the initial and final states of the ϕ and ϕ † .Hence, in the sum over abcd of Eq. ( 15), it is counted four times, whereas the ϕ-only or ϕ † -only processes are counted once.
The 3 → 2 matrix elements squared can be obtained by crossing Eq. ( 18), too.For instance, if we cross the final-state ϕ † with momentum P 2 into an initial-state ϕ with momentum P ′ 3 we have By further crossing one finds The processes with two ϕ and with two ϕ † in the initial state are counted three times each in the sum over abcd.
Putting everything together we then find Details on the seven-dimensional numerical integration of the phase space can be found in Appendix A.
The O κ 4 component is sourced by 2 → 2 processes: where 1/(8k) is the product of the 1/(2k) from the Lorentz phase space measure, a factor of 1/2 for the two polarizations of the graviton and a factor of 1/(2!) for possible identical initial state particles.The matrix element squared arises from four diagrams.Three of them are shown in Fig. 3, and the fourth comes from the u-channel analog of the second diagram.We have computed the matrix element squared by using the previously described FeynRules, FeynArts and FormCalc machinery: where s, t and u are the standard Mandelstam invariants.Note that the matrix element squared can also be extracted from Refs.[59, 60]2 and we have checked that our results agree.Accounting for a factor of 2 from the sum in Eq. ( 24) yields: We refer again to Appendix A for details on the reduction of the phase space integration to a two-dimensional integral that we evaluate numerically.The production rates for (i, j) = (2, 2) and (i, j) = (0, 4) can be written in a compact form as: where we have defined dimensionless ψ-functions and used the convention that n B (x) = 1/(e x − 1) if the argument of n B is dimensionless.The ψ-functions are shown in Fig. 5 for the 2 → 3 (dotted black) and 3 → 2 (dashed black) processes.We also show the sum of the 2 → 3 and 3 → 2 processes which is labeled as 2 ↔ 3 and shown as a FIG.5: We plot the ψ-functions that have been obtained by numerically integrating the phase space integrals.The results are shown for the 2 → 3, 3 → 2 and 2 → 2 processes respectively.The sum of the 2 → 3 and 3 → 2 contributions is denoted as 2 ↔ 3. 2 ↔ 3 processes come from single graviton production processes that have a lowest order matrix element squared of the order O λ 2 κ 2 = (2, 2). 2 → 2 processes arise in graviton pair production processes with a lowest order matrix element squared of the order O κ 4 .red dot-dashed line.The 2 → 2 processes are shown as a solid blue line.Note that while the ψ-function for the 2 ↔ 3 processes has only a relatively mild k/T dependence around k/T ≃ 1 this is not the case for the 2 → 2 processes.In Appendix A we have derived an asymptotic form for ḟ (0,4) h , i.e., ψ (0,4) in the limit k > T .We find ψ (0,4) = 32/(15π)(k/T ) 2 for k > T what confirms the results in Fig. 5. ψ (2,2) and ψ (0,4) differ in their functional form due to the fact that the underlying phase space structure and matrix elements are different for single and graviton pair production.The single gravitons are produced in 2 ↔ 3 processes, while graviton pair production is a 2 → 2 process.In addition, the non-renormalizable nature of gravity leads to a matrix element squared for single graviton production that is proportional to κ 2 ∼ 1 m 2 p times an expression with dimensions (energy) 2 .Similarly the matrix element squared for graviton pair production is proportional to κ 4 ∼ 1 m 4 p times an expression with dimensions (energy) 4 .In the case of graviton pair production the 2 → 2 phase space integration for large k translates this into a quadratic dependence on momentum, which explains the form of ψ (0,4) for large k.
As Fig. 5 shows, the 2 → 3 contribution to ḟ (2,2) which implies a naively IRdivergent contribution to the number density of gravitons (n gw ∝ dkk 2 f h ) and a finite but enhanced contribution to the energy density ρ gw ∝ dkk 3 f h from the IR domain k ≪ T .This IR contribution is an artifact of treating the external and intermediate scalar states as massless.If we would include their thermal mass m ϕ T , then such a behavior would go away, as the matrix element would no longer behave like 1/k 2 at small k.Scalars have the nice property that their thermal mass behaves like an ordinary local mass term in the Lagrangian, unlike gauge fields and fermions.This in turn would make thermal mass resummation relatively straightforward.For this paper we limit ourselves to unresummed (massless) results, and consider the result for ḟh to be a proper leading-order determination in the regime k ≳ m ϕ T = λ/12 T .Furthermore, for smaller k, k ≪ λ 2 T , the quasi-particle description breaks down completely and gravitational waves are sourced from hydrodynamic fluctuations [26].

V. GRAVITATIONAL WAVE SPECTRUM
In this section we embed the graviton production rate into cosmological evolution.Our main assumption is that, after the hot Big Bang, a thermal plasma of ϕ-particles with temperature T max is present.Throughout the expansion of the universe the thermal plasma produces GWs.From the definition of Eq. ( 4) it follows that the GW differential energy density is , where a flat space-time metric has been assumed, and the factor of 2 takes the two polarization states into account which contribute to the energy density.We can rewrite the equation for the energy density as dρgw dt dlnk = k 4 π 2 ḟh .Generalizing to an expanding universe, the GW energy density evolves as [26] ( where H is the Hubble parameter and we have defined R(t, k) ≡ 2 k ḟh (t, k).Note that for the O λ 2 κ 2 and O κ 4 contributions, which we shall discuss here, ḟh has no explicit time dependence.We will therefore treat ḟh without explicit time dependence in the following derivation.Now that we consider a thermal plasma in an expanding universe, the temperature decreases over time.Therefore we have an implicit time/temperature dependence.We further note that, in a radiation-dominated universe, with ρ(T ) = g * ρ (T )π 2 T 4 /30 and g * ρ (T ) = 2 is the effective number of energy density degrees of freedom.The scalars remain in thermal equilibrium 3 as long as their interaction rate, which is on the order λ 2 T , is at least as fast as the Hubble rate H ∼ T 2 /m p .Therefore we obtain the equilibrium condition, λ 2 ≳ T /m p , which has to be understood as an order of magnitude estimate.We can integrate Eq. ( 28): where we have used that the entropy density s fulfills ṡ+3Hs = 0 and we have made the temperature/time dependence explicit.We assume that at the beginning of the GW production no GWs are present, i.e., ρ gw (t 0 ) = 0.At t 0 the thermal plasma was first in thermal equilibrium and had its maximum temperature T max .The time t 1 is the time when the mass of the scalar field cannot be neglected anymore.The temperature corresponding to the time t 1 is referred to as T ϕ in the following.In Eq. ( 30) we can only integrate to t 1 since the production rates that we have calculated are only valid for temperatures above T ϕ .The time integral in Eq. ( 30) can be transformed into an integral over the temperature by using the relation [61]: where g * s and g * c are the effective degrees of freedom for the entropy density and heat capacity, which are defined as s(T ) ≡ g * s (T ) 2π 2 45 T 3 and c(T ) ≡ g * c (T ) 2π 2 15 T 3 .In the following we use the assumption of isotropy under which we can simplify the d 3 k integral: T, k).From Eq. ( 30) we can then read off the GW energy density per logarithmic momentum interval at T ϕ , normalized to the total energy density: Ω gw = 1 ρ dρgw dlnk .Redshifting all corresponding quantities to today [28] yields: where f g is the current day GW frequency, h 2 0 Ω γ = 2.473 × 10 −5 is the present fractional photon energy density, h 2 0 a factor that eliminates the experimental uncertainty that is coming from measurements of the Hubble constant, 3 By this we mean that the zeroth order term of f ϕ is a massless Bose-Einstein distribution with the current temperature.
T today = 2.7254 K the current day temperature [62] and g * s (T today ) = 3.931 are the effective entropy degrees of freedom today [63].
With the parametric form of ḟh from Eq. ( 27) we can write R as: where the dots denote higher order terms.We plug Eq. ( 33) in Eq. ( 32) and, in order to get an analytical result, we approximate g * s (T ) = g * ρ (T ) = g * c (T ) = g * (T max ) in the region of temperatures above T ϕ .We thus obtain where we have assumed T max ≫ T ϕ and we have defined y max ≡ .
Note that models which can describe the entire thermal history of the early universe have g * (T max ) > g * s (T today ).We work with a model that includes only one complex scalar field and therefore, g * (T max ) < g * s (T today ).Nonetheless the features in the GW spectrum that we work out here will hold in general even with a more realistic model that can describe the thermal history of the universe consistently.We comment on this aspect further below.The terms that are represented by the dots in Eq. ( 34) include processes which encode quantum gravity effects, cf.Fig. 4.These effects arise at the order O λ 2 κ4 and would appear as a term ψ (2,4) λ 2 (T max /m p ) 2 in the parantheses in the second line in Eq. ( 34).The single graviton production processes have a matrix element squared of the order O λ 2 κ 2 , i.e., it is proportional to 1/m 2 p .Since the GW production happens on a time scale that is comparable with m p , the single graviton production processes are proportional to T max /m p in the GW spectrum.A similar argument applies to the contribution from the graviton pair production.The matrix element squared is of the order O κ 4 , i.e., proportional to 1/m 4 p , and hence the contribution to the GW spectrum is of the order (T max /m p ) 3 .We refer to this production mechanism as GW freeze-in, since the GWs are produced from the thermal plasma throughout the expansion of the universe.Note that while the GW production mechanism is conceptually very similar to the production of dark matter from the thermal plasma, the final GW spectrum is very ultraviolet sensitive in the sense that it depends on the maximum temperature.
In Fig. 6 we plot the GW spectrum coming from single graviton production processes (red dot-dashed lines) and from graviton pair production processes (blue dashed lines).The sum of both contributions, i.e., the total GW spectrum, is shown as a black solid line.The quartic coupling is always set to λ = 10 −1 . 4When we evaluate the GW spectrum we also have to evaluate the ψ-functions.Since these functions have only been calculated reliably for arguments that are larger than √ λ = 0.31, cf.Sec.IV, we only show the GW spectrum in the corresponding frequency regime.Note that the GW spectrum from the single graviton production mimics to a good approximation a black body spectrum, since the function ψ (2,2) is very flat in the regime T /m p > √ λ.The graviton pair production contribution has a significantly different shape since the function ψ (0,4) is not constant in the frequency interval of interest.
Fig. 6 (left) shows a scenario where the maximum temperature is set to T max /m p = 10 −3 .In this case the single graviton production contribution dominates over the graviton pair production contribution.The total spectrum has therefore mostly the form of the single graviton production spectrum.
Fig. 6 (right) shows the GW spectrum for a slightly larger maximum temperature, T max /m p = 10 −2 , which requires a non inflationary scenario.Note that the chosen maximum temperature is consistent with the equilibrium condition for the scalar fields, i.e. λ 2 ≳ T /m p .In the case at hand both contributions are parametrically equally important and around the peak frequency the graviton pair production contribution is even substantially larger than the single graviton production contribution.This can be seen explicitly from Eq. ( 34) by comparing the two terms in the second line.The first term corresponds to single graviton production processes, while the second one describes the contribution due to graviton pair production.Comparing both terms we find as an order of magnitude estimate, that the contribution from graviton pair production processes are equally important or even larger than the contribution from the single graviton production processes if 10 T max /m p ≳ λ.Therefore, the relevance of the graviton pair production processes depends crucially on the size of the coupling and the maximum temperature.The GW spectra associated with single graviton and graviton pair production processes peak at slightly different frequencies and have a distinct functional form.As a result, the total spectrum takes a very characteristic form that is substantially different from the single graviton production spectrum, i.e., an approximate black body spectrum.
In the following we derive analytic expressions for the peak frequencies of the single graviton and graviton pair production GW spectra.For the single graviton contribution we find from Eq. ( 34) that the GW spectrum peaks at: 2 × 10 11 Hz (2/g * s (T max )) 1/3 , where we have assumed that ψ (2,2) = const.which is a good approximation in the frequency interval of interest.The peak frequency of the graviton pair production curve lies at a slightly higher frequency: 3.5 × 10 11 Hz (2/g * s (T max )) 1/3 , where we have used the asymptotic form ψ (0,4) = 32/(15π)(k/T ) 2 that was derived in Appendix A. The higher order terms which are depicted by the dots in Eq. ( 34) are suppressed further by powers of T max /m p and λ.For values of T max /m p and λ that are relatively close to unity one has to check in detail that the higher order ψ functions are not larger than the leading order ψ functions such that the suppression from the additional powers of T max /m p and λ is not spoiled.The values of T max and λ, that we consider in Fig. 6 are much smaller than unity and therefore we do not expect such a scenario to happen.For example the ψ (0,6) function would have to be four orders of magnitude larger than the ψ (0,4) function for T max /m p = 10 −2 .The higher order ψ functions will be either phase space suppressed or have the same phase space as the leading order processes.In both cases we do not expect an enhancement of the higher order ψ functions by orders of magnitude.
The maximum of the CGMB spectrum is bounded from above, h 2 0 Ω gw ≲ 10 −6 [28], due to constraints on the additionally allowed amount of dark radiation.Both scenarios that are shown in Fig. 6 do not saturate this bound and are therefore not excluded.In the complex scalar model the dark radiation bound is saturated for T max /m p ≃ 0.5.Note that in this case the main contribution to the GW spectrum is coming from graviton pair production processes, which illustrates the importance to include these processes at high temperatures.The SM predictions for the GW production from the thermal plasma includes currently only single graviton production processes, cf.Refs.[26][27][28], and yields h 2 0 Ω gw = 5 × 10 −7 for T max = m p .Therefore, adding graviton pair production processes to the SM calculation can already lead to a violation of the dark radiation bound for T max < m p , because the current prediction which includes only single graviton production processes, is already very close to the dark radiation bound.In conclusion, this might then be used to constrain the maximum temperature of the universe.We plan to work out the details in a follow-up study where we will also address the aspect of thermal equilibrium at ultra high temperatures in the SM and beyond the SM theories.

VI. CONCLUSIONS AND OUTLOOK
The thermal plasma in the early universe produced a guaranteed stochastic GW background through thermal fluctuations.At each time the emitted GW spectrum peaks at the respective temperature.Due to the temperatureredshift relation, the peak frequencies of the GW spectra are all redshifted to the same frequency today and therefore add up.Conceptually, the GW production from the thermal plasma has many similarities with the so-called dark matter freeze-in production from the thermal plasma.The GWs are produced out of equilibrium and their distribution function is small at all times, f h ≪ 1.Furthermore, the f h distribution function evolves much slower than the Hubble rate.We therefore dubbed the GW production from the thermal plasma GW freeze-in production.The GW freezein scenario is ultraviolet dominated in the sense that it depends on the maximum temperature of the universe, as expected from a non-renormalizable coupling.

total
(2,2) single graviton prod.(0,4) graviton pair prod.FIG.6: Gravitational wave spectrum with respect to the present day GW frequency.The single graviton production processes which are of order O λ 2 κ 2 = (2, 2) are shown as a red dot-dashed line.The graviton pair production processes are shown as a dashed blue line and are of order O κ 4 = (0, 4).The total GW spectrum is shown as a solid black line.On the left we show a scenario where the maximum temperature is limited to 10 −3 m p .In the right figure we set the maximum temperature to T max = 10 −2 m p .In this case the graviton pair production processes yield an even larger contribution compared to the single graviton production processes.The GW spectra are calculated for a complex scalar model with g * (T max ) = 2.
In this paper we use a Boltzmann-like formalism to study the microscopic particle collision processes that contribute to the CGMB spectrum.We have done all calculations in a model with a complex scalar field and quartic selfinteraction.Our basic assumption is that after the hot Big Bang a plasma of scalars with temperature T max is present and this plasma produced the CGMB spectrum.First, we considered the contribution of single graviton production processes to the CGMB spectrum.In a scalar theory with quartic interaction, single graviton production processes are 2 ↔ 3 processes, which have not been calculated before.Our calculation is motivated by the fact that a quartic coupling exists in the Higgs sector of the SM and in many BSM theories.The second class of processes that we investigate are graviton pair production processes.These are 2 → 2 processes and have not been considered before in the context of GWs from the thermal plasma.We show that their contribution to the CGMB spectrum can be larger than the contribution from the single graviton production processes.As an order of magnitude estimate, graviton pair production processes dominate the GW spectrum if 10 T max /m p ≳ λ.Note however that the maximum temperature is also bounded from above by an equilibrium requirement for the scalar particles: λ 2 ≳ T max /m p , which has to be seen as a parametric estimate.Therefore, the degree to which graviton pair production processes contribute significantly to the CGMB spectrum depends on the values of the coupling coefficient and the maximum temperature.As an example we show the two different scenarios in Fig. 6.On the left single graviton production processes dominate (T max /m p = 10 −3 and λ = 10 −1 ).When increasing T max by one order of magnitude (Fig. 6, right) graviton pair production processes yield a significant contribution to the total GW spectrum.
The single graviton and graviton pair production processes are the lowest order contributions which can easily be incorporated into our Boltzmann-like formalism.We have also discussed the first steps and problems that would arise if one would add real and virtual quantum gravity corrections to the presented results.While finite real corrections can be incorporated in our formalism, virtual corrections depend on the distribution functions themselves and this complicates their inclusion into the Boltzmann-like approach that we use here.A possible future direction is to derive a quantum Boltzmann equation from the Wigner function by performing a systematic ℏ expansion.This would allow one to explicitly identify the quantum corrections.
The results that we have worked out for a scalar model are qualitatively also valid for more general theories.In this case the coupling coefficient λ would have to be replaced with the heat bath couplings in the more general theory, which we generically refer to as g.Then the contribution from graviton pair production processes to the GW spectrum dominates over the single graviton contribution if X × T max /m p ≳ g, where X is a model dependent constant.A value X < 10 in the SM would indicate that graviton pair production dominates at higher temperatures compared to our scalar model.In a follow-up study we will answer this point with a full SM and BSM calculation.Confronting the graviton pair production calculation in the SM and BSM theories with existing dark radiation constrains can therefore already lead to constraints on T max .Constraints on T max can be used to test different models of our universe, cf.Ref. [28].For example non-standard inflationary cosmological models would be required if a T max ≳ 10 −3 m p would be inferred from the GW spectrum.Furthermore, the CGMB can be used to constrain non-standard cosmological histories, cf.Ref. [64].The authors of Ref. [64] have only considered single graviton production processes.It would be interesting to add graviton pair production processes to their calculation since it could lead to stronger constraints on non-standard cosmological histories.
A detection of the CGMB with Earth-based detectors will be challenging, cf.Ref. [65], however there exist detector proposals with sensitivities comparable to the dark radiation bound, cf.Ref. [66].Future experimental work will have to show if the proposed detectors can be realized and if their foreseen sensitivity can even be improved.Our results motivate further work on high-frequency GW detection since a detection of the CGMB in the future would pave the way to probe our understanding of particle physics and cosmology at ultra high energies.
since massless 1 ↔ 2 processes are kinematically forbidden.Finally, ḟ (2,0) ϕ vanishes because of detailed balance arguments.As an example consider the two terms