{'fraction_non_alphanumeric': 0.037439525441859335, 'fraction_numerical': 0.015972576763182846, 'mean_word_length': 4.095255790907965, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 maximum temperature of the universe in units of Planck mass 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 also show that quantum gravity effects arising in single graviton production are smaller than the leading order result by a factor of the square of the ratio between the maximum temperature and the Planck mass. 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 which have been detected so far were emitted in the late-time universe, GWs can also be produced in the early universe. These GWs are expected to be stochastic in nature, i.e., they have no fixed phase relation and come from all directions. A 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], first-order phase transitions [9,10], cosmic defects such as cosmic strings [11,12], noisy turbulent motion [13][14][15][16][17] and equilibrated gravitons [18,19]. For a review on early universe GW sources see Ref. [20]. The full GW spectrum for a specific particle physics model which can describe the entire cosmological history was worked out in Ref. [21].
In this paper we consider GWs that were produced from the thermal plasma [22][23][24][25] in the early universe. Every plasma, even in thermal equilibrium, produces GWs due to microscopic particle collisions -with GW momenta on the order of the temperature, k ∼ T -and macroscopic hydrodynamic fluctuations in the plasma, for k T [23]. 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 [23].
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 [26], however in inflationary scenarios it cannot be much higher than 10 −3 m p [25]. The maximum temperature of the thermal plasma is also bounded from below. The most conservative estimates set a lower limit around a few MeV [27][28][29][30][31], shortly before Big Bang nucleosynthesis took place. However, most scenarios require temperatures reaching well above the electroweak scale such that, e.g., sphalerons may be active in leptogenesis scenarios [32].
In previous works the CGMB spectrum was calculated by taking into account only single graviton production processes. The calculation has been done for the Standard Model (SM) [23,24] and for Beyond Standard Model (BSM) theories [25,[33][34][35]. The resulting GW energy density per logarithmic momentum interval Ω gw is in this case 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 which 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 for. 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 high energy 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 λ which is the internal coupling in the thermal bath. Note that in previous works [23][24][25] such a coupling has not been considered even though the SM has such a coupling in the Higgs sector. That is because Refs. [23,24] 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 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 simplifies a solution of the otherwise 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. [36] 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.

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 forṁ 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 p 1 , · · · , p m and k 1 , · · · , k n , respectively. The momenta of the outgoing φ's and gravitons are p 1 , · · · , p i and k 1 = k, k 2 , · · · , k j . 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 equations (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 non-zero contributions. The evolution equations (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,ḟ processes are kinematically forbidden. Finally,ḟ (2,0) φ vanishes because of detailed balance arguments. As an example consider the two terms 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ḟ φ such that overallḟ φ = 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 . With kinematic and detailed balance arguments one can prove that the other terms that are not shown in Eqs. (10) and (11) must also vanish.
Next we discuss the non-zero fourth order terms for the graviton distribution function. Theḟ 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 O λ 2 κ 2 . The contributions toḟ h at order O κ 4 are sourced by processes which 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ḟ 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.
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 which can be obtained by crossing the final state graviton.
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 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 which 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 treelevel 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 [37]. 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 [38,39]. 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.
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. [40] for a recent work on this problem in a non-gravitational setting.
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.
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 [41][42][43][44][45]. 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. [24] for the specific case of GWs and more generally in Ref. [46] 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ḟ [23] 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)  φ , 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 |M ab cdh (p 1 , p 2 ; p 1 , p 2 , k)| 2 and |M abc dh (p 1 , p 2 , p 1 ; p 1 , k)| 2 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. [24]. We first used FeynRules [47] to derive Feynman rules for the Lagrangian in Eq. (3). Using the appropriate interface [48], FeynRules generates a model file for FeynArts [49]. This package and its companion FormCalc [50] were then used to generate, evaluate and square all amplitudes. Tensor boson polarization sums had to be implemented following the method discussed in Ref. [24]. 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 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 M φφ † φ φh (p 1 , p 2 , p 3 ; p 1 , k) 2 = κ 2 λ 2 2 (P 1 · P 2 ) 2 P 1 · K P 2 · K + (P 1 · P 3 ) 2 P 1 · K P 3 · K + (P 2 · P 3 ) 2 P 2 · K P 3 · K − (P 1 · P 1 ) 2 P 1 · K P 1 · K − (P 1 · P 2 ) 2 P 1 · K P 2 · K − (P 1 · P 3 ) 2 P 1 · K P 3 · K − 2 .
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 finḋ 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, the fourth comes from the u-channel analogue 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. [51,52] 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 is 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 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. As Fig. 5 shows, the 2 → 3 contribution toḟ which implies a naively IR-divergent 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 were to include their thermal mass m φ T , then such a 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 which 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 .
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 [23].

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 maximum 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 dρ gw (t, k) = 2kf h (t, k) d 3 k (2π) 3 , where a flat space-time metric has been assumed, and the factor of 2 takes into account the two polarization states that 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 [23] 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 , cf. Sec. III, is at least as fast as the Hubble rate H ∼ T 2 /m p . Therefore we obtain the equilibrium condition λ 2 T /m p , where this condition is 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, at t 0 , 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. We refer to the corresponding temperature as T φ . In Eq. (30) we can only integrate to t 1 , i.e. T φ , 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 [53]: 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: d 3 k R(T, k) = 4π dk k 2 R(T, k) = 4π dln(k)k 3 R(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 [25] 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, T today = 2.7254 K the current day temperature [54] and g * s (T today ) = 3.931 are the effective entropy degrees of freedom today [55].
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. The single graviton production processes have a matrix element squared of the order O λ 2 κ 2 , i.e. it is 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 . 4 When 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) changes rapidly. 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 O λ 2 κ 2 dominates over the graviton pair production O κ 4 contribution. The total spectrum has therefore mostly the form of the single graviton production spectrum, with the largest deviation from it at k/T > 1.
In Fig. 6 (right) we show a scenario with T max /m p = 10 −2 . Note that the chosen maximum temperature is still 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 O κ 4 contribution is even substantially larger than the O λ 2 κ 2 contribution. This can be seen explicitly from Eq. (34) by comparing the O λ 2 κ 2 and O κ 4 contributions in the parentheses. As an order of magnitude estimate, we find 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 contributions from the O λ 2 κ 2 and O κ 4 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. Note that a maximum temperature larger than T max /m p = 10 −3 would require a non-inflationary scenario. The maximum of the CGMB spectrum has an upper bound h 2 0 Ω gw 10 −6 [25], which is coming from dark radiation constraints. Both scenarios that are shown in Fig. 6 do not saturate this bound and are therefore not excluded.
Quantum gravity effects arise at the order O λ 2 κ 4 and would be -even in the scenario where the temperature is as high as T max /m p = 10 −2 -small corrections to the total GW spectrum.

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 temperature- 4 In a more realistic model that can describe the entire history of our universe one would have to use renormalization group equations to run the parameters up to very high energy scales.

total
(2,2) single graviton prod. (0,4) graviton pair prod.  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 /m p = 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.
redshift relation, the GW spectra, which are emitted at different times, 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 freeze-in scenario is ultraviolet dominated in the sense that it depends on the maximum temperature of the universe, as expected from a non-renormalizable coupling.
In this paper, using a Boltzmann-like formalism, we have studied the microscopic particle collision processes that contribute to the CGMB spectrum. We have done all calculations in a model with a complex scalar field with a quartic self-interaction. 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 the CGMB. 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 crucially on the values of the coupling coefficient and the maximum temperature.
In Fig. 6 we considered two different scenarios. Fig. 6 (left) shows a scenario where single graviton production processes dominate, i.e., we set T max /m p = 10 −3 and λ = 10 −1 . In this case the contribution from graviton pair production is a small correction to the single graviton production. However increasing T max /m p by one order of magnitude to T max /m p = 10 −2 leads to a scenario where both single graviton and graviton pair production processes are equally important, cf. Fig. 6 (right). The total spectrum has now a very characteristic form, where the contribution of the graviton pair production dominates around the peak frequencies. At lower frequencies the contribution from single graviton production is larger than the contribution from graviton pair production. The O λ 2 κ 2 contribution (from single graviton production) is to a good approximation a black body spectrum while the O κ 4 contribution (from graviton pair production) is not.
The results that we have worked out for a scalar model are qualitatively also valid when one considers more general theories. In this case the coupling coefficient λ would have to be replaced with the coupling coefficients that appear in the more general theory. In the future, it would be interesting to work out the details for the SM and BSM theories.
The discussed 2 ↔ 3 and 2 → 2 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 strongest constraints on stochastic GW backgrounds such as the CGMB come from dark radiation constraints: h 2 0 Ω gw 10 −6 . A future detection of the CGMB with Earth based detectors seems to be more challenging [56], see also Ref. [57] for a review on high frequency GW detection prospects in general. However, a possible future detection of the CGMB would bring overwhelming implications for fundamental physics and therefore motivates further effort and research.