Gravitational wave background from Standard Model physics: Qualitative features

Because of physical processes ranging from microscopic particle collisions to macroscopic hydrodynamic fluctuations, any plasma in thermal equilibrium emits gravitational waves. For the largest wavelengths the emission rate is proportional to the shear viscosity of the plasma. In the Standard Model at T>160 GeV, the shear viscosity is dominated by the most weakly interacting particles, right-handed leptons, and is relatively large. We estimate the order of magnitude of the corresponding spectrum of gravitational waves. Even though at small frequencies (corresponding to the sub-Hz range relevant for planned observatories such as eLISA) this background is tiny compared with that from non-equilibrium sources, the total energy carried by the high-frequency part of the spectrum is non-negligible if the production continues for a long time. We suggest that this may constrain (weakly) the highest temperature of the radiation epoch. Observing the high-frequency part directly sets a very ambitious goal for future generations of GHz-range detectors.


Introduction
Gravitational waves offer a possible way to observe phenomena taking place very early in the history of the universe. Famously, long-wavelength waves are produced during a period of inflation (cf. e.g. refs. [1,2]). However at higher frequencies gravitational waves can also probe post-inflationary non-equilibrium phenomena, such as preheating [3]- [9], topological defects [10]- [18], bubble dynamics related to a first-order phase transition [19]- [24], or noisy turbulent motion [25]- [29]. Recent numerical simulations start to account for both bubble dynamics and the subsequent motion [30,31]. For a review concerning post-inflationary sources and the associated observational prospects, see ref. [32].
It is well-known that gravitational waves are being produced in thermal equilibrium as well (for an example, see ref. [33]). In thermal equilibrium particles scatter on each other, which implies the presences of forces and accelerations. However, for a physical momentum k > T , the thermal production rate is suppressed by e −k/T , because the energy carried away by the graviton must be extracted from thermal fluctuations. Since typical particle momenta are ∼ 3T and scatterings are proportional to the coupling strengths responsible for the interactions, it may be assumed that the rate is suppressed by ∼ αT 3 e −k/T /m 2 Pl , where α is a fine-structure constant and m Pl is the Planck mass. In a weakly-coupled system such as the Standard Model, this rate is small.
On the other hand, a thermal system experiences also long-wavelength fluctuations not associated with single particle states. At the smallest momenta k ≪ T these can be called hydrodynamic fluctuations [34]. We are not aware of an estimate of a corresponding contribution from equilibrium Standard Model physics to the gravitational wave background, and one purpose of the present note is to provide for one. (Very similar physics, although boosted by a conjectured turbulent cascade, has recently been discussed in ref. [29].) In contrast to the non-equilibrium phenomena mentioned above, the "advantage" of a thermal contribution is that it is guaranteed to be present. Another purpose of our study is to roughly estimate the production rate at k > T , and to motivate the need for a complete computation.

JCAP07(2015)022
In order to be more specific, consider the closely analogous case of the production rate of photons from a plasma which is neutral but has electrically charged constituents. Even though the expectation values of the electromagnetic charge density and current vanish, n em = 0, J em = 0, thermal fluctuations do induce charge fluctuations which have a non-zero root mean square value: where V denotes the volume, χ em a susceptibility, and . . . a thermal expectation value. The susceptibility is non-zero even without interactions, for instance for a plasma of free massless Dirac fermions representing electrons and positrons it reads χ em = e 2 T 2 /3, e 2 ≡ 4πα em . Because of diffusion, the charge fluctuations induce electromagnetic currents, and currents in turn source photons. Currents can also directly originate from fluctuations. Assuming that the photons produced do not equilibrate as fast as the plasma, which is the case for instance for the plasma generated in heavy ion collision experiments, the thermal average of their production rate can be evaluated. A text-book computation shows that the rate per unit volume can be expressed as [35,36] dΓ γ (k) µ,k denote polarization vectors. For k = k e 3 , the polarization sum only couples to the transverse components J 1,2 em . For small k ≪ T , operator ordering plays no role in eq. (1.2), and the fluctuations are also uncorrelated in space and time [34]. Their amplitude is related to diffusion or, equivalently, to conductivity (σ). This yields finally where we inserted the parametric form of the conductivity of a QCD plasma [37,38]. For large k 3T , in contrast, the rate originates from particle scatterings rather than hydrodynamic fluctuations, and has the parametric form [39,40] dΓ γ (k) In the following we show that results analogous to eqs. (1.3) and (1.4) apply to gravitational waves, just with the replacements α em → T 2 /m 2 Pl and α s → α. Our presentation is organized as follows. After deriving an expression for the gravitational wave production rate in section 2, we analyze the structure of the energy-momentum tensor correlator for k ≪ T in section 3. The quantity parametrizing this structure, the shear viscosity, is briefly discussed in section 4. In section 5 we turn to the other case k 3T and compute the logarithmically enhanced terms in this regime. The results are embedded in a cosmological background in section 6 and compared with a well-studied non-equilibrium source in section 7. Section 8 offers some conclusions and an outlook. 2 Production rate of gravitational waves from thermal equilibrium As a first ingredient, we consider the rate at which energy density is emitted in gravitational waves. The derivation can be carried out in two different ways: by treating gravitons as quantized particles, or through a purely classical analysis. We start with the first method, leading to a result analogous to eq. (1.2). We work first in Minkowskian spacetime, adding cosmological expansion in section 6.
The linearized equation of motion for the metric perturbation h ij in the traceless transverse gauge readsḧ where G = 1/m 2 Pl . The right-hand side of this equation plays the role of the electromagnetic current in the photon case. The classical energy associated with gravitational waves reads where V is a volume. It is well-known that the corresponding energy density cannot be localized. However, if we express a free h TT ij as a usual linear combination of forward and backward-propagating plane waves, and omit fast oscillations exp(±2iωt), then eq. (2.2) can be re-interpreted as a Hamiltonian with a familiar canonical form: Here . . . denotes an average over an oscillation period. From eq. (2.3) canonically normalized fields can be identified asĥ TT ij ≡ h TT ij / √ 32πG. According to eq. (2.1) they are sourced as ∂ 2 We can now directly overtake eq. (1.2) for the production rate of gravitons, by replacing the polarization vectors accordingly. Subsequently, weighting the production rate by the energy carried by individual quanta, we obtain The sum over the polarization vectors yields where P ij ≡ δ ij − k i k j /k 2 . Choosing k = k e 3 and rotating subsequently the diagonal In order to be convinced that eq. (2.4) is correct, let us repeat the analysis on a purely classical level. Fourier transforming eq. (2.1) with respect to spatial coordinates and denoting the retarded Green's function related to the time evolution by ∆(t, k), its time derivative readṡ . Dividing by volume, the energy density corresponding to eq. (2.2) can then be expressed as 3 . Following now a standard argument, let us assume that the sources switch off before the observation time t. Then the upper bounds of the time integrals can be treated as independent of t, and we can average over fast oscillations within the integrand: Taking a time derivative and assuming that φ is a function of the time difference 1 yieldṡ Going finally back to configuration space and making use of translational invariance in spatial and temporal directions we get Given that eq. (2.5) defines a projection operator to the TT modes and that in the classical limit operator ordering plays no role, eq. (2.10) indeed agrees with eq. (2.4) for k ≪ T .

Correlation function in the tensor channel
Having obtained eq. (2.6), the next task is to determine the shape of the energy-momentum tensor correlator in momentum space. Here we do this for small light-like four-momenta (ω, k α 2 T ), returning to the regime k 3T in section 5. Consider hydrodynamic fluctuations associated with a local flow velocity v i around an equilibrium state at a temperature T . To first order in gradients and in v i , 2 the energymomentum tensor has the form where e, p, ζ, η are the energy density, pressure, bulk viscosity, and shear viscosity, respectively. The equation for energy-momentum conservation asserts that ∂ 0 T 0j + ∂ i T ij = 0 ∀j ∈ {1, 2, 3}. Let us consider a plane wave perturbation with a momentum vector k = k e 3 . Then the equations of motion for the transverse velocity components (v ⊥ · k = 0) decouple from the equations relating v 3 and ∂ 3 p. The resulting system is immediately integrated to obtain

JCAP07(2015)022
We now consider the 2-point correlator where the operator ordering is only relevant in the quantum theory. This correlator is symmetric in t → −t and has a classical limit. Therefore equations (3.1) and (3.3) lead to a hydrodynamic prediction for the transverse components where we returned to configuration space for the equal-time correlator. Let us take k to be very small, and look for the leading term in this limit. In the equal-time correlator we can send k → 0. Then it equals the susceptibility related to the total momentum in the i ′ -direction: Even though the average momentum is zero, its susceptibility is non-zero, in analogy with eq. (1.1): Despite including an integral over operator correlations at short distances, this exact equation is ultraviolet finite just like eq. (1.1) (for a rigorous discussion see ref. [41]). Now, a Ward identity related to energy-momentum conservation asserts that Therefore, eqs. (3.5)-(3.7) can be re-expressed as Taking lim ω→0 lim k→0 from here yields the well-known Kubo formula for η.
Of interest to us is not the correlator of eq. (3.9) (known as the "shear channel") but the corresponding correlator for the spatial components transverse to k (known as the "tensor channel"). It can be argued, however, that its functional form is closely related to that in eq. (3.9). The tensor components also experience hydrodynamical fluctuations; but, in our coordinate system with k = k e 3 , these are related to the velocity gradients ∂ 1 v 2 or ∂ 2 v 1 , cf. eq. (3.2). These components are decoupled from the equations of motion following from energy-momentum conservation. Therefore, they are not represented by smooth differentiable functions responsible for the transfer of hydrodynamic information from one point and time to another; rather, nearby points are uncorrelated, as is the case for generic thermal fluctuations [34]: Consequently a Fourier transform like in eq. (3.9) is independent of ω, k. Putting finally ω = k and sending k → 0 so that the distinction between spatial directions disappears, we can fix the coefficient Φ through a comparison with eq. (3.9): This is the main result that is needed below. We remark that, apart from the physical arguments discussed above, the same expression can be derived more formally by a linear response analysis related to a metric perturbation (cf. e.g. ref. [42]).

Estimate of shear viscosity at T > 160 GeV
Shear viscosity (η) is a macroscopic property of a plasma that originates from the microscopic collisions that its constituents are undergoing. It is inversely proportional to a scattering cross section and therefore large for a plasma in which there are some weakly interacting particles. In the Standard Model above the electroweak crossover, right-handed leptons are the most weakly interacting degrees of freedom, changing their momenta only through reactions mediated by hypercharge gauge fields.
Omitting for the moment all particle species which equilibrate faster than right-handed leptons, the shear viscosity can be extracted from refs. [37,38]: where m D1 = 11/6 g 1 T is the Debye mass related to the hypercharge gauge field. Inserting g 1 ∼ 0.36 for the gauge coupling we obtain We use this value for order-of-magnitude estimates below. If we increase the temperature above 160 GeV, the hypercharge coupling g 1 grows and the weak and strong couplings g 2 , g 3 decrease. Presumably, the top Yukawa coupling h t and the Higgs self-coupling λ are also of a similar magnitude. In this situation the analysis of refs. [37,38] should be generalized to include a scalar field and a more complicated set of reactions. Even though conceptually straightforward, implementing and solving numerically the corresponding set of rate equations is a formidable task and beyond the scope of the present investigation. We note, however, that the shear viscosity is likely to decrease with increasing g 1 , so that eq. (4.2) should represent the most "optimistic" estimate from the point of view of detecting a thermally emitted low-frequency gravitational wave background.
5 Leading-logarithmic production rate at large momentum Before turning to numerical estimates we wish to complete the qualitative picture concerning the thermal graviton production rate by considering the case of "hard momenta", k ∼ 3T . A full computation of the rate in this regime represents a complicated task, similar to the full computation of the shear viscosity when all couplings are of the same order of magnitude, and is postponed to future work. In contrast to the shear viscosity, for hard momenta the result is dominated by the largest couplings, in particular the strong gauge coupling. If we Figure 1. Processes leading to a logarithmically enhanced graviton production rate. Wiggly lines denote gauge bosons; arrowed lines fermions; dashed lines scalars; and a double line a graviton. By k, p ∼ 3T we denote typical momenta of the scattering particles, whereas the filled blob indicates that the vertical rung carries a soft spacelike momentum transfer (t ∼ −q 2 ⊥ ∼ −g 2 T 2 , where q ⊥ · k = 0) so that the gauge boson needs to be Hard Thermal Loop resummed. restrict to logarithmically enhanced terms (cf. eq. (1.4)) then it can be shown that only the gauge couplings (g 1 , g 2 , g 3 ) contribute at leading order.
An elegant way to determine the logarithmically enhanced terms has been discussed in ref. [43], section 4.2. Scatterings experienced by soft space-like gauge bosons, the vertical rung in figure 1, correspond to Landau damping, and can be represented within the Hard Thermal Loop (HTL) effective theory [44,45]. Computing the 2-point correlator of T 12 within the HTL theory and noting that only one of the gauge bosons attaching to the graviton vertex can be soft at a time, yields (for q ≪ k ∼ 3T ) where f B is the Bose distribution; q ≡ q · k/k; m 2 D is a Debye mass squared; Λ indicates that this treatment only applies to soft modes q ⊥ ≪ 3T ; and ρ T/E are spectral functions corresponding to the "transverse" and "electric" polarizations, respectively. In the last step we made use of a sum rule for the HTL-resummed gluon propagator that has been derived in refs. [46,47].
The integral in eq. (5.1) happens to be identical to that appearing in the context of the jet quenching parameterq in QCD [47]. Carrying it out and inserting the result into eq. (2.6), we obtain where d 1 ≡ 1, d 2 ≡ 3, d 3 ≡ 8; m Di is the Debye mass corresponding to the gauge group U(1), SU(2) or SU(3), respectively; g 2 ∈ {g 2 1 , g 2 2 , g 2 3 , h 2 t }; and the ultraviolet scale within the logarithm has been (arbitrarily) taken over from eq. (4.1). The Debye masses read m 2 D1 = 11g 2 1 T 2 /6, m 2 D2 = 11g 2 2 T 2 /6, and m 2 D3 = 2g 2 3 T 2 . Because of the largest multiplicity, the result is dominated by the QCD contribution. We note that a similar computation for t-channel fermion or Higgs exchange does not lead to logarithmic enhancement. This applies in a local Minkowskian frame. The function φ, is quantitatively correct at k α 2 T whereas at k 3T it only represents the qualitative structure (in particular the coefficient "5" inside the logarithm is but a convention, and there could be substantial non-logarithmic contributions from h 2 t or from O(g)-suppressed effects like in the case of the jet quenching parameterq [47]). We would now like to re-express the result in an expanding cosmological background, and subsequently obtain numerical estimates. As our reference temperature we take that corresponding to the electroweak crossover in the Standard Model, T 0 ≡ 160 GeV (cf. e.g. refs. [48,49]).
The basic equations needed from cosmology are (for a flat spatial geometry) -8 -

JCAP07(2015)022
where H is the Hubble rate, a is the scale factor, and s(T ) is the entropy density. Combining the two equations in eq. (6.3), the relation between time and temperature can be expressed as where c s is the speed of sound, c 2 s (T ) = p ′ (T )/e ′ (T ). The energy density carried by gravitational waves is of the form ρ GW (t) = k k f (t, k), where f is a phase space distribution. Making use of the known evolution equation for f in an expanding background, the energy density can be seen to evolve as where R(T, k) = 32πηT φ(k/T )/m 2 Pl in the notation of eq. (6.1). Given that (∂ t + 3H)s = 0, the factor 4H can be taken care of by normalizing ρ GW by s 4/3 . Subsequently the equation can be integrated, by assuming that at an initial time t min (corresponding to a maximal temperature T max ) there were no (thermally produced) gravitational waves present: .

(6.6)
Taking into account that momenta redshift as k(t) = k 0 a(t 0 )/a(t) and expressing the momentum space integrals in terms of k 0 finally yields where we also inserted H from eq. (6.3). If we approximate c 2 s ≈ 1/3; assume all thermodynamic functions to scale with their dimension (s =ŝT 3 , η =ηT 3 , e =êT 4 , withŝ,η,ê roughly constants at T > T 0 ); and consider T max ≫ T 0 , then This result is plotted in figure 2, after a redshift of k 0 /T 0 to a current-day frequency. For a given mode k 0 , production starts at a maximal temperature when the argument of φ is of order unity; since the entropy density roughly scales with T 3 for T > T 0 , this poses no particular constraint if we restrict to k 0 T 0 . The horizon radius of a given period redshifts, and we are interested in causal physics taking place within the horizon. For instance, a temperature leading to a horizon radius comparable to the planned eLISA [50] arm length, ∼ 10 6 km, is T max ∼ 10 6 GeV. Let us take this as an example. Insertingη ∼ 400,ê ∼ 35, T max ∼ 10 6 GeV into eq. (6.8), we thus get

JCAP07(2015)022
The Hubble radius (H −1 ) of the electroweak epoch (T = T 0 ) corresponds to ∼ 10 10 km today, so if we consider wavelengths extending up to ∼ 10 6 km, then k 0 ∼ αH with α ∼ 2π × 10 4 . In this situation k 0 /T 0 can be estimated as Inserting into eq. (6.9) we find a very small energy fraction. However the fraction is larger if we consider the total amount of energy in gravitational waves, which originates dominantly from k 0 ∼ T 0 ; this energy is constrained to be below that corresponding to one equilibrated relativistic degree of freedom [51,52], ln k 0 Ω GW ≪ 1/100. We return to this consideration around eq. (8.2), but first compare the infrared part with non-equilibrium processes.
7 Order-of-magnitude comparison with a non-equilibrium source In order to get a qualified impression about the magnitude of the thermal background, consider the well-studied case of a first-order phase transition at the electroweak epoch. As before Ω GW denotes the ratio of the energy densities of gravitational waves and radiation at T 0 ∼ 160 GeV. Today, the energy density in radiation corresponds to Ω rad ∼ 5 × 10 −5 , and the projected eLISA sensitivity is Ω eLISA ∼ 10 −11 . So, in order to be detectable, we hope to find a signal in the range Ω GW ∼ 2 × 10 −7 at the electroweak epoch. Apart from the overall magnitude, an important feature of any observatory is that its sensitivity peaks in a certain frequency range. For eLISA, this is f ∼ (10 −3 . . . 10 −2 ) Hz. This corresponds to a distance scale ℓ B ∼ (10 −3 . . . 10 −2 ) ℓ H in terms of the horizon radius of the electroweak epoch, where we have defined ℓ H ≡ H −1 . Now, the overall signal from a first-order transition is traditionally argued to be of the form [19]- [24] Ω GW ∼ v n w κ 2 L e where v w is a bubble wall velocity; n > 0; L is the latent heat released in the transition; κ < 1 parametrizes the efficiency at which energy is converted into gravitational waves [22]; and ℓ B is the typical bubble separation. The spectrum peaks at momenta corresponding to the bubble separation scale, k max ∼ 2π/ℓ B . For k < k max , a behaviour ∼ k 3 has been found, whereas at k > k max the spectrum falls off, perhaps as 1/k or 1/k 2 [24]. A recent numerical study shows that the complicated dynamics following bubble collisions continues for a long time and thereby boosts eq. (7.1) by a factor ∼ ℓ H /ℓ B , with the price of a more rapid decay of the power spectrum at large k [30]. If we optimistically insert ℓ B ∼ 0.01ℓ H into eq. (7.1), it still remains a challenge to get Ω GW ∼ 10 −7 out. Basically, a very large latent heat L/e ≫ 0.01 would be needed. This requires a drastic modification of the Higgs sector, which normally carries just a handful of degrees of freedom in comparison with ∼ 100 contributing to e. However, with somewhat less drastic assumptions, and including a boost by ∼ ℓ H /ℓ B from ref. [30], numbers like Ω GW ∼ 10 −10 could be obtained, which is also of interest for future generations of observatories.
In any case, inserting α ∼ 2πℓ H /ℓ B ∼ 10 3−4 and T max ∼ 10 6 GeV into eqs. (6.9), (6.10), the thermal background would be ∼ 40 orders of magnitude below the desired level at the peak eLISA frequency (cf. also figure 2). What is significant about the thermal background, though, is that it continues to grow with k for another more than 10 decades, and therefore -10 -

JCAP07(2015)022
eventually overtakes the decaying non-equilibrium signal at short distance scales. In fact, originating as it does from fluctuations at the scale k ∼ 3T max and red-shifting as dictated by entropy conservation, the peak power is in the range k ∼ 3T dec (3.9/106.75) 1/3 ∼ T dec at the time of photon decoupling, and in the microwave range today. Therefore it falls in the range of recently conceived high-frequency experiments [53]- [57].

Conclusions and outlook
We have estimated the magnitude and shape of the gravitational wave background that is produced by Standard Model physics during the thermal history of the universe until the temperature T 0 ≈ 160 GeV corresponding to the electroweak crossover, cf. eqs. (6.2), (6.8) and figure 2. The infrared part could have been of potential interest in that the forthcoming eLISA experiment is probing sub-Hz frequencies with unprecedented precision. Unfortunately, we have found that in this range the thermally produced gravitational wave signal is many orders of magnitude below the observable level, cf. section 7.
In general, the thermally produced gravitational wave background resembles a bit the blackbody spectrum of photons and neutrinos. Its shape is not the same because gravitational waves never equilibrate. Therefore the shape has to be determined by a dynamical computation which has not been carried out even at full leading order. Nevertheless, it is already clear from a leading-logarithmic estimate that the peak of the power today lies in the same microwave domain as for photons and neutrinos (cf. figure 2). Therefore, the best observational prospect lies with high-frequency (0.1 -4.5 GHz) experiments [53]- [57]. At the current stage it seems challenging to reach a sensitivity below Ω GW ∼ 10 −5 [56], whereas an optimistic theoretical expectation would be Ω GW ∼ Ω rad × 1 100 where 1/100 corresponds to the maximally allowed fraction in gravitational waves at the electroweak epoch when all degrees of freedom are relativistic, and 100 GHz to the frequency associated with generic blackbody radiation. So, there is surely a long way to go till detection. There is, however, one consideration which can already be carried out. Indeed, unlike neutrinos, gravitational waves must not carry as much energy density as one relativistic degree of freedom [51,52]. This constrains the total energy density stored in them and, given that the production rate peaks at the maximal temperature, the maximal temperature reached. The total energy density corresponding to eq. (6.8) can be estimated as where we varied φ between two limits: the factor 8 originates if we adopt the form of φ appearing on the second line of eq. (6.2), setting the unknown constant to zero and the couplings to values mentioned in the caption of figure 2, whereas the factorη/3 originates if we use the first line of eq. (6.2) and cut off the integral at k 0 = T 0 : T 0 0 dk 0 k 2 0 = T 3 0 /3. According to Planck data [58] only a small fraction of a relativistic degree of freedom beyond those in the Standard Model can be permitted, so at T 0 ∼ 160 GeV we must require Insertingê ∼ 35,η ∼ 400 we obtain T max 10 17 . . . 10 18 GeV. This is not a very strong constraint, 3 but the estimate could be sharpened with more knowledge about the function φ.
To summarize, a determination of the function φ, defined by eq. (6.1), beyond the leading-logarithmic terms that we have obtained here, seems to pose an interesting problem. This computation represents a well-defined challenge in thermal field theory, analogous to that for the photon production rate from a QCD plasma [39,40] or the right-handed neutrino production rate from a Standard Model plasma [43,59]. It is technically more challenging, because every single particle species carries energy and momentum, and therefore we leave the practical implementation to future works. (A somewhat related computation, but for the off-shell kinematics k = 0, ω T , has been presented in ref. [60].)