Gravitational Waves as a Big Bang Thermometer

There is a guaranteed background of stochastic gravitational waves produced in the thermal plasma in the early universe. Its energy density per logarithmic frequency interval scales with the maximum temperature $T_{\rm max}$ which the primordial plasma attained at the beginning of the standard hot big bang era. It peaks in the microwave range, at around $80\,{\rm GHz}\,[106.75/g_{*s}(T_{\rm max})]^{1/3}$, where $g_{*s}(T_{\rm max})$ is the effective number of entropy degrees of freedom in the primordial plasma at $T_{\rm max}$. We present a state-of-the-art prediction of this Cosmic Gravitational Microwave Background (CGMB) for general models, and carry out calculations for the case of the Standard Model (SM) as well as for several of its extensions. On the side of minimal extensions we consider the Neutrino Minimal SM ($\nu$MSM) and the SM - Axion - Seesaw - Higgs portal inflation model (SMASH), which provide a complete and consistent cosmological history including inflation. As an example of a non-minimal extension of the SM we consider the Minimal Supersymmetric Standard Model (MSSM). Furthermore, we discuss the current upper limits and the prospects to detect the CGMB in laboratory experiments and thus measure the maximum temperature and the effective number of degrees of freedom at the beginning of the hot big bang.


Introduction
Standard hot big bang cosmology provides a successful description of the evolution of the universe back to at least a fraction of a second after its birth, when the primordial plasma was radiationdominated, with temperatures around a few MeV. It nicely explains the Hubble expansion, the cosmic microwave background (CMB) radiation, and the abundance of light elements. But it does not predict the maximum temperature, T max , which the thermal plasma had at the beginning of the radiation-dominated era. It could be arbitrarily high, although there are arguments that the maximum temperature is bounded from above by the Planck scale, T max M P ≡ 1/ √ 8πG 2.435 × 10 18 GeV [1]. At temperatures higher than that quantum gravity effects become very important and we simply do not know what happens in that regime. Nevertheless, it is natural to assume that for T max > M P the gravitons would reach thermal equilibrium and acquire a blackbody spectrum which would decouple at T dec ≈ M P [2]. After decoupling, the blackbody spectrum would simply redshift with the expansion of the Universe, ending up with an effective temperature around 0.9 K [106.75/g * s (T dec )] 1/3 , where g * s (T dec ) is the effective number of entropy degrees of freedom at decoupling.
For T max < M P gravitons are not expected to thermalize, as the Planck-suppressed gravitational interaction rates will remain below the expansion rate of the universe. Nevertheless, out-ofequilibrium gravitational excitations can still be produced from the thermal plasma, and remarkably T max can be probed by gravitational waves (GWs) and bounded by corresponding limits [3,4]. In fact, any plasma in thermal equilibrium emits GWs produced by physical processes ranging from macroscopic hydrodynamic fluctuations to microscopic particle collisions. The magnitude and spectral shape of the corresponding stochastic GW background that is produced during the thermal history of the universe -from the beginning of the thermal radiation dominated epoch after the big bang, at a temperature T max , until the electroweak crossover, at a temperature T ewco 160 GeVhas been calculated in Refs. [3,4]. As the thermal emission always peaks at energies of the order of the temperature, and as the frequency of the emitted waves redshifts in correlation with the temperature, the spectral shape of the ensuing gravitational wave background resembles a bit the blackbody spectrum of photons and neutrinos, its power peaking today in the same microwave domain -that is, for frequencies around ∼ 100 GHz -as the ones for photons and neutrinos. We dub it therefore the Cosmic Gravitational Microwave Background (CGMB), similar to the Cosmic Microwave Background (CMB).
Even though at small frequencies, in the sub-10 kHz range, where all the ongoing and near-future planned GW detectors operate, this stochastic background is many orders of magnitude below the observable level and tiny compared with that from astrophysical and other, more speculative, nonequilibrium sources, the total energy density carried by the microwave part of the spectrum near the peak frequency is non-negligible if the production continues for a long time, that is if T max T ewco . This is due to the fact that, although the thermal rate of production is Planck suppressed, peak emissions at different times add up constructively because of the correlated redshifting of frequency and temperature leads to an approximate linear relation between the total energy emitted in gravitational waves around the peak frequency and the temperature T max . Observing this part directly sets an ambitious but worthwhile goal for future generations of GW detectors, allowing to probe properties of the primordial thermal plasma at the beginning of the hot big bang era, such as its maximum temperature and, as we will see, its effective number of degrees of freedom.
This paper is organized as follows. In Section 2 we determine -based on the work of Refs. [3,4] -the frequency spectrum of the CGMB in a general theory, and subsequently focus on the Standard Model (SM) and several of its extensions. As minimal extensions we choose the Neutrino Minimal Standard Model (νMSM) [5,6] and the SM -Axion -Seesaw -Higgs portal inflation model (SMASH) [7,8]. Both explain neutrino masses and mixing, the non-baryonic dark matter (DM) abundance, the baryon asymmetry of the universe (BAU), and eventually also solve the horizon and flatness problems of the standard hot big bang cosmology. As an example of nonminimal extension of the SM we focus on the Minimal Supersymmetric Standard Model (MSSM) [9][10][11], which is motivated by the Higgs naturalness problem, gauge coupling unification and dark matter. In Section 3 we confront the CGMB predictions with upper limits on the total energy density of any extra relativistic radiation field at the time of big bang nucleosynthesis (BBN) or of decoupling of the CMB photons. In Section 4 we compare the predictions with current limits from direct laboratory searches for GWs and we discuss laboratory experiments which may ultimately probe sub-Planck-scale values of T max . Finally, we summarize our findings and give an outlook for further investigations in Section 5.

GW background from primordial thermal plasma
In this section we exploit the results from Refs. [3,4] concerning the CGMB produced in the primordial thermal plasma at sub-Planckian temperatures, T max < M P . While the former references focused mainly on the SM case, we provide when possible expressions generalized to an arbitrary theory with gauge fields, real scalars and Weyl fermions. The fields are treated as massless, which is a good approximation for temperatures much above the masses of particles in the vacuum. For temperatures below the mass threshold of a given particle, the former decouples from thermal plasma and one can work in an effective theory in which the heavy particle has been integrated out. Therefore, the general results given below can be applied at different temperature ranges when using the appropriate effective theories for the light excitations.
In the following we will start with the general formulae for the production rate of gravitational waves from the primordial plasma and derive expressions for the current energy fraction of gravitational waves per logarithmic frequency interval. Next, we will focus on the predicted spectrum for the SM, to be followed by calculations for three different theories Beyond the SM (BSM): the νMSM, SMASH, and the MSSM.

Production rate of GWs from a general thermal plasma
Here we revise the state-of-the art results for the rate of emission of gravitational waves from a thermal plasma in a generic quantum field theory coupled to gravity. We draw from the results of refs. [3,4]. While the explicit expressions in the former references were tailored for variations of the SM with different numbers of Higgs doublets, generations and colors, we will rewrite the results in a way that facilitates the application to arbitrary quantum field theories with gauge fields, real scalars and Weyl fermions. For our general theory we consider n = 1, . . . N g gauge groups with coupling constants g n ≡ √ 4πα n and Lie algebras of dimension N n spanned by generators T a n , a = 1, . . . , N n . We further assume real scalar fields φ i , i = 1, . . . , N s and Weyl spinors ψ α , α = 1, . . . , N f . The real scalars transform under each group n under a direct sum of irreducible representations r n,î , which can include several copies of the same representation. For each irreducible representation of each gauge group we consider the Dynkin index T n,î defined from the identity Tr r n,î T a n T b n = δ ab T n,î . Analogously, we define fermion representations r n,α with Dynkin indices T n,α . The Dynkin indices of the adjoint representations of the gauge fields themselves will be denoted as T n,Ad .
Regarding the interactions of the fields, it turns out that scalar quartic couplings do not contribute to gravitational wave production at leading order, and thus we will focus on gauge and Yukawa interactions. For the latter we use the convention L ⊃ − i,α,β y i αβ φ i ψ α ψ β + c.c.. (2.1) With the representations of the matter fields defined as above, one may recover the Debye thermal masses of the gauge fields in the plasma from the following expression: m 2 n (T ) = g 2 n (T )T 2 1 3 T n,Ad + 1 6 î T n,î + 1 6 α T n,α ≡ T 2m2 n (T ). (2.2) In the equation above we included a temperature dependence of the couplings g n (T ), arising from choosing a renormalization scale proportional to the temperature. This is expected to provide optimal accuracy for the computations of thermal effects, as they involve excitations whose typical momenta are of the order of T . This choice of renormalization scale implies that the dimensionless quantitym 2 n (T ) inherits a logarithmic temperature dependence which has been explicitly indicated.
Within the gauge interactions, we will assume that hypercharge is the weakest. This has an impact in the production of gravitational waves with low frequencies, which as will be seen later is related to the plasma's shear viscosity [3], which is known to be dominated by the effect of the weakest gauge interaction [12]. Due to this, following the former reference it is convenient to define a number N species given by one half of the sum over the hypercharge Dynkin indices of the real scalar and Weyl fermion representations. Analogously one can define N leptons as one-half of the squared hypercharges of the Weyl fermions that interact with no other SM gauge group than hypercharge. Assigning k = 1 to the hypercharge group, one has We note that for the expression for N leptons it was assumed that the only fields that interact exclusively with hypercharge are fermions. We expect that scalar fields with similar properties will also contribute to N leptons ; however, the estimates of transport coefficients in Ref. [12] did not account for such fields. For the MSSM studies in Section 2.5 we will assume a contribution to N leptons coming from the supersymmetric partners of the right-handed leptons, obtained by adding to N leptons in Eq. (2.3) the analogous sum over representations of real scalar fields. In a homogeneous and isotropic universe, with scale factor a and Hubble expansion rate H = a/a, the energy density ρ CGMB carried by the CGMB, which was generated in a thermal plasma with temperature T , evolves in cosmic time t as The former equation assumes a very small energy density of gravitational waves, so that one can neglect the backreaction contribution from gravitational excitations annihilating or decaying back into the plasma. As emphasized in the introduction, this is expected to be a good approximation for temperatures below the Planck scale. For momenta lower than the temperature, the dimensionless source termη T, k T can be understood from long-range hydrodynamic fluctuations, while for momenta comparable or greater than the temperature,η T, k T is dominated by contributions from quasi-particle excitations in the plasma [3,4]. While the hydrodynamic contribution is known to leading-log order in the gauge couplings, recently the quasi-particle contribution has been estimated to full leading order in the gauge and Yukawa couplings [4]. The results for temperatures above the critical T of the electroweak phase transition can be written as: In the equations above we defined a dimensionless momentumk = k/T and introduced a coefficient η and functionsη HTL (T,k), η gg (k), η sg (k), η f g (k), η sf (k) which will be described next.
First, the hydrodynamic contribution for k α 2 1 T coincides with the shear-viscosity of the plasma divided by T 3 [3],η The shear viscosity is inversely proportional to a scattering cross section and therefore large for a plasma in which there are some weakly interacting particles. Under our assumption that hypercharge is the weakest gauge force, right-handed leptons (or additional fields only charged under U (1) Y ) are the most weakly interacting degrees of freedom, changing their momenta only through reactions mediated by hypercharge gauge fields above the electroweak crossover. The results of Ref. [12] give then the following value for the coefficientη in Eq. (2.5): where ζ is Riemann's zeta function. For k max {α 2 n }T , GWs are dominantly produced via microscopic particle scatterings, despite the fact that their rates are suppressed by the coupling strengths responsible for the interactions and a Boltzmann factor e −k/T , which takes into account that the energy carried away by the graviton must be extracted from thermal fluctuations. These contributions were first estimated at leading-log accuracy in Ref. [3] and then calculated at full leading order in Ref. [4]. In the latter calculation, the production rate is obtained from the imaginary part of the two-point correlator of the stress-energy momentum tensor at two-loop order. Some of the loop integrals involved turn out to be infrared divergent when treating the fields as massless. A resummation of the effects of the thermal masses resolves the singularities, and Ref. [4] implemented this in the following manner: a contribution containing the divergence was added and subtracted; the subtracted part was used to define infrared-finite loop integrals, while the added piece was rendered finite by implementing the resummation of the thermal masses. In this way, the thermal resummation is only performed in the region of phase space near the singularities, but the procedure guarantees full leading-order accuracy. The functionη HTL (T,k) in Eq. (2.5) corresponds to the regulated divergence (with "HTL" alluding to the hard thermal loop resummation [13]) and is given bŷ .

(2.8)
This is to be contrasted with the initial leading-log estimate ofη,η LL , that was computed in Ref. [3]: .
Finally, the remaining functions in Eq. (2.5) correspond to the infrared-finite two-loop integrals mentioned before, for diagrams involving only gauge fields (η gg (k)), scalars and gauge fields (η sg (k)), fermions and gauge fields (η f g (k)), and scalars and fermions (η sf (k)). The infrared subtraction has to be performed in the gg, sg and f g sectors. Explicit expressions for the loop functions are given in Appendix A (see Eqs. (A.1),(A.2)). From the latter equations and from Eqs. (2.5), (2.8) it is clear that, as advertised earlier, all contributions to the production rate are suppressed by powers of the gauge and Yukawa couplings, as well as Boltzmann factors.

Current stochastic GW background from a primordial thermal plasma
In this section we relate the thermal production rate discussed above with the stochastic background of gravitational waves in the present universe. At every value of the temperature, the Boltzmann suppression factor ensures a peak emission for momenta of the order of the temperature. The expansion of the Universe redshifts temperature and momenta by approximately the same amount, and as a consequence of this the peak emission at a given temperature overlaps with the redshifted peak emissions of the previous history of the universe. Thus, the energy density of gravitational waves at the peak frequency is sensitive to the entire history of the hot primordial plasma, so that the weak Planck-suppressed production rates can be partly compensated. Similarly to the case of the CMB, the peak of the thermal spectrum lies currently in the microwave region -simply because the current CMB temperature T 0 ≈ 2.73 K is associated with frequencies in the 100 GHz regimeleading to the CGMB. Given that (∂ t + 3H(t))s = 0, where s is the entropy density, the factor 4H(t) in Eq. (2.4) can be taken care of by normalizing ρ CGMB by s 4/3 . Subsequently, the equation can be integrated from the time t hbb , when the hot big bang era starts with a temperature T max (< M P ), to the time t ewco , when the electroweak crossover takes place, by assuming that at t hBB there were no (thermally produced) gravitational waves present: Here we have used that time and temperature are related as [14] dT are the effective degrees of freedom of the energy density, ρ(T ), the entropy density, s(T ), and the heat capacity, c(T ), respectively. The above result can be written in the form from which we can read off the spectrum of the GW energy density fraction per logarithmic wave number interval at the time of the electroweak crossover, . We have obtained this by taking into account that momenta redshift as (2.15) and expressed the momentum space integral in (2.13) in terms of k ewco . The spectrum of the GW energy density fraction per logarithmic wave number interval at the present time is then 3.931 is the number of effective degrees of freedom of the entropy density after neutrino decoupling and Ω γ ≡ ρ and For order of magnitude estimates and an analytic understanding of the result, one may exploit the fact that the temperature dependence of the effective degrees of freedom of the energy density, entropy density, and heat capacity approximately agree amongst each other above the electroweak crossover, g * ρ (T ) ≈ g * s (T ) ≈ g * c (T ), for T T ewco , see Appendix B. Therefore, Eq. (2.19) can be approximated by Furthermore, assuming that g * s andη are almost independent of temperature above the electroweak cross-over, up to possible steps when new BSM degrees of freedom get relativistic, and exploiting the fact that the T dependence ofη(T,k) is only logarithmic (from the running of the coupling constants with temperature), one may obtain an approximate analytic expression by ignoring the temperature dependence of the integrand. Doing so the energy density of the CGMB per logarithmic frequency interval, valid for T max T ewco , can be approximated as: Its overall magnitude scales approximately linearly with the maximum temperature of the hot big bang. Therefore, it can play the role of a hot big bang thermometer. Moreover, a measurement of the peak frequency of Ω CGMB provides a measurement of the relativistic degrees of freedom at T max . In fact, the peak frequency is, according to Eq. (2.21), approximately determined by the peak ofk 3η (T max ,k), which in turn can be estimated from its leading-log behaviour,k where W is the Lambert W function, leading to Up to the factor g * s(fin) g * s(Tmax)

1/3
, it coincides with the peak frequency, of the present energy fraction of the CMB per logarithmic frequency, In fact, around the peak and to leading-log accuracy, the spectral form of the CGMB resembles the one of the CMB with an effective temperature Direct detection bounds or projected sensitivities are often expressed in terms of a characteristic dimensionless GW amplitude h c (f ), which is related to the cosmic energy density fraction of stochastic GWs per logarithmic frequency interval as (we use the conventions of Ref. [15]) For the CGMB, its overall magnitude scales approximately linearly with the square root of the maximum temperature of the hot big bang. Moreover, its peak frequency is approximately determined by the peak ofkη(T max ,k), which in turn can be estimated from its leading-log behaviour, become maximal in an arbitrary theory, taking into account the complete leading-order expression fork 3η (T max ,k) andkη(T max ,k), can be obtained as follows. First, one can determine numerical approximations for each loop function appearing in Eq. (2.5), multiplied byk 3 ork, around its critical point, using a Taylor expansion of order 2 with coefficients determined numerically. While the functionsη HTL ,η sf have maxima, the loop functionsη gg ,η sg ,η f g have minima, and all the critical points have similar values ofk. As the latter three functions give only relatively small corrections toη HTL , it turns out that the location of the maxima ofk 3η andkη are close to all the maxima and minima associated with the different loop functions, and hence doing Taylor expansions around the critical points allows for accurate determinations of the peak frequencies.
Using the notation of Sec. 2.1, we first define shorthands for the coefficients of the loop functions appearing in Eq. (2.5) and evaluated at T = T max : (2.33) From this one may estimate the values ofk Ω peak ,k hc peak as: (2.35) The aforementioned quantities depend onm 2 n , and we estimated the dependency with a numerical fit. The sum in (xy) goes over the pairs (gg), (sg), (f g), (sf ). The quantities f Ω/hc xy are numbers related to the second derivatives of the functionsk 3 η xy (k) andk η xy (k) at their respective extrema, while thek Ω/hc xy are the locations of these extrema. We have computed the above numbers and functions as:   To obtain the peak frequencies, as follows from Eq. (2.21) one simply has to use: where the values ofk Ω peak (T max ),k hc peak (T max ) can be computed from Eqs. (2.34), (2.36) and (2.37). In Eq. (2.38) we normalized the values ofk Ω peak (T max ) andk hc peak (T max ) by their values in the SM at T max = 10 16 GeV.
The above formulae, obtained by exploiting the approximation (2.21), reproduce with an accuracy better than 1% the peak frequencies of the different models considered (see later in Table 2) computed according to Eq. (2.19) from the integration over temperature of the full production rate. In the models to be analyzed in the following sections, we find that the differences in values ofk Ω/hc peak are smaller than the variation of g * s (T max ). For example,k Ω peak only changes by 5% between the SM and the MSSM, whilek hc peak changes by 15%, and g * s (T max ) varies by 30%. Thus the effect of g * s (T max ) dominates in Eq. (2.38), which predicts that peak frequencies decrease with growing values of g * s (T max ). The reduced variability ofk Ω/hc peak can be understood from (2.34) by noting that the model dependence enters through the values of the Debye massesm 2 n appearing in the HTL contributions, and through the coefficients c xy . In general, the HTL contributions dominate both the numerator and denominator, so that there is an approximate cancellation of the leading model dependence throughm 2 n in Eq. (2.34). Similarly, we can also give approximate formulae for the values of Ω CGMB and h CGMB c at their respective peak frequencies: (2.40) In the formulae above,m 2 n are meant to be evaluated at T = T max , and aside from the functions and constants of Eqs. (2.36) and (2.37), we introduced the quantities g xy , which correspond to the peak values of k 3 η xy (k) andk η xy (k). These functions and constants are given next: Again, we find that the approximate expressions (2.39) and (2.40) reproduce with an accuracy better than 3% (for Ω CGMB ) and 1% (for h CGMB c ) the results for the different models (see later in Table 2) to be analyzed next. Aside from the model-dependence coming from g * s (T max ), we note that, as opposed to the case ofk Ω/hc peak , there is no approximate cancellation of the dependence on the values ofm 2 n . Thus, we anticipate more variability of the peak values of h 2 Ω CGMB and h CGMB c across models for a fixed T max . Nevertheless, in weakly coupled extensions of the SM it is expected that the variations inm 2 n and c xy will be subleading with respect to changes in g * s (T max ). Under this assumption, the leading model dependence in the peak frequencies and in the peak values of h 2 Ω CGMB and h CGMB c can be captured by g * s (T max ) alone, so that for a general model one has: as implied by Eqs. (2.39), (2.40) when ignoring changes inm 2 n , c xy . Throughout our explorations of minimal and non-minimal extensions of the SM in the next subsections, we find that Eq. (2.43) is satisfied with better than 15% accuracy, while the relation for h CGMB c in Eq. (2.44) works with better than 30% accuracy. We find that by simply comparing with the SM predictions, a measurement of both the peak frequency and the maximum value of h CGMB c could be used to estimate the hot big bang temperature T max and the value of g 1/3 * s (T max ) up to systematic effects below 15% for g 1/3 * s (T max ) and below 40% for T max . These numbers correspond to the MSSM, while for example in SMASH the accuracies would be around 1% for g 1/3 * s (T max ) and 5% for T max .
We also note that Eq. (2.44) predicts that, in weakly coupled theories, the SM is expected to give the highest power of thermally produced gravitational waves. Naively, in SM extensions one would expect an increase in the instantaneous rate of production of gravitational waves, as the coefficients of the loop functions in Eq. (2.5) will receive additional contributions. However, the presence of additional degrees of freedom also implies that the waves produced at early times will suffer a larger red-shifting, as the latter is proportional to [g * s (T emission )/g * s (fin)] 4/3 (note that g * s (T ) −4/3 is indeed present in Eq. (2.19)). As long as the models are weakly coupled, the effect of redshifting will be dominant, and this is captured by the approximate relations in Eq. (2.44).
To end this section, let us recall that the previous results are valid for T max < M P , for which the gravitons never reach thermal equilibrium and one may ignore the backreaction from graviton annihilations and decays in the rates of Eq. (2.4) and (2.5). For gravitons that were in thermal equilibrium for temperatures above the Planck scale, the prediction for the relative abundance of gravitational waves is that of a blackbody spectrum with an effective temperature obtained by redshifting the decoupling temperature M P by the expansion of the universe between the decoupling time and the present [2]. This follows simply from noting that after decoupling, gravitons would stop interacting and start to propagate freely, with their momenta redshifting due to the expansion. As a consequence of this one has For the equilibrated spectrum the peak frequencies are given by Eqs. (2.24) and (2.32), with T max in these expressions replaced by the decoupling temperature M P , without any further modeldependence than the one coming from g * s . The maxima of the spectra however scale differently with g * s than in Eqs. (2.44). For the equilibrated spectrum one has (2.46) Figure 1: Function (k/T ) 3η SM (T, k/T ) determining the background of stochastic gravitational waves produced in the thermal SM plasma, showing the hydrodynamic contributions (straighter solid lines for smallish k/T ), the microscopic contributions at full leading order (curved solid lines for higher k/T ), and their leading-log approximations (dashed lines). The lines are colored for the scales in which the calculations can be trusted, i.e. k < α 2 1 T for the hydrodynamic contributions, and k > m 3 (T ) for the microscopic ones. From top to bottom, the temperatures correspond to T = 10 3 GeV (black), T = 10 8 GeV (red), T = 10 13 GeV (green) and T = M P (violet). The gauge and Yukawa couplings were evaluated at the renormalizationμ = 2πT using the respective two-loop renormalization group equations.

CGMB in the SM
In the SM we assign n = 1, 2, 3 to the gauge groups U(1) Y , SU(2) L , SU(3). With the SM matter content one has (2.47) With this one can fixη as well as theη HTL contribution of Eqs. (2.5) and (2.8). Computing as well the coefficients of the loop functions in terms of the representations and couplings in the SM leads to: In the previous equations we omitted for simplicity the logarithmic T -dependence of the couplings g i , y i and the rescaled Debye massesm i . We have also ignored the Yukawa couplings of the lightest fermions.  [16], supplemented with the g 3 -dependent three-loop contribution to the running of g 3 [17], evaluated at a temperature-dependent renormalization scale. The couplings were fixed at low scales using values m t = 172.9 GeV, m h = 125.10 GeV for the physical top and Higgs masses. For the determination of y t from the top mass we used oneloop electroweak and three-loop QCD threshold corrections [18][19][20], while for computing the Higgs couplings we used the full two-loop effective potential plus appropriate momentum corrections, as in Ref. [21]. We made different choices of the renormalization scale, µ = κ2πT with κ = 1/2, 1, 2, in order to estimate theoretical uncertainties. Figs. 1 and 2 illustrate how the leading-log result captures the leading-order result quantitatively for k/T 10, whereas it overestimates the latter by a factor around two in the phenomenologically most interesting region k/T ∼ 4. The peak positions ofk 3η SM (T,k), on the other hand, are shifted slightly, by less then 10 % from the generic valuê k Ω CGMB The bands show the uncertainties coming from changing the RG scale within a factor of 2 and from shifting the unknown parameter q c of the three-loop QCD corrections to the pressure in Ref. [22] between −3000 GeV and 3000 GeV. The calculation is expected to lose accuracy near the electroweak phase transition around 160 GeV.
For the computation of the spectrum of thermally produced gravitational waves one has to use Eq. (2.19) and carry out the numerical integration. This requires knowledge of the functions g * ρ (T ), g * s (T ) and g * c (T ). As reviewed in Appendix B, all these quantities can be derived from the thermal corrections to the effective potential -which correspond to minus the pressure of the thermal plasma-and the use of thermodynamical relations. In our calculations we use the full one-loop contributions to the thermal potential, supplemented with three-loop QCD contributions [22]. As we consider gravitational wave production before the electroweak crossover, we can use perturbative results; for lower temperatures one requires more sophisticated techniques [23]. The values used for g SM * ρ (T ), g SM * s (T ) and g SM * c (T ) are shown in Fig. 3. The resulting spectrum of thermally produced gravitational waves is shown in Fig. 4 for different values of the maximum temperature, together with the predicted sensitivities of upcoming gravitational wave experiments like the Big Bang Observer (BBO) [24], the Cosmic Explorer (CE) [25], the Deci-hertz Interferometer Gravitational Wave Observatory (DECIGO) [26], the Einstein Telescope (ET) [27], and LISA [28]. The sensitivity projections were taken from [29]; for ultimate DECIGO we use the curve in Ref. [30] based on Ref. [31].
In Fig. 5 we show the spectra rescaled by T max and compare these results with the analytic approximation (2.21), which predicts a value of h 2 Ω CGMB /T max independent of T max aside from variations in g * s (T max ). The figure shows that the analytic prediction gives results with an accuracy better than 3% near the peak for T max 10 5 GeV. 1 Within this uncertainty, in accordance with the expectation from Eq. (2.21) the absolute value of Ω CGMB (f ) scales approximately linearly with T max /M P . Therefore, a measurement of it determines the maximum temperature of the hot big bang. The peaks in the spectra, for different T max , occur around 80 GHz, less than 10 % higher than the generic estimate (2.24) based on the analytic approximation (2.21) and the leading-log result forη(T,k), while they are reproduced with an accuracy of the order of 3% or better (1% or better for h CGMB c ) by the formulae (2.38) and (2.34). To end this section, let us note that the theoretical uncertainty of the above results for h 2 Ω CGMB in the SM is of the order of 0.1%. This has been estimated by considering the effect of varying the renormalization scale by a factor of 2, and by considering values between -3000 GeV and 3000 GeV of the unknown parameter q c appearing in the three-loop contributions to the QCD pressure of Ref. [22]. Note that the final uncertainty is one order of magnitude lower than the maximal theoretical uncertainties found forη; this is due to cancellations between the variations ofη and the effective numbers of degrees of freedom.

CGMB in minimal BSM models explaining neutrino masses, DM, and the BAU
So far, our predictions were based on the assumption that the SM is valid up to the Planck scale, and the value of the temperature T max was left unspecified. However, there is a strong case for BSM physics. It is definitely required to explain neutrino masses and mixing, the origin of the nonbaryonic DM, and the BAU. Therefore, we consider now two minimalistic extensions of the SM which solve also these problems. In addition to the latter issues, these models also accommodate realizations of the inflation mechanism, which can address the flatness and horizon problems asso- ciated with the observed lack of curvature and the striking homogeneity of the universe. As such, they have all the necessary ingredients to explain the cosmic history of the universe from inflation until the present time; in particular, this means that they give rise to concrete predictions of T max as a result of the post-inflationary reheating dynamics. In turn, this gives refined predictions for the spectrum of gravitational waves originated from the thermal plasma.
The νMSM [5,6] extends the SM by three right-handed SM singlet neutrinos, which have a GeV scale Majorana neutrino mass and mix with the three active left-handed neutrinos via Yukawa interactions with the SM Higgs field. This model may be valid up to the Planck scale. Neutrino masses and mixing are generated by the type-I seesaw mechanism [32][33][34][35]. DM is comprised by a keV-scale neutrino mass eigenstate, and the BAU is produced by a low-scale leptogenesis mechanism involving neutrino oscillations [36]. Chaotic inflation can be provided by the Higgs field when allowing for a non-minimal gravitational coupling  [40,41], yet arguments against this have been given e.g. in [42][43][44]. In any case, one may have to pay the prize of uncertain predictions due to unknown non-perturbative corrections to the tree level results. Nevertheless, ignoring this caveat, the tensor-to-scalar ratio and the maximum temperature of the universe in the νMSM after reheating from Higgs inflation have been determined as [38] r 0.0034, These temperatures are much below the absolute upper bound following from the CMB constraint on the tensor-to-scalar ratio, r = P T /P S < 0.058, 2 , and the unphysical assumption of instantaneous and maximally efficient reheating to a radiation dominated universe, cf. Appendix C. The thermal plasma of the νMSM differs from the thermal plasma of the SM only slightly -from the subleading effects of the Yukawas of the singlet neutrinos, that contribute to the η sf term in Eq. (2.5) -and therefore the rate of production can be approximated with that in the SM, Eq. (2.48). In regards to the calculation of the present day spectrum using Eq. (2.19), one has to use the values for g * ρ , g * s , g * c appropriate for the νMSM. We assume that the singlet neutrinos remain in thermal equilibrium above the electroweak crossover, so that the values of the effective degrees of freedom can be obtained from those of the SM by adding 3 units. An alternative minimal extension of the SM explaining the origin of DM and the BAU is SMASH [7,8]. A SM singlet complex scalar field σ, which features a spontaneously broken global U (1) PQ Peccei-Quinn (PQ) symmetry [46], and a vector-like coloured Dirac fermion Q are added to the field content of the νMSM. Exploiting the PQ mechanism, this model solves the strong CP problem. DM is comprised by the axion [47][48][49] -the pseudo Nambu-Goldstone boson of the U (1) PQ breaking [50,51] -provided that the PQ breaking scale is in the range 1.3×10 9 v σ /GeV 2.2 × 10 11 [52]. The right-handed neutrinos get their Majorana masses also from spontaneous PQ symmetry breaking. The generation of the BAU proceeds via high-scale thermal leptogenesis [53]. Finally, inflation can be accommodated for perturbative values of the non-minimal gravitational couplings [7,8]. Allowing for a non-minimal coupling ξ σ of the PQ field to the Ricci scalar, S ⊃ − d 4 x √ −g ξ σ σ * σ R, a mixture of the modulus of the complex PQ field, ρ = √ 2 |σ|, with h, the neutral component of the SM Higgs doublet in the unitary gauge, is a viable inflaton candidate. Fitting the inflationary predictions to the observed fluctuations in the CMB relates the size of the non-minimal coupling and the quartic coupling; for the latest Planck data [45] this gives [30] 7 × 10 −3 ξ σ 4 × 10 4 λ σ 1. (2.51) The above window was obtained after ensuring a consistent post-inflationary history in which Planck's CMB pivot scale was matched to the appropriate mode during inflation. The lower bound, ξ σ 7 × 10 −3 , arises from taking into account the upper limit on the tensor-to-scalar ratio, r < 0.058, while the upper bound, ξ σ 1, arises from perturbativity and unitarity requirements 3 . It corresponds to a lower limit on the tensor-to-scalar ratio, r 4 × 10 −3 . As a consequence, the quartic coupling should be in the range 7 × 10 −13 λ σ 5 × 10 −10 . The initial conditions for the standard hot big bang cosmology following inflation, non-perturbative preheating and perturbative reheating can be predicted from first principles in SMASH. The maximum temperature of the thermalized SMASH plasma after reheating is obtained as [8] 8 × 10 9 GeV T SMASH max 2 × 10 10 GeV . (2.52) Again, this is significantly below the upper bound on T max following from the assumption of instant reheating and the CMB constraint r < 0.058, which in this case gives In order to calculate the CGMB in SMASH we can use the general expressions of Section 2.1. This requires knowing the BSM Yukawa couplings in SMASH. Stability in the σ direction demands small couplings for the RH neutrinos [8], whose effect will be ignored as in the νMSM; this leaves the Yukawa couplings of the exotic vector quark. Assuming a small interaction with the down quarks (we consider the SMASH realization with hypercharge 1/3 for Q, in which such mixing allows the Qs to decay before nucleosynthesis), one has L ⊃ y Q σQP L Q + c.c..
,k α 2 1 , η HTL,SMASH (T,k) + (3g 2 2 + 12g 2 3 )η gg (k) For the computation of g * ρ , g * s , g * c in SMASH we proceed as in Ref. [30]. In order to reliably follow the change in degrees of freedom across the PQ phase transition, one has to use an improved daisy resummation of thermal effects compatible with thermal decoupling. At low temperatures, the SMASH theory is matched to the SM plus the real part of σ and the nearly massless axion; we include again three-loop QCD corrections plus corrections from the loss of chemical equilibrium of the axion due to its feeble interactions, which imply that the axion population has a different effective temperature than the rest of the plasma. For details, see Ref. [30]; a summary is given in Appendix B. Figure 6 shows results for g SMASH * ρ (T ), g SMASH * s (T ) and g SMASH * c (T ) for two benchmark points with r = 0.0037 and r = 0.048, taken from Ref. [30].
With the previous results forη, g * ρ , g * s and g * c , and taking into account the ranges of T max in Eqs. (2.49) and (2.52), one can use Eq. (2.19) to calculate the predictions for the CGMB in the νMSM and SMASH. As the results forη assume massless fields, in the SMASH case we use (T ) (red) for benchmark points with r = 0.0037 (left) and r = 0.048 (right). The bands show the uncertainties coming from changing the RG scale within a factor of 2 and from shifting the unknown parameter q c of the three-loop QCD corrections to the pressure in Ref. [22] between −3000 GeV and 3000 GeV. The calculation is expected to lose accuracy near the electroweak phase transition around 160 GeV.
Eq. (2.56) at high temperatures, and the SM result of Eq. (2.48) below the temperature at which the axion field decouples. The latter was estimated as in Ref. [30] by finding the temperature near the critical temperature of the PQ phase transition at which the trace of the stress-energy momentum tensor has a local maximum. The results are shown in Fig. 7, together with the inflationary Cosmic Gravitational Wave Background (iCGWB) due to the tensor modes generated by quantum fluctuations during inflation. For frequencies f 10 −16 Hz, which re-entered the horizon during radiation domination, the iCGWB can be calculated as [23], In the equation above P T (f ) is the power spectrum of gravitational waves generated during inflation expressed in terms of the present frequency, where a inf , H inf are the scale factor and Hubble constant during inflation, which we have computed assuming non-critical inflation for the νMSM and without resorting to the usual slow-roll approximation, but rather by numerically solving the equation of motion for the inflationary background as a function of the number of efolds [8,54]. T hc (f ) in Eq. (2.57) is the temperature at which the mode corresponding to the frequency f re-entered the horizon during radiation domination. It can be obtained by solving [23] T hc (f ) = 10 8 GeV × f 1.2 Hz × g * s (fin) 3.91 The iCGWB has an upper cutoff corresponding to frequencies that never exited the horizon during inflation. We have approximated this by a sharp feature, yet for frequencies near this threshold our calculations beyond the slow-roll approximation already show a drop in the power spectrum, as can be be seen in Fig. 7. 4 Note that, for both SMASH and the non-critical νMSM, the iCGGWB does not overlap with the peak in the CGMB, so that both sources of gravitational waves -inflationary perturbations and thermal processes -become distinguishable if experiments reach the appropriate sensitivity.
In the case of SMASH, the theoretical uncertainty in the calculations of h 2 Ω CGMB /T max is of the order of 2%. As before, this quantity corresponds to the variations of the results under changes of the renormalization scale and the unknown three-loop contributions to the QCD pressure.
In our calculations we find that at sufficiently high temperaturesk 3η SMASH peaks neark Ω peak ∼ 4.2, and within 1% range of the SM results at the same temperatures. Similarly we getk hc peak ∼ 2.1-2.2, within 1.5% of the corresponding SM results. This means that Eq. (2.44) is satisfied with better than ∼ 1% accuracy. On the other hand, we find that the relation for h c in Eq. (2.44) is satisfied with better than 4% accuracy. Under the assumption of a SMASH plasma and given a measurement of the peak frequency and the maximum value of h CGM B c , then Eqs. (2.43) and (2.44) would allow to estimate g 1/3 * s (T max ) and T max within errors around 1% and 5%, respectively. Some details of the peak frequencies and amplitudes for several SMASH benchmark points are given in Table 2.

CGMB in the MSSM
The MSSM goes beyond the SM by adding a fermionic (scalar) superpartners for all the SM bosons (fermions). Additionally, the model contains an extra scalar Higgs doublet with its corresponding fermionic partners.
Given the matter content of the MSSM, and assuming that the scalar partners of the righthanded leptons contribute to N leptons in analogous manner to the usual leptons (i.e. adding the contribution 1/2 î :T n,î =0,n>1 T 1,î to Eq. (2.60) Computing the coefficients of the loop functions in Eq. (2.5) leads to: Note how in the sg and f g contributions the coefficients in front of the gauge couplings squared are larger than in the previous models, due to the extra matter fields charged under the gauge interactions. Additionally, one has gauge coupling contributions in the coefficient of the loop function η sf , as a consequence of the fact that supersymmetry implies a relation between the Yukawa couplings of the gauge superpartners and the usual gauge couplings. Analogously, the coefficients of the usual Yukawa couplings are larger than before because supersymmetry relates the usual Yukawa couplings to those of additional interactions involving scalar superpartners. For our estimates of gravitational wave spectra in the MSSM, we have used the naive value of the effective number of relativistic degrees of freedom, g * ρ,MSSM (T ) ≈ g * s,MSSM (T ) ≈ g * c,MSSM (T ) ≈ 228.75. (2.62) The reason for this simplification is that we lack knowledge of the QCD corrections to the pressure coming from scalar superpartners. For our numerical estimates we consider a simple scenario in which the dimensionful parameters in the MSSM that are not present in the SM are assumed to lie around the 2 TeV scale. We further assume that the lightest neutral Higgs state is SM-like, which can be realized with a small neutral Higgs mixing angle α (we take cot α = 10) and a heavy

Dark radiation constraint on the CGMB
The CGMB acts as an additional dark radiation field in the universe. Any observable capable of probing the expansion rate of the universe, and hence its energy density, has therefore the potential ability to constrain the CGMB energy density ρ CGMB present in that moment. BBN and the process of photon decoupling of the CMB yield a very precise measurement of H, when the universe had a temperature of T BBN ∼ 0. Since the energy density in the CGMB must satisfy ρ CGMB (T ) ≤ ∆ρ rad (T ), one finds a constraint on the CGMB energy density redshifted to today in terms of the number of extra neutrino species, where ρ (0) c = 3H 2 0 M 2 P , and we have used g * s (T = MeV) ≈ 10.75 and g * s (fin) 3.931. This bound corresponds roughly to a direct bound on the CGMB energy fraction per logarithmic frequency interval, because it has a large width of order the peak frequency itself. The latest BBN constraints on ∆N ν can be found in Ref. [59]. The 4 He alone is not very constraining due to degeneracies with the baryon-to-photon ratio η B , so that the best constraints come from combining BBN measurements of 4 He and deuterium abundances with CMB results. In this case Ref. [59] finds ∆N ν < 0.3 at 95% implying, from Eq. (3.2),  Table 1: Values of the maximal temperature allowed by the dark radiation constraints of Eqs. (3.6)-(3.5), as well as the bound assuming ∆N ν equal to the theoretical uncertainty of 10 −3 . In the case of SMASH the numbers between brackets reflect the change in the last significant digit coming from choosing benchmark scenarios with different values of the tensor-to-scalar ratio. The results for trans-Planckian temperatures are not physically meaningful, as for T > M P one expects early-time equilibration of gravitons and a spectrum as in Eq. (2.45). We list the trans-Planckian temperatures simply to illustrate the reach of current dark radiation bounds.
A similar bound is obtained from other inferences from CMB [60][61][62]. In particular, the analysis of Ref. [62] uses Planck data, together with CMB lensing, baryon acoustic oscillations and also deuterium abundances, and finds a constraint that goes down to Not surprisingly, this is comparable to what is obtained from the BBN analysis in Ref. [59], which also uses CMB data to pin down the baryon to photon ratio η B . However, Ref. [62] only analyses adiabatic initial conditions. From the results of Refs. [60,61], one can infer that there is a gain when imposing homogeneous initial conditions, due to the breaking of degeneracies with neutrino parameters [15]. This has been confirmed by Ref. [63], which under the hypothesis of GW with homogeneous initial conditions finds Finally, it should be noted that the current theoretical uncertainty of ∆N ν is of the order of 10 −3 [4]. If experiments were to reach this level of precision, one would obtain an upper bound of h 2 ρ (0) CGMB /ρ (0) c < 5.7 × 10 −9 . We have turned the above limits into upper bounds on T max for the SM, the νMSM, SMASH and the MSSM, cf. Table 1. The observational limits of Eqs. (3.5) and (3.6) give bounds of the order of 10 19 GeV. These correspond to temperatures above the Planck scale, for which the gravitons can be expected to enter thermal equilibrium and the calculations based on Eq. (2.5) cannot be applied. Thus the previous dark radiation bounds cannot reliably constrain T max . For trans-Planckian temperatures one has to use the equilibrium form (2.45) for the CGMB spectrum; integrating over the frequency so as to obtain the total energy fraction gives Intriguingly, this just about saturates the current dark radiation bound obtained assuming homogeneous initial conditions, Eq. (3.6). Note that, in the case of early time equilibration of gravitational waves, one expects in fact homogeneity, and thus the the relevant dark radiation bound is indeed given by Eq. (3.6) instead of (3.5). Thus the current dark radiation bound is just on top of the value that corresponds to the contribution from gravitational waves that were in equilibrium at  Table 2: Values of peak frequencies and peak power spectra for the gravitational waves produced from the thermal plasma, for different models and maximum temperatures related to dark radiation and inflationary bounds, or to direct estimates of the reheating temperature. For a given model, the upper row corresponds to temperatures above the Planck mass, for which gravitons are expected to reach thermal equilibrium at early times, leading to an energy fraction around the dark radiation bound of Eq. (3.6). The temperature in the second row is the upper bound corresponding to a hypothetical constraint ∆N ν = 10 −3 . In the third row one has the temperature bound of Eq. (C.4), which applies under the assumption of slow-roll inflation. If present, the range of T max in the fourth row corresponds to direct estimates of the reheating temperature.
early times. Taking the significant digits of the bound of Eq. (3.6) seriously, then the result of Eq. (3.7) would imply that current dark radiation bounds are compatible with a CGMB with early time equilibrium in an extension of the SM in which g * s is augmented by a few degrees of freedom (∆g * s > 2.8 taking the naive value g * s,SM = 106.75, ∆g * s > 4 when including additional radiative corrections as summarized in Appendix B (see Fig. 3). The next generation of CMB experiments is expected to improve the sensitivity on ∆N ν by one order of magnitude. Correspondingly, the upper bound on T max may decrease by a factor of ten and thus reach the reduced Planck scale in the next decade. If future experiments were to reach the theoretical uncertainty ∆N ν ∼ 10 −3 , then one would probe T max at scales of the order of 10 17 GeV, as was already emphasized in Ref. [4]. Note that T max bounds increase for models with more degrees of freedom, as expected from the scaling of Eq. (2.44). More details for the peak frequencies and values of Ω CGMB , h CGMB c for the maximal temperatures that follow from the dark radiation bounds are given in Table 2.

CMB Rayleigh-Jeans tail constraint on the CGMB
In the presence of magnetic fields, GWs are converted into electromagnetic waves (EMWs) and vice versa. This is called the (inverse) Gertsenshtein effect [64][65][66][67][68]. Recently, it has been shown that this conversion results in a distortion of the CMB, which can act therefore as a detector for MHz to GHz GWs generated before reionization [69]. The measurements of the radio telescope EDGES have been turned into the bound These upper bounds are displayed in Fig. 8 as green exclusion regions. Future advances in radio astronomy and a better knowledge of cosmic magnetic fields are required in order that this method can get competitive with the dark radiation constraint.

Laboratory searches for the CGMB
In this section we will discuss current constraints on the CGMB from direct experimental GW searches in the laboratory and future possibilities to search for a stochastic GW background in the frequency range around the peak.

Current direct bounds from GW experiments
Current large-size ground-based laser-interferometric GW detectors, such as GEO, KAGRA, LIGO, and VIRGO [70][71][72][73], are sensitive in the frequency range from about 10 Hz to 10 kHz. Their technology is not necessarily ideal for studying very-high-frequency (VHF: 100 kHz − 1 THz) and ultra-high-frequency (UHF: above 1 THz) GWs. Several other small-size experiments have performed pioneering searches for stochastic GWs in the VHF and UHF range and put corresponding upper bounds which are confronted in Fig. 8 to the CGMB prediction: • A cavity/waveguide prototype experiment searched for polarization changes of electromagnetic waves, which are predicted to rotate under an incoming GW [74]. It provided an upper limit, h c < 1.4 × 10 −10 , on the characteristic amplitude of stochastic GWs at 100 MHz.
• Two laser interferometers with 0.75 m long arms have been set-up as a so-called synchronous recycling interferometer [75]. They provided an upper limit, h c < 1.4×10 −12 , on the amplitude of a stochastic GW background at 100 MHz [76].
• • Planar GWs induce resonant spin precession of electrons [78,79]. The same resonance is caused by coherent oscillation of hypothetical axion dark matter [80]. Recently, searches for resonance fluorescence of magnons induced by axion dark matter have been performed and upper bounds on the axion-electron coupling constant have been obtained [81,82]. These bounds can be translated to bounds on the amplitude of stochastic GWs: h c 1.3 × 10 −13 at 14 GHz and h c 1.1 × 10 −12 at 8.2 GHz [78,79]. T SMASH max 2 × 10 10 GeV (predicted by (pre-)heating in SMASH [7,8]).
• As mentioned earlier, in an external magnetic field, GWs partially convert into EMWs [64][65][66][67][68], which can be processed with standard electromagnetic techniques and detected [83], for example, by single-photon counting devices at a variety of wavelengths, cf. Fig. 9. The authors of Ref. [84] used data from existing facilities that have been constructed and operated with the aim of detecting axions or axion-like particles by their partial conversion into photons in magnetic fields: the light-shining-through-walls (LSW) experiments ALPS [85,86] and OSQAR [87,88], and the helioscope CAST [89,90]. They excluded GWs in the frequency bands from (2.7 − 14) × 10 14 Hz and (5 − 12) × 10 18 Hz down to a characteristic amplitude of h c < 6 × 10 −26 and h c < 5 × 10 −28 , at 95% confidence level, respectively. Using suitable EMW detectors sensitive to h c around its peak value at ∼ 40 GHz one may exploit such axion experiments also for the search of the CGMB, as we will show in the next subsection.
In summary: all the current upper bounds on the characteristic amplitude of stochastic GWs from direct experimental searches are many orders of magnitude above the CGMB predictions.

Prospects of EM detection of the CGMB in the laboratory
In this subsection, we will discuss the prospects of magnetic GW-EMW conversion experiments to probe the CGMB. We will first concentrate on GW-EMW conversion in available static magnetic fields in vacuum with dedicated detectors appropriate for the tens of GHz range and then proceed to a proposal exploiting an additional VHF EM Gaussian beam in order to generate a conversion signal which is first order in h c .

Magnetic GW-EMW conversion in vacuum
In this subsection, we consider experiments exploiting the pure inverse Gertsenshtein effect [64], cf. Fig. 9. To this end, we assume that stochastic GWs of amplitude h c propagate through a transverse and constant magnetic B in an evacuated tube of length L and cross-section A for a time ∆t. Then the average power of the generated EMW, per logarithmic frequency interval, at the terminal position of the magnetic field (z = L in Fig. 9) is obtained as [65,67,83,84] f dP (4.1) The index "2" denotes here the fact that the generated EMW power is second order in h c . The associated expected average number of generated photons, per unit logarithmic frequency interval, is given by These expressions are valid as long as the GWs and the generated EMWs are in phase coherence throughout their propagation in the magnetic field region. Under the assumption that the external B-field is surrounded by a circular beam tube of diameter d, coherent EMW generation is guaranteed if (see Appendix D):  Table 3: Parameters of the magnetic field regions of ALPs IIc [92,93], MADMAX [94,95], BabyI-AXO and IAXO [96], used to estimate the minimum detectable GW amplitude through magnetic conversions of GWs (gravitons) to EMWs (photons) in vacuum: B is the magnetic field magnitude, L is the magnetic field length, d is the diameter of the magnetized tube, and BLA 1/2 , with A = n tubes πd 2 /4, is the figure of merit for GW detection by magnetic conversion into EMWs. Also shown are the effective lower frequency cut-off Eq. (4.3) and the projected CGMB sensitivities around f = 40 GHz, exploiting these magnetic field regions and the benchmark values in Eqs. (4.6) and (4.10). ALPS IIc is not sensitive to the CGMB, because the lower frequency cut-off of its magnetic field region is around 5 THz.
where c 11 = 1.8 and d is the diameter of the beam tube. This effective lower frequency cut-off arises from the fact that the evacuated beam tube acts as an EM waveguide, in which the phase velocity of the EMW is higher than the phase velocity of light in vacuum Around the peak frequency of the h c spectrum, f h CGMB c peak ∼ 30 − 40 GHz (see Table 2), one may either use heterodyne (HET) radio receivers or single photon detectors (SPDs) to search for an EM signal that was generated from magnetic conversion of the CGMB.
The sensitivity of the HET technique is limited by thermal noise in amplifiers and mixers (for an introduction, see Ref. [91]). In this context, it is useful to introduce an effective signal noise temperature T S equal to the power of the generated EMW in a frequency bin ∆f around the peak frequency, Exploiting linear amplifiers with system noise temperature T sys , the signal-to-noise ratio is determined then by [91] where ∆f is the pre-detection bandwidth of the receiver, ∆t is the measurement time, and K rec is a receiver-system dependent dimensionless constant of order one 5 From this, we obtain the sensitivity of a magnetic GW-EMW conversion experiment with a heterodyne radiowave receiver to the CGMB as h CGMB c HET sens 9.65 × 10 −21 S/N 2 (4.6) The figure of merit of the magnetized region for conversion of GWs into EMWs is BLA 1/2 , cf. Eq. (4.6). This is shared also by LSW experiments exploiting optical cavities at the generation and regeneration side of the experiment and helioscopes searching for the magnetic conversion SM ΔN =10 -3 Figure 10: Characteristic amplitude h c of the CGMB (2.29) as in Fig. 8, but focusing on the frequency range near the peak and showing the predicted reach for experiments detecting GW-EMW conversions using heterodyne or single photon detectors in a static magnetic field, or resorting to an additional high-power 40 GHz EM Gaussian beam. The dashed line shows an optimistic projection for the Gaussian beam technique accounting for future technological improvements, see main text. of axions into photons or vice versa. In Table 3 we show the parameters of the magnetic field region of the next generation of axion experiments: the LSW experiment ALPS IIc [92,93], the haloscope MADMAX [94,95], and the helioscopes BabyIAXO and IAXO [96]. Unfortunately, the prospects to probe the CGMB exploiting these magnetic conversion facilities appear to be rather slim. For example, collecting the signal from all eight magnetized tubes of IAXO with a heterodyne radio receiver in a one year CGMB-EMW conversion experiment, the projected sensitivity given in Table 3, [h CGMB c (f ≈ 40 GHz)] HET sens ≈ 1.10 × 10 −22 , is about ten orders of magnitude above the CGMB predictions with early time equilibration, corresponding to initial temperatures above the Planck mass and an approximate saturation of the bound of Eq. (3.6), cf. Fig. 10.
The prospects are slightly better if progress is made on single photon detection at photon energies around ω = 2πf = 1.65 × 10 −4 eV[f /40 GHz]. The signal-to-noise ratio is given in this case by Here, denotes the number of signal counts in a time interval ∆t and an energy interval ∆ω = 2π∆f (cf. Eq. (4.2)), with being the single photon detection efficiency and ∆N D counts Γ D ∆t (4.9) the number of dark counts, in terms of the dark count rate Γ D . The sensitivity of a magnetic GW-EMW conversion experiment with an SPD detection system is then (4.10) If the experimental benchmark values chosen in Eq. (4.10) can be reached 6 , the SPD sensitivity is about three orders of magnitude better than the HET sensitivity. However, it is fair to say that, from today's perspective, vacuum magnetic GW-EMW conversion experiments will fail to beat the dark radiation constraint on h CGMB c by more than six orders of magnitude, cf. Table 3 and Fig. 10.

Magnetic GW-EMW conversion in a VHF EM Gaussian beam
The signal for magnetic GW-EMW conversion in vacuum, such as the generated EM power (4.1) or the number of generated photons (4.2), is of second order in the tiny amplitude h c of the passing GWs. A number of modified schemes have been proposed which introduce in the magnetic conversion region certain powerful auxiliary EM fields oscillating at the frequency of the gravitational wave, such as plane EMWs [101] or EM Gaussian beams (GBs) [102], to obtain GW-induced EMWs which are first order in h c . For h c 1, their signal strength overwhelms the one from the second order EMWs induced by the inverse Gertsenshtein effect. However, this does not mean automatically that the sensitivity of these modified magnetic conversion experiments is much larger than the one of the experiments based on the inverse Gertsenshtein effect in vacuum, because the powerful auxiliary EMWs tend to increase the noise floor and consequently to decrease the signal-to-noise ratio.
The arguably most promising of these modified magnetic GW-EMW conversion detection proposals exploits a VHF EM GB to induce a first order signal in magnetic GW-EMW conversion [102][103][104][105][106][107][108][109][110][111][112][113][114][115]. A continuous traveling wave EM GB with frequency f 0 , propagating in the zdirection with linear polarization along the x-direction, passes through a transverse static magnetic field, cf. Fig. 11. If a GW of frequency f = f 0 propagates along the z-direction, the resonant interaction of the GW with the EM fields of the GB and the static magnetic field will not only generate a longitudinal first order photon flux (denoted by n (1) z in Fig. 11), which will be swamped by the background EM flux n (0) z from the GB, but also a transverse first order photon flux (denoted by n (1) x in Fig. 11) in the direction perpendicular to the GB, which reads, for z ≥ L, In this context it is interesting to note that a quantum dot detector at 50 mK has achieved already a dark count rate of order mHZ in the photon energy range from 6.0 to 7.1 meV [97]. SPD with even lower dark count rates may be realized with Graphene-based Josephson junctions [98]. A research and development program on dedicated SPD at sub-THz frequencies is also motivated by future axion experiments, such as the LSW experiment STAX [99] and the haloscope TOORAD [100].
+ n x (1) x Ref Figure 11: Illustration of the proposal [102] to exploit an EM GB to induce a GW-EMW conversion signal which is first order in h c .

Δy Δz
Here E 0 is the amplitude of the electric field of the GB at the center of the beam at its waist, w 0 its waist radius, z R = πw 2 0 f 0 = 1.05 m[f 0 /(40 GHz)][w 0 /(0.05 m)] 2 its Rayleigh range, δ the relative phase between the GW and the GB, and [106] ψ (1) Depending on the overall sign of ψ (1) x , that is on δ, this flux points either in the positive or negative x-direction. The idea is then to place at x = ±x Ref reflectors (e.g. fractal membranes) parallel to the y-z plane, see Fig. 11, which could reflect and focus a portion of this flux to receivers and detectors placed at positions x = ±x Det with x Det > x Ref which are further away from the GB and therefore expected to suffer less from noise [103]. The number of signal photons within a bandwidth ∆f 0 around f 0 passing in a time interval ∆t through a detector surface element ∆S = ∆y∆z in the y-z plane, which extents from y = 0 to y = ∆y in the y-direction and from z = L to z = L + ∆z in the z-direction 7 , is then where 0 < η < 1 is the reflectivity of the reflector and F (1) x (x) = 1 4π (4.14) Numerical results for F (1) x (x, w 0 ), for particular values of f 0 , ∆y, ∆z, and L, are displayed in Fig. 12 (top panels). Right to the red lines, the Gaussian beam amplitude has dropped by a factor more than 1/e. We assume that placing the reflector in this region will cause only minor disturbances of the signal photon flux. Based on this assumption we find a benchmark value of 10 −5 for F 1 x which we have taken in Eq. (4.13). As benchmarks for the amplitude E 0 and the relative bandwidth ∆f 0 /f 0 of the GB we have taken in Eq. (4.13) values which can be achieved with a state-of-the-art free-running high-power (MW scale 8 ) gyrotron in this frequency range 9 [116,117]. 8 The total power of a GB is given by P0 = (π/4)w 2 0 E 2 0 = 1.30 MW [w0/(0.05 m)] 2 [E0/(5 × 10 5 V/m)] 2 . 9 It is interesting to note that a similar gyrotron has been proposed as the photon beam source of the axion LSW experiment STAX [99]. Therefore, in principle, one could extend STAX to a multi-purpose facility to search not only for axions, but also for GWs.
The zeroth order flux in the x-direction propagates radially out from the GB's axis, where α p 1 is the ratio of the y to x components of the GB electric field and [106] ψ (0) The corresponding number of background photons passing in a time interval ∆t through the detector surface is then Numerical results for F (0) x (x, w 0 ), for particular values of f 0 , ∆y, ∆z, and L, are displayed in Fig. 12 (bottom panels). As its benchmark we have taken in (4.17) a value of 10 −40 which is appropriate when the detectors are put at x Det = 0.45 m or 1 m, for L = 1 m or 5 m, respectively, cf. Fig. 12 (bottom panels). Therefore, this direct background from the GB is expected to be quite small if the receiver is placed sufficiently away from the beam. Moreover, it should occur simultaneously at the two detectors located at x = ±x Det , while the signal (n (1) x propagating towards x = 0, see Fig. 11), for fixed phase difference δ, occurs only at one of the two detectors. Nevertheless, this consideration neglects the possibility that the radiation from the GB is perturbed by the presence of the reflectors which in turn could disturb the signal photon flux. Furthermore the reflectors can be a noise source if the GB interacts with them. This kind of noise can be minimized by placing the reflectors at least right to the red lines in Fig. 12 where the GB amplitude has fallen off by 1/e [108]. However, the exact noise level, which is introduced by the interaction of the GB with the reflectors, has to be evaluated in a future study. If all these sources can be dealt with and the apparatus can be designed in such a way that finally the dark count rate (4.8) in SPD is the dominating background, then the sensitivity is given by Here we have taken current state-of-the-art benchmarks for the various experimental parameters. This is still three orders of magnitude above the CGMB predictions with early time equilibration, corresponding to T max > M P and an approximate saturation of the dark radiation bound (3.6), cf. Fig. 10. However, the latter may be reached and eventually surpassed by progress in the development of gyrotrons, SPD, and magnets. In fact, one may gain an order of magnitude in h c sensitivity by increasing the total power of the gyrotron by two orders of magnitude to ∼ 100 MW (and thus E 0 by one order of magnitude) and another order of magnitude in h c sensitivity by increasing the stable running time of the gyrotron by two orders of magnitude to ∆t ∼ 10 6 s. The remaining order of magnitude one may gain by developing SPD with a dark count rate of order 10 −5 s. Further improvements can be obtained by increasing the reflector size and by developing magnets with higher magnetic field and length. A sensitivity corresponding to T max = M P seems to be reachable in the not-so-distant future. In case of the detection of a signal, one may explore the frequency region around the nominal frequency of the gyrotron in a range ∆f 0 ∼ 10 −3 f 0 = 40 MHz [f 0 /(40 GHz)] by changing the acceleration voltage of the gyrotron.

Summary and outlook
Based on the pioneering work of Refs. [3,4], we have provided general formulae for the production rate of GWs from a primordial thermal plasma with sub-Planckian temperatures in an arbitrary theory with gauge fields, real scalars and Weyl fermions (cf. Sect. 2.1 and Appendix A) and derived general expressions for Ω CGMB (f ), the current energy fraction of those GWs per logarithmic frequency interval (see Sect. 2.2). It is found to peak around f Ω CGMB peak 80 GHz [106.75/g * s (T max )] 1/3 (cf. Eqs. (2.38), (2.43) and Table 2) -hence we chose CGMB (for Cosmic Gravitational Microwave Background) as the acronym for this stochastic background. Its overall magnitude scales approximately linearly with the maximum temperature which the primordial plasma attained at the beginning of the standard hot big bang era (cf. Eq. (2.21)) while also depending on g * s (T max ). For weakly coupled theories, the peak emission satisfies the approximate scaling of Eq. (2.44), (see also Eqs. (2.39) and (2.40)) implying that for a given T max the SM will typically maximize the CGMB with respect to its value in weakly coupled extensions. With the leading behaviour of the peak frequency and the magnitude of the CGMB being determined by T max and g * s (T max ), the CGMB can therefore act as a hot big bang thermometer and, additionally, allow a measurement of the number of thermalized BSM degrees of freedom, g * s (T max ) − g * s,SM (T max ). As special cases, we have determined the CGMB spectrum for the cases of a SM (cf. Sect. 2.3), a νMSM, a SMASH (cf. Sect. 2.4), as well as an MSSM (cf. Sect. 2.5) plasma. We confirmed that the leading model dependence is indeed captured by the effects of g * s (T max ), so that within a broad class of weakly coupled SM extensions, a simple comparison of a hypothetical measurement of the CGMB peak with the SM prediction would allow to estimate g 1/3 * s (T max ) and T max in a model-independent manner and with respective theoretical accuracies that should be better than the MSSM results of 15% and 40% (e.g. ∼ 1% and 5% in SMASH).
The previous features of the CGMB apply for T max < M P , while for larger early-time temperatures one expects gravitons to thermalize and lead to the blackbody spectrum of Eq. (2.45), with peak frequencies and maxima scaling with g * s (M P ) as in Eqs. (2.24) and (2.46), and independent of the concrete value of T max and of any additional model details. Here, a possible detection would allow a precise determination of g * s (T max ).
We have found that current dark radiation constraints from BBN and CMB on the total energy density fraction in GWs cannot yet probe the out-of-equilibrium gravitational wave emission with T max < M P , as a naive application of the constraints implies a trans-Planckian upper bound on T max around 10 19 GeV, cf. Table 1. Nevertheless we find the intriguing result that the CGMB background with early time equilibration (i.e. with T max > M P ) just about saturates the dark radiation constraint of Eq. (3.6), to be compared with the prediction of Eq. (3.7). The former bound corresponds to homogeneous initial conditions for the gravitational waves, as appropriate under the assumption of thermal equilibrium at early times. Applying the bound of Eq. (3.6) and using the determination of g * s,SM of Appendix B (illustrated in Fig. 3) would discard a CGMB with early time equilibration in models in which g * s is augmented by less than ∼ 4 units.
Further improvements on the dark radiation constraints would discard early time equilibration of the gravitational waves in a large classes of models, and start constraining the out-of-equilibrium CGMB for sub-Planckian T max . In case that future CMB constraints on dark radiation reach the theoretical uncertainty from the pure SM expectation (∆N ν = 0.001), the upper bound of T max can be improved to sub-Planckian values around 2 × 10 17 GeV, cf. Table 1.
Further progress should come from direct detection of GWs. However, all the current upper bounds from direct searches for stochastic GWs and also the projected sensitivities of planned GW detectors are at least nine orders of magnitude away from the prediction of the maximally allowed characteristic amplitude h c of the CGMB respecting the dark radiation constraint (corresponding to the CGMB with early time equilibration), cf. Fig. 8.
Conversion of GWs into EMWs in a static magnetic field has been identified as a promising search technique for stochastic GWs at frequencies around f  Table 2. We investigated the prospects of GW-EMW conversion in state-of-the-art superconducting magnets used in present and near future axion experiments and the detection of the generated EMWs/photons with dedicated detectors appropriate around 40 GHz. The projected sensitivity of this technique, cf. Table 3 and Fig. 10, turned out to fail to beat the dark radiation constraint on h CGMB c by about six orders of magnitude. We then investigated the prospects of a proposal exploiting an additional EM Gaussian beam, delivered by a MW-scale 40 GHz gyrotron, propagating along the magnetic conversion region in order to generate a transverse EMW conversion signal which is first order in h c . Assuming state-of-the-art benchmarks for the gyrotron, the detector performance and the magnetic field strength and length, the projected sensitivity in h c is still three orders of magnitude above the the maximum amplitude of the equilibrated CGMB (see Eq. (2.46)) which saturates the dark radiation bound, cf. Fig. 10. However, the latter may be reached by progress in the development of gyrotrons towards higher power and stable run time, single photon detection towards lower dark count rate, and superconducting magnets towards higher magnetic fields. The direct detection of the CGMB at a level corresponding to T max ∼ M P by such a magnetic conversion experiment seems possible, although challenging. In this connection, it should be emphasized that the search for the CGMB is truly a critical endeavour. Any measurement of T max above 6.6 × 10 15 GeV would be ground-breaking, since it would rule out standard inflation as a viable pre hot big bang scenario.
It would be very interesting to investigate the CGMB also in other BSM models with a complete and consistent cosmological history and thus giving a prediction of T max , such as for example the model in Ref. [118]. Furthermore, it seems worthwhile to explore more deeply the possible synergies between axion and GW experiments which we have been touching upon in this paper.
Universe" -390833306. CT acknowledges financial support by the DFG through SFB 1258 and the ORIGINS cluster of excellence.
A Loop functions for the rate of GW production from the primordial plasma at full leading order Here we collect formula for the loop functions η gg (k), η sg (k), η f g (k), η sf (k) appearing in Eq. (2.5).
We use the results of Ref.
[4] with a simplified notation. One can define a set of six integrals I ± (k), J ± (k), K, L, in terms of which the above functions are expressed as: The loop integrals are given next: In the above expressions, one has while the functions L ± i (x, y), M ± i (x, y), with i = 1, 2, 3 are defined below: , M ± 2 (x, y) = Li 2 (±e − 1 2 (x+y) ) + Li 2 (±e − 1 2 (x−y) ), As explained in the main text, the integrals above include a subtraction of infrared divergences, so that they remain finite. The subtraction corresponds to the last term in the integrand of I ± (k).

B Effective number of degrees of freedom of the SM, the νMSM and the SMASH plasma
In the primordial plasma, the thermal contributions to the effective potential correspond to the free-energy density. Thermodynamic relations imply that the latter is equal to minus the plasma's pressure p: In the equation above, φ i denote the scalar field backgrounds, and the subtraction of the second term above guarantees that the pressure is zero in the vacuum (T = 0). We will assume that the system relaxes to the minimum of the effective potential, so that we will evaluate the backgrounds at the configurations φ i =φ i that extremize V eff : The effective potential with its finite-temperature correction at one-loop plus higher-order QCD effects can be written as: Above, V is the tree-level potential, V CW the usual one-loop Coleman-Weinberg contribution at zero temperature, V T the one-loop thermal correction to the potential, and V QCD includes two and three-loop thermal corrections induced by QCD effects, which are known for zero field backgrounds. Note that we have included a temperature-dependence in the Coleman-Weinberg vacuum piece. The reason is that achieving accuracy in the finite-temperature effects near a phase transition requires a resummation of thermal contributions to the 2 point functions -daisy resummation -which can be implemented by replacing the effective masses in the vacuum propagators by temperature corrected masses; this explains the T dependence in the Coleman-Weinberg piece. The daisy resummation is usually performed by keeping the leading terms of the corrections to the two point functions in a high-temperature expansion. However, as noted in Ref. [30], achieving accuracy across a phase transition, i.e. for temperatures below the masses gained by particles during the transition, requires to go beyond the high-T expansion. As we consider temperatures above the SM's electroweak crossover, for the SM and νMSM we use a traditional daisy resummation. For SMASH, for which the PQ phase transition lies much above the electroweak scale, we implement the improved resummation of Ref. [30] for the corrections involving fields that acquire masses during the PQ transition.
In the Landau gauge and in the the MS scheme, V CW (φ i , T ) and V T (φ i , T ) are given by

(B.4)
In the previous formulae, V, S, F and G correspond to massive gauge bosons, real scalars, Weyl fermions and ghosts. It is assumed that for gauge bosons one has to sum over the three polarizations that propagate in the Landau gauge, while for fermions one has to sum over two helicities, and for ghosts one should count one degree of freedom per generator of each gauge group. m 2 V /S/F (φ i , T ) are the field-dependent masses in the background of φ i , including the thermal corrections from the (improved) daisy resummation. µ is the renormalization scale, while the thermal loop functions J B and J F are The last missing piece in the effective potential in Eq. (B.3) corresponds to higher loop QCD contributions, which have been computed in Ref. [22] up to three-loop order in a theory with arbitrary massless QCD flavours. This allows to apply the QCD corrections for all the models of interest.
Once the pressure is calculated from the minimization of the effective potential that follows from the above equations, the energy density ρ, standard thermodynamic relations allow to compute the entropy density s and the specific heat capacity c = 1/V (∂U/∂T )| V as From the former quantities one defines the effective numbers of degrees of freedom g * ρ , g * s , g * c by writing In our calculations of the effective potential, we use a renormalization scale proportional to the temperature, µ = κ2πT with κ ∈ {1/2, 1, 2}. It should be noted that for the SM our perturbative calculations are not trustworthy for temperatures below the electroweak phase transition, as the daisy resummation does not capture the decoupling of the degrees of freedom acquiring masses during the transition. Furthermore, for lower temperatures the QCD interactions become nonperturbative, and other techniques are necessary [23]. In this paper we neglect the gravitational wave production at temperatures below the electroweak phase transition.
For SMASH we introduce an additional correction coming from the loss of chemical equilibrium of the axion, which implies a separate conservation of the entropy of the axion and the rest of the thermal bath, and thus separate temperatures for both, to be denoted as T axion and T . Labeling the results for g SMASH * ρ , g SMASH * s , g SMASH * c under the assumption of chemical equilibrium for the axion with the superfix "eq", the quantities accounting for decoupling can be written as [30] g SMASH * ρ = g SMASH,eq * ρ

C Upper bound on T max in inflationary cosmology
The hypothesis of inflation postulates a period of accelerated expansion,ä > 0, in the very early universe [119,120], preceding the standard radiation-dominated era. It offers a physical model for the origin of the initial conditions of the HBBC. In slow-roll inflationary cosmology [121,122], the energy scale at the end of inflation, ρ I , can be inferred from the amplitude A s of scalar perturbations generated during inflation and the ratio r ≡ A t /A s of the amplitudes of tensor, A t , and scalar perturbations via where H I is the Hubble expansion rate during inflation. The upper limit on r obtained from CMB observations of the BICEP2/Keck Array and Planck Collaborations provides an upper bound on ρ I [45], ρ I < 1.6 × 10 16 GeV 4 (95% CL) .

(C.2)
This may be turned into an upper bound on the maximum temperature of the post-inflationary universe by assuming instantaneous and thus maximally efficient reheating to a radiation dominated universe with energy density In SMASH [8], g * ρ (10 15 GeV) 124.5.

D 3D effects in GW-EMW conversion experiments D.1 Waveguide effect
If the magnetic conversion volume is surrounded by a waveguide the phase velocity inside the waveguide will be different from the speed of light. We illustrate this for a circular waveguide which is infinite in z-direction and has radius a. For TE modes we have E z = 0 everywhere and the normal derivative of B z vanishes at the boundaries. After the application of all boundary conditions we obtain for the longitudinal B-field: where F, φ 0 are constants, ν is an integer, r, φ and z are the three coordinates in the cylindrical coordinate system, J ν is the Bessel function of the first kind and k c = pνn a is the transversal momentum, where p νn is the nth zero of J ν (x). The total wave vector is given by k 2 = k 2 c +k 2 z = ω 2 . From the longitudinal field B z we can derive all other electromagnetic fields in the cylindrical waveguide. The phase velocity is: If the radius of the waveguide is much larger than the wavelength we can expand the square root and get: v p ≈ 1 + 1 2 p νn λ 2πa 2 .

(D.3)
Therefore the phase velocity is larger than the speed of light. For a coherent and thus efficient conversion of gravitons into photons inside a waveguide, they should be in phase all along the waveguide. Therefore the phase difference between gravitons and photons has to be smaller than π at the end of the waveguide: where k G = ω is the graviton dispersion, k z is the photon wave vector in z-direction and L the length of the external B-field. For the dominating T E 11 mode we find p 11 = 1.8 and therefore in this case the prefactor is 2π pνn 2 ≈ 12.

D.2 Diffraction in open systems
If no waveguide encloses the magnetized region we only have to take into account diffraction effects. For a magnetized region of length L along the z-direction with a circular shape in the xy-plane we obtain the diffraction angel θ in the far field as: where d is the diameter of the magnetized region. For IAXO we get θ = 0.6 • and therefore diffraction is negligible. In the case of ALPS IIc the diffraction angle is much larger and due to the huge length of ALPS IIc it cannot be neglected.