Thermal Emission of Gravitational Waves from Weak to Strong Coupling

We study the production of gravitational waves by a thermalized plasma of $\mathcal N$=4 Supersymmetric Yang Mills matter. We focus on the large number of colors limit, $N_c\rightarrow \infty$, and compute the spectrum of gravitational waves both for infinitely large and infinitesimally small values of the 't Hooft coupling constant $\lambda$. In the $\lambda\rightarrow \infty$ limit we employ the gauge/gravity duality to compute the emission rate via the analysis of Energy-Momentum tensor thermal correlators. In the $\lambda \rightarrow 0$ limit we employ state-of-the-art perturbative analyses to calculate the complete leading order emission rate. By comparing these extreme limits, we bracket the magnitude of the spectrum induced by this source of gravitational waves. Embedding our results in a cosmological evolution model, we find qualitative and quantitative similarities between the strong coupling spectrum and the extrapolation of the perturbative results up to an intermediate value of the coupling, after an appropriate re-scaling of the effective number of degrees of freedom. We comment on how our results can help better understand the contribution of thermalized matter to the stochastic spectrum of gravitational waves.


Introduction
The detection of gravitational waves (GW) by the LIGO and VIRGO collaboration [1] has initiated a new era in the exploration of the Universe. The measurements performed by those experiments have had a tremendous impact on our knowledge of compact objects in the Universe and provided tests of general relativity [2]. In addition, gravitational waves at different frequencies carry information of many different cosmological processes at very different energy scales. These probes not only can help us better understand the early stages of the Universe, but can also provide new information on particle physics beyond the Standard Model (SM) [3].
Going beyond our current capabilities for detection, the high-frequency region f ≳ 30 kHz has recently attracted a lot of interest (see for example [4]). Gravitational waves at those frequencies are radiated by several different sources active both in the early and the late Universe (see [5] for a review). While accessing this region requires the development of new technologies, the absence of known astrophysical sources at those frequency bands provides a strong motivation to study gravitational wave production in this challenging part of the spectrum. And one of the unavoidable sources that must be present in this frequency region is thermal emission of gravitational waves.
The fact that fluctuations in a thermalized system lead to the emission of gravitational waves is a well-known phenomenon that can be found in textbooks (see for example [6]). From a quasi-particle point of view, this emission may be understood as the result of the collision of the plasma constituents. In analogy with bremsstrahlung photon radiation, those collisions lead to a gravitational wave spectrum with a characteristic scale given by the temperature T of the plasma. Note that, as long as T is small compared to the Planck mass, T ≪ M Pl , the emitted gravitational waves interact weakly with matter and therefore do not re-scatter with the plasma. Accordingly, the emitted GW radiation does not equilibrate at those temperatures. As a consequence, the spectrum of GW is not thermal, but contains information about the thermalized medium that produced those waves. This source of GW has been recently re-analysed in [7,8] by employing state-of-theart perturbative analysis to compute the complete leading-order GW spectrum within SM. Building on those works, the authors of Ref. [9] performed a survey of theories Beyond the Standard Model (BSM) and argued that this thermal emission may be used to constrain the maximum temperature and the effective number of degrees of freedom of the early Universe. For the SM as well as for the extensions studied in [9], this spectrum peaks in the vicinity of f ∼ 100 GHz and is negligible at the much lower frequencies explored by current or already planed GW experiments.
A natural question to ask is how the features of this spectrum change by increasing orders in perturbation theory. However, going beyond leading-order in the calculation of this rate is a formidable task. As it is well known, the perturbative expansion in thermal field theory exhibits many complications which require performing different types of resummations (see [10] for a recent review). Nevertheless, the next-to-leading order analysis for a similar process (the thermal emission rate of photons by a QCD plasma) has been performed in [11]. The analysis of that photon rate indicates that even for large values of the strong coupling constant, the main qualitative and quantitative features of the leading-order calculation remain the same at next-to-leading order. It would be interesting to complete the analysis for this channel.
A complementary way to assess the effect that matter interactions have on the production of GW is to study how much the spectrum changes when varying the magnitude of the coupling constant. A particularly interesting regime is the large coupling limit, when interactions among the plasma constituents make them lose their individual identity and the thermal excitations of the plasma cannot be described by quasi-particles. To perform such analysis we need to compute the emission rate by other means than perturbation theory. While we do not expect to find such regime within the SM, some of its extensions, for example composite models (see [12] for a review) may possess a strongly-coupled phase.
While the strong coupling limit is challenging in general, there is an interacting gauge theory, N = 4 supersymmetric Yang Mills (SYM), for which we can access this regime. Since this theory is conformal, the 't Hooft coupling λ ≡ g 2 N c does not run and becomes a parameter of the theory which we can choose. In the large λ limit, the so-called gauge/gravity duality [13] provides us with an alternative tool which allows us to address dynamical questions in the large number of colors N c → ∞, large 't Hooft coupling λ → ∞ limits.
In this paper we use the gauge/gravity duality to perform the first calculation of the spectrum of thermally produced gravitational waves at strong coupling. We will contrast this strong coupling calculation with a complete leading-order analysis of the same quantity in perturbative N = 4 SYM. For this perturbative computation we will exploit the perturbative results obtained for a general SU (N ) gauge theory coupled to an arbitrary number of bosons and fermions presented in Refs. [8,9]. This paper is organised as follows. In Section 2, we review the weak-field approach to gravitational waves and the relation between the emission rate and a particular contraction of the thermal Energy-Momentum tensor correlators. In Section 3, we use the gauge/gravity duality to compute this correlator and determine the thermal emission rate of gravitational waves at strong coupling. In Section 4, we perform the perturbative computation and compare it with the strong coupling analog for different values of the coupling constant. We later embed these two computations into a cosmological model in Section 5 to compare the spectra after expansion. Finally, in Section 6, we summarize our main findings and we put our results into the general context of characterising this particular component of the stochastic GW background.

Gravitational waves from a thermal source
In this section, we briefly review the existing literature on the emission of gravitational waves from a thermal source. The main physics of this process may be found in textbooks (see for example [6]). To make a direct connection with the strong coupling calculation of Section 3, we will express the production rate in terms of a particular thermal correlator of the energy-momentum tensor. This discussion is based on [7], which we follow closely.
As it is well known in general relativity, gravitational waves at far-away distances from the emission source may be described under the weak-field approximation by small perturbations of flat spacetime η µν , such that The gauge freedom enjoyed by metric fluctuations can be removed by sticking to some particular gauge. In the so-called harmonic gauge, h µν satisfy the condition Also, physical excitations that propagate far away from the source can be expressed in terms of the transverse traceless (TT) components of the fluctuating fields, h TT µν : with k µ = (ω, k) the four-momentum of the wave. Under these gauge conditions, the spatial Einstein equations become where □ ≡ −∂ 2 t + ∇ 2 and T TT ij are the spatial components of the transverse traceless part of the energy-momentum tensor. After a Fourier transform, T TT µν may be obtained as where we have introduced a projection operator onto the transverse traceless modes with P ij the spatial components of the transverse projector A solution to Eq. (2.4) can be easily found in momentum space as where the iϵ-prescription selects the retarded solution. Fourier transforming the frequency, this solution may be written as Since the h TT ij components possess all the physical information of the gravitational field excitation, the classical energy carried away from the source by gravitational waves may be expressed in terms of h TT ij as Substituting the time-derivative of Eq. (2.9) into this expression yields (2.11) As it is standard in the classical computation of the radiation, the energy carried by gravitational waves is obtained after averaging the instantaneous energy E GW over an observation period that is long as compared to the frequency of the wave: After this observation-time average, the radiated energy is given bȳ This expression is valid for any source; in particular, it is valid for any of the individual configurations of the thermal ensemble. To determine the energy emitted by the thermal system, we average over those configurations with the thermal weight. Assuming that the theory is space and time translation invariant, we express the emitted energy in terms of the correlator where the brackets represent the thermal average and V is the volume of the system. After averaging over all configurations and taking a time-derivative in Eq. (2.13), the mean radiated energy per unit time becomes Introducing the definition of C(k, τ ) given by Eq. (2.14) and using the time translation invariance of the average ensemble, we express the energy density ρ GW of emitted gravitational waves per unit time as where {A, B} ≡ AB+BA. As anticipated, this classical analysis relates the rate of emission of gravitational waves to a correlation function in medium. However, the symmetrized appearance of the correlator is an artifact of the classical approximation, since in this limit the order of operators does not matter. As argued in [7], and in analogy with the thermal emission rate of (weakly-coupled) photons, the full calculation is expressed in terms of a particular Wightman function. Recalling the definition of T TT ij given in Eq. (2.5), and since Λ ijmn is a projector, the energy density of gravitational waves per unit time is given by The expression for the emission rate that we have just derived may be apparently formal, but it has a convenient form for our purposes: it explicitly separates the weaklyinteracting dynamics of the gravitational waves emitted by the thermal plasma from the self-interaction of the matter in the plasma. While the calculation has been performed in the weak-field approximation for the gravitational disturbance, we have made no assumption on the dynamics of the plasma, other than some general symmetries, such as spacetime translation invariance. In the next section we will evaluate the thermal correlator in different approximations to the plasma dynamics.

Thermal emission rate at strong coupling
In this section, we evaluate the thermal correlator appearing in Eq. (2.17), and therefore the emission rate, in a strongly-coupled theory. In particular, we will focus on interacting gauge theories with a large value of the coupling strength. This is a complicated limit to analyse theoretically, since there are not many computational methods available. Certainly, perturbative methods are not generally applicable in this limit. Similarly, lattice field theory computations, even for theories with a well understood lattice implementation, are also of limited guidance, since the emission rate equation (2.17) is controlled by a light-like correlator. For this reason, in this paper we focus on N = 4 SYM at strong 't Hooft coupling λ → ∞ (SCSYM), since for this theory we have a framework with which to address the strongly-coupled regime: the gauge/gravity duality. In this section we use this framework to compute the emission rate.

Holographic analysis of the thermal emission rate
The gauge/gravity duality, also known as holography, is a relation between a quantum field theory and a theory of gravity in at least one extra dimension (see [14] for a review). The earliest example of this duality, the so-called AdS/CFT correspondence [13], establishes a relation between N = 4 SU (N c ) SYM in the large number of colors N c → ∞ and the large 't Hooft coupling λ → ∞ limits, for λ ≡ g 2 YM N c , and type IIB supergravity on an AdS 5 ×S 5 background, whose metric is given by with R the radius of AdS and dΩ 2 5 the metric of a unit five-sphere. For the vacuum of N = 4, the function f (r) = 1. For a thermal system, the function f (r) transforms into f (r) = 1 − r 4 0 /r 4 and Eq. (3.1) therefore accounts for the black brane metric, with r 0 representing the position of the event horizon. In the latter case, the temperature of the dual field theory is the Hawking temperature of the black brane, T = r 0 /πR 2 . Since the observables we describe in this paper are not charged under the R-current, we will not consider excitations of the five-sphere and, therefore, it will suffice to restrict ourselves to the five-dimensional part of the background. Introducing a new radial coordinate u = r 2 0 /r 2 , the metric may be rewritten as According to the dictionary of the gauge/gravity duality, bulk fields in the gravity theory couple to local operators in the dual quantum field theory. More specifically, for each possible operator O(x) and its corresponding point-dependent source ϕ(x), we can assign a dual bulk field Φ(x, u) whose value at the boundary identifies with ϕ(x) [14][15][16].
One of the important implications of this relation between bulk fields and operators is the possibility of calculating the n-point function of local operators in the gauge theory in terms of a gravity description [14][15][16]. For instance, the general strategy for computing the retarded two-point function of a local, gauge-invariant operator O in terms of the gravity prescription is as follows [17]: identify the bulk field Φ dual to the operator O; find the string action and the corresponding equation of motion satisfied by Φ and solve it subject to the in-falling boundary condition; express the solution Φ c in the basis of two local solutions at the boundary; the connection coefficients of this linear combination determine the retarded Green's function for O.
This statement connects directly with the result we have derived in Eq. (2.17): the production rate of gravitational waves from a thermal source depends on the two-point correlation function of the energy-momentum tensor. In this way, the gauge/gravity duality allows us to simplify the calculation of the energy production rate by a thermal ensemble of strongly-coupled matter by simply transforming it into a general relativity problem. We only need to apply the correspondence to the concrete case of the energy-momentum tensor two-point function.
According to the duality, the energy-momentum tensor is dual to the metric of the gravitational theory. Thermal energy-momentum tensor correlators are determined by considering small gravitational perturbations, h µν , over the black brane metric. Given that we are using the same notation, and to avoid possible confusions, it is important to point out at this point that these metric perturbations are different from the gravitational wave fluctuating fields described in Section 2. In the latter case, those are fluctuations over the four-dimensional metric on the field theory side and represent dynamical fields that interact weakly with the strongly-coupled matter; on the former, the metric fluctuations encode modifications of the 5-dimensional metric dual to a state in the field theory and, from the point of view of the field theory, may be viewed as a calculational tool.
Before we describe the dynamics of the metric fluctuations in AdS 5 , let us classify the possible excitations that we may encounter. Without loss of generality, we shall assume these fluctuations to propagate along the z-axis and to be of the form with momentum k = (0, 0, k). With this assumption, metric perturbations can be classified according to their transformation properties under rotations in the xy-plane, whether they transform as vectors (shear channel), as tensors (tensor channel) or remain invariant (sound channel) [18]: with α, β = x, y. These excitations are governed by Einstein equations to linear order in h. Since Einstein equations and the black brane state are both invariant under rotations, the equations obeyed by each of these three sets of fluctuation components decouple from each other, making it easy to evaluate them separately. Each of the symmetry channels has its own contribution to the energy-momentum tensor correlation function. In this way, in a conformal quantum field theory, the retarded two-point function of the energy-momentum tensor in a thermal rotation-invariant state is given by a sum of three independent components [18] with S µναβ , Q µναβ , L µναβ appropriate orthogonal projectors providing three different Lorentz index structures and G 1 , G 2 , G 3 three scalar functions determining the shape of the correlator. Those Lorentz structures may be constructed in terms of the transverse and longitudinal parts of the projector (2.7), defined as The three relevant, mutually orthogonal combinations of P T µν and P L µν are the following: with H µναβ a projector onto traceless symmetric tensors The first combination (3.7a) projects onto the vector modes (shear channel), whereas the second one (3.7b) projects onto the scalar modes (sound channel) and the third one (3.7c) onto the tensor modes (tensor channel). The corresponding scalar functions G 1 , G 2 , G 3 therefore couple to the spin 1, spin 0 and spin 2 channels, respectively. Finding a solution to the second-order ODEs obeyed by the variables in each of these channels allows us to compute the scalar functions G i and thus also the retarded energy-momentum tensor correlation function. However, one should notice that the energy production rate (2.17) we are ultimately interested in does not depend on the retarded correlator, but rather on the spatial components of the Wightman function, which we shall denote as This correlator may be determined by making use of the Kubo-Martin-Schwinger (KMS) relations [19], which allow us to write it in terms of the spectral function which, in turn, is proportional to the imaginary part of the retarded correlator as Using the KMS relations, the Wightman function is related to the spectral function as where n B = 1/ (exp[k/T ] − 1). Writing Eqs. (3.11) and (3.12) together yields the relation between both correlation functions: Introducing this expression into Eq. (2.17) and computing the projectors for k = (0, 0, k), the energy production rate reduces to 14) or, equivalently expressed per unit logarithmic wave number interval, Here we have used that Λ ijmn S ijmn = Λ ijmn Q ijmn = 0 and Λ ijmn L ijmn = 2, which implies that the energy-momentum tensor in a thermal gravitational-wave system only couples to fluctuating fields in the tensor channel (3.4c). For this reason, to compute the GW production rate, it will suffice to restrict ourselves to the study of the tensor channel and keep only off-diagonal perturbations h xy .

Computation of the emission rate
In this section, we study the dynamics of tensor fluctuations over the black brane solution (3.2). This type of fluctuations has been extensively studied in the past (see for example [18]). Here we will follow the general analysis of [20,21] and use it to numerically determine the spectral function ρ, defined in Eq. (3.11), for light-like momenta. The dynamics of tensor fluctuations about the black brane metric are governed by the quadratic approximation to the gravitational action with G 5 the effective five-dimensional Newton constant, Λ ≡ 6/R 2 the cosmological constant of the AdS spacetime and K the extrinsic curvature over the induced boundary metric ∂g. As a consequence of the duality, G 5 is related to the parameters of the field theory as where N c is the number of colors of the dual SU (N c ) gauge theory. Either by extremizing the gravitational action (3.16) or by linearizing the Einstein equations, the equations for fluctuations h µν of the form (3.3) may be written as µν the Ricci tensor to linear order in the perturbation h. As it is well known, diffeomorphism invariance implies that h possesses a gauge symmetry and, therefore, not all solutions to these equations correspond to physical excitations. For this reason, it is convenient to remove this gauge freedom by working with gaugeinvariant variables constructed as linear combinations of the fields and their derivatives in the corresponding symmetry channel (3.4). Since we are only interested in spin 2 variables, we will restrict ourselves to the gauge-invariant variable in the tensor channel [18]: Introducing this expression into the linearized equation of motion (3.18) for µ = x, ν = y, yields the following second-order differential equation in terms of the new variable Z: with the dimensionless parameters w and κ defined as w = ω/2πT and κ = k/2πT , with ω and k the frequency and momentum of the excitation. Since gravitational waves are controlled by light-like correlators, in the rest of the analysis we will set w = κ. The general strategy for computing the two-point correlation function follows the steps we have outlined in Section 3.1. Following that discussion, we search for a solution of this equation satisfying the incoming-wave boundary condition. However, given that it is not possible to find analytic solutions to the equation (3.20), we will solve it through numerical methods by performing a series expansion near the boundary and the horizon.
The above ODE (3.20) has two singular points at u = 0, 1. Following [22], we can work out a Frobenius series solution about each of these singularities. Near the horizon u = 1, the series solution takes the form with F H (u) a power series about u = 1. The minus sign corresponds to the incoming-wave solution, while the plus sign is associated to the outgoing-wave solution. Since we demand Z to satisfy the incoming-wave boundary condition, we will reject the solution with a positive sign in the exponent. The expression for F H (u) can be inferred by demanding Eq. (3.21) to satisfy the equation of motion (3.20). In this way, we may find a power series solution near the horizon: From this approximate solution we can extract a numerical value for Z and its first derivative Z ′ at any point very close to the horizon. These initial conditions, together with the above equation of motion, define an initial value problem. We can then implement a multiple shooting method that allows us to construct a numerical solution over the whole interval between the horizon and the boundary, u ∈ (0, 1). A similar analysis can be done near the boundary, where the two series solutions acquire the form

25)
As before, these solutions are used to define an initial value problem and to implement a multiple shooting method to construct a numerical solution over the whole interval between the horizon and the boundary. Since the ODE (3.20) is second order, the three solutions that we have discussed are not independent. The in-falling solution Z 1 (and its canonical momentum) can be expanded in the basis of the two local solutions at boundary Z 1 , Z 2 , which implies for any point u 0 between the boundary and the horizon. We can therefore use the numerically computed solutions to determine the connection coefficients A and B for different values of κ.
Once A and B are known, the retarded correlator at light-like momentum condition is given by [18] with the constant coming from the prefactor 1/16πG 5 of the action (3.16). However, as discussed in Section 3.1, in order to implement this result into the computation of the gravitational-wave spectrum, we need to apply the KMS relations and express the retarded propagator in terms of the Wightman function. The relation between both thermal correlators allows us to compute the energy production rate per logarithmic frequency interval by simply introducing Eq. (3.28) into Eq. (3.15): for Planck mass M 2 Pl = 1/G, where we have defined a dimensionless function which we shall denote as the source term. The energy spectrum with respect to the dimensionless frequency of the waves is shown in Fig. 1. In next sections, we will extend the discussion about this result and its cosmological implications.

Thermal emission rate at weak coupling
So far in this paper, we have derived the thermal spectrum produced by a strongly-coupled thermal source. In this section we will focus on the opposite limit, that is, small interaction between the plasma constituents λ → 0. Even though perturbative analysis has already been made within thermal emission in different gauge theories [7][8][9], if we want to estimate the effect of the coupling on the shape of the GW spectrum it is important to compare the strong and weak coupling computations in the same theory. For this reason, in this section we will determine the thermal production rate of gravitational waves in N = 4 SYM theory at weak 't Hooft coupling λ → 0 (WCSYM), while keeping the large N c → ∞ approximation. For our perturbative analysis, we will take advantage of the generalization of the perturbative calculations of [7,8] to a theory with arbitrary gauge group, Weyl fermion and scalar content performed in [9]. As it is well known, computing thermal correlators in perturbation theory for all possible values of the external momentum is a complicated task that requires different approximations for different momentum regions. Here, as in [7][8][9], we will only focus on the region where the momentum of the correlator is much larger than the Debye thermal mass, k ≳ m(T ), since it is in this region where the perturbative thermal emission rate is maximal. In this momentum range, the complete leading-order result of [8] for the spin-2 projection of the energy-momentum tensor 1 two point function appearing in Eq. (2.17) was extended in [9] for a general SU (N c ) gauge theory 2 to givê η k T ≡k ≃η HTL (k) where T Ad = N c , T i and T α are the Dynkin indices of the Adjoint representation of the gauge group, and of the individual gauge representations of the scalar and fermion fields, respectively. The four two-loop functions η gg (k), η sg (k), η f g (k), η sf (k), involve diagrams with only gauge fields, scalar and gauge fields, fermions and gauge fields, and scalars and fermions, respectively (their full description and mathematical definition can be found in [8]). In addition, the "hard thermal loop" (HTL) resummation functionη HTL is given byη which depends on the Debye screening length of the plasma, given at leading order by Finally,η depends on the Yukawa couplings among fermions and scalars, which in [9] are parametrized as where ϕ and ψ are the scalar and Weyl fermion fields of the theory. In this contribution of the Lagrangian, all internal indices are traced over.
To obtain the thermal emission rate in N = 4 SYM, we particularize those general results to this precise case. This theory contains four Weyl fermions ψ α , α = 1, . . . , 4, and six real scalar fields Φ i , i = 1, . . . , 6. Both these two types of matter fields are in the adjoint representation of the gauge group, and therefore their Dynkin indices are T i = T α = T Ad = N c . Introducing these properties of the matter fields into the general expression for the Debye mass (4.3), we obtain the leading-order expressionm 2 WCSYM = 2λ, with λ the 't Hooft coupling. This knowledge of the matter content is sufficient to determine all contributions to the production rate in Eq. (4.1) but those proportional to the Yukawa couplings. The N = 4 SYM Lagrangian does not possess a term of the form (4.4); however, it does possess interacting terms involving two fermions and one scalar field. To write this term explicitly, we use the representation of the N = 4 Lagrangian used in [25,26], where the scalar fields are expressed as a multiplet where X p and Y q , p, q = 1, 2, 3, denote scalar and pseudoscalar fields. The interaction Lagrangian of the scalar with two fermions is then where we have introduced a Majorana fermionΨ α ≡ ψ α ,ψ α and A p and B q are 4 × 4 matrices that satisfy A p αγ A p γβ = B p αγ B p γβ = −3δ αβ . The explicit form of these matrices is given in terms of the 2 × 2 Pauli matrices as The main difference between the interaction Lagrangian (4.6) and the Standard-Model-like Yukawa Lagrangian (4.4) lies in the gauge-index structure of the interaction. However, for leading-order computation, the gauge indices in Eq. (4.6) may be regarded as internal indices, since at this order the fields interacting via this vertex do not suffer color exchanges.
Expanding the matter fields in the basis of color generators as Φ i = Φ a i t a and ψ α = ψ a α t a , with t a the SU (N c ) generators in the fundamental representation, 3 we can identify where the couplings y iabc αβ are combinations of the A and B matrices and the structure constants of SU (N c ). Using these relations, we find iαβ |y iabc αβ | 2 = We have already identified the values of all parameters appearing in the N = 4 SYM source term (4.1). Taking the large-N c limit to facilitate the comparison with the strong-coupling calculation, the final leading-order expression for this conformal field theory iŝ Note that the λ-dependence of this expression is non-trivial. The two-loop functions η gg , η sg , η f g and η sf depend solely on the dimensionless ratiok = k/T and their contribution to the correlator is proportional to λ, as expected from a leading-order calculation. However, this is not the case for theη HTL function. The reason for this non-trivial dependence on the coupling is that this contribution resums a set of higher-order diagrams that become large for momentak ∼m, via the HTL approximation. Note also that an important assumption of the HTL framework is that the temperature and the Debye mass scales are well separated, i.e., m ≪ T . We therefore expect this result to hold for the λ ≪ 1/2 regime in N = 4 SYM theory. Now that we have introduced all the features and quantities characterizing our theory, we can derive the spectrum of stochastic thermal waves emitted within weakly-coupled N = 4 SYM. Analogously to Eq. (3.29), the energy density carried by gravitational waves per unit time and unit logarithmic wave number interval is given by As we can see, the source term and, consequently, the shape of the spectrum emitted by a weakly-coupled thermal plasma in SYM depend on the coupling constant λ of the theory. In order to visualize this coupling-dependence of the emission rate, the full spectra of thermally produced gravitational waves in the perturbative regime is shown in Fig. 2 for convenient values of λ: λ = 0.1 (dotted green), λ = 0.5 (dashed orange) and λ = 1.8 (dotdashed blue). Note that only the first one of these values satisfiesm < 1 and may therefore lie within the regime of validity of this calculation. Based on these criteria, the second value λ = 0.5 marks in fact the limit beyond which this analysis is no longer valid. Nevertheless, it will also be useful to extrapolate the perturbative calculation to intermediate values of λ in order to compare its behaviour with the strong coupling computation.
To better understand the main features of the perturbative calculation, we study the λ-dependence of two key characteristics of the spectrum: the value of the spectrum at the extrema (peak energy density) and the location of the latter. These two quantities are shown in Fig. 3. As we can see in the left panel, the position of the peak increases monotonically with the coupling. On the contrary, the behaviour of the maximum of the spectrum is not monotonic: at small coupling, where the perturbative computation is valid, the energy radiated in gravitational waves increases as the coupling grows. This behaviour is expected, since as the coupling grows the plasma constituents suffer more violent collisions, leading to an increase in the emitted radiation. However, when we extrapolate the perturbative calculation to λ ∼ 1, the maximum of the spectrum starts to decrease, as shown in the right panel of this figure. This maximum in fact occurs at λ ≃ 1.8; the leading-order perturbative spectrum extrapolated to this value is shown in the dot-dashed blue curve of Fig. 2. We will come back to this non-monotonic behaviour when comparing this calculation with the strong-coupling analysis of Section 3.

Cosmological evolution
The analysis we have performed so far within the framework of N = 4 SYM theory has been restricted to the emission from a static thermal source. However, in any conceivable situation in which the thermally-emitted gravitational waves might be detected, matter does not remain static but expands with the expansion of the Universe. Certainly, no form of matter in the Universe is governed by N = 4 SYM. Nevertheless, if we want to use Figure 2: Comparison between the energy density of gravitational waves produced by quasi-particle excitations in a weakly-coupled thermal plasma for different values of the coupling constant λ. The lines are colored for the values of the wave number in which the computation can be trusted (k > 2λ). The spectra correspond to λ = 0.1 (dotted green), λ = 0.5 (dashed orange) and λ = 1.8 (dot-dashed blue). The peaks of the spectra can be seen to occur around k ∼ T for different values of λ. Figure 3: Left: position of the peak wave number with respect to the coupling constant λ; the perturbative analysis predicts that an increase in the interaction strength of the thermal plasma shifts the peak energy towards higher frequencies. Right: peak energy density with respect to the coupling constant λ; contrary tok peak , the peak energy does not increase monotonically with λ, but rather increases until reaching a maximum value at λ max ∼ 1.8, then decreases. the N = 4 calculations to build intuition about the coupling-dependence of the emission rate, we must include the cosmological history in our analysis, since this process strongly modifies the magnitude and position of the gravitational-wave signal.
For this study, we will assume that over a certain time in the early, very hot history of the Universe, matter dynamics could be understood in terms of N = 4 SYM interactions. We will also assume this period to extend from a very early time t 0 when the temperature of the Universe was T max to a later time t end with much lower temperature T end ≪ T max . Later than that, matter dynamics stops being close to N = 4 SYM. Even though matter may still be interacting below T end , we will assume the output of gravitational waves to be negligible for times larger than t end . This assumption is justified since the spectrum of gravitational waves scales with the temperature of the plasma [9]. Therefore, the energy emitted by thermalized matter is dominated by the range T end < T < T max . In this region of temperatures, the expansion of the Universe is also taken to be dominated by the energy density of the thermalized matter e(T ), that is, This expansion rate implies that the temperature of the system decreases with time. Since the theory we are considering is conformal, the ratio T (t)/a(t) remains constant, leading to the following relation between time and temperature: This expression may be further simplified by exploiting the conformal invariance. Following the notation in Ref. [9], we define the effective number of degrees of freedom of the theory as where s accounts for the entropy density. In general, g * s (T ) is a temperature-dependent function. However, for the particular case of a conformal theory, it becomes a couplingdependent constant. Conformal invariance also implies that the energy density is given by e = 3g * s T 4 /4. Therefore, in a conformal theory, the relation (5.2) becomes For N = 4 SYM in the large-N c limit, the effective number of degrees of freedom in the λ → 0 limit may be obtained directly from the matter content of the theory and is given by 4 g * s,WC = 15N 2 c . Meanwhile, for the λ → ∞ limit, the holographic computation of the entropy [27] yields g * s,SC = 45 4 N 2 c . The cosmological expansion characterized by H induces a dilution of the energy density carried by gravitational waves. As long as H ≪ T , we can identify the microscopic calculation at fixed temperature as a rate of energy injection per unit volume to this component of the energy density. As the Universe cools down, this rate decreases as time progresses.
To simplify the discussion, in this section we will denote by R the energy spectrum of gravitational waves per unit three momentum which, as we know, is a function of the coupling constant in N = 4 SYM. Compactifying Eqs. (3.29) and (4.11) into one single expression: The energy injected by this rate adds to the total energy density of gravitational waves, which, taking into account the expansion rate, satisfies the differential equation [7] ( This differential equation may be easily integrated from t 0 to t end . This way, we can express the energy density of gravitational waves at t end as where in the second line we have used the relation (5.4) and the definition of the effective number of degrees of freedom (5.3).
In addition to diluting the energy density, the expansion of the Universe also changes the frequency of the emitted waves. Therefore, to determine the spectrum at a fixed value of k end , the wave momentum at t end , and since momenta redshift as k(t) = k end a(t end )/a(t), we must evaluate the emission rate at different values as time changes. However, in a conformal theory, the ratio k/T remains constant during the expansion. Taking this into account, the GW energy density per logarithmic wave number interval at t end turns into where we have introduced the identity (5.3). Using Eq. (5.5), we observe that the integrand in this expression does not depend on T , and the integral can be easily performed. Assuming T max ≫ T end we neglect the lower limit of integration and express Although we assume that the production of gravitational waves below T end is negligible, the expansion of the Universe continues redshifting and diluting the spectrum. The effective number of degrees of freedom also changes after t end , since at the end of the evolution the Universe contains the current matter content. We will assume that these changes are fast and are followed by long periods in which the equation of state is approximately conformal. With this picture in mind, we can approximate the dilution of the spectrum from t end up to the final time t fin by ρ(fin) ≃ ρ end a 4 (t end )/a 4 (fin), which implies where g * s (fin) = 3.931 ± 0.004 [28] is the number of entropy degrees of freedom after neutrino decoupling and Ω γ ≡ e γ (fin) e cr (fin) = 2π 2 30 is the present energy density of the CMB photons, with temperature T 0 = 2.72548(57) K. Finally, we need to express the spectrum in terms of the present-day frequency. Assuming once again close-to-conformal expansion, the current-day frequency of the waves is related to k end as leading to a suitable expression Using these relations, the present-day spectrum is given by Substituting either the strong coupling (3.30) or the weak coupling source term (4.10) into this expression, we can draw a numerical estimate of the gravitational waves energy production rate in the context of an expanding cosmological model for the different values of the coupling constant. Note that both in the weak and strong coupling calculations the number of effective degrees of freedom is proportional to N 2 c . The magnitude and the position of the peak energy density are thus expected to depend on the number of colours of the theory, which we do not specify. This sets an additional degree of freedom in the shape of the spectrum. Finally, one should also notice that g * s at strong and weak coupling differ by a 3/4 factor, which shifts the spectrum to somewhat different frequencies as a consequence of the expansion.
The result of these two extreme calculations are shown in the top panel of Fig. 4. In this plot, the strong coupling result corresponds to the solid black line while the weak coupling calculations for λ = 0.1, 0.5, 1.8 correspond to the dotted green, dashed orange and dot-dashed blue, respectively. Note that, for all the values of the coupling, the strong coupling result has a larger magnitude in the spectrum of gravitational waves. Since, as we have already seen, λ = 1.8 corresponds to the maximum value of the peak energy density we can obtain from the extrapolation of the perturbative, leading-order result, we can ensure that the strong coupling calculation produces more gravitational waves than the extrapolated perturbative calculation for any value of the coupling. Note, however, that the main features of the strong coupling spectrum are well-captured by the maximum perturbative result at λ = 1.8.
This qualitative agreement between the strong coupling and λ = 1.8 computations is, in fact, the result of two different contributions. On the one hand, the thermal sourcê η at the peak position is larger for the perturbative result than for the strong coupling calculation, as it may be inferred from inspection of Fig. 1 and Fig. 2. On the other hand, the effective number of degrees of freedom is reduced at strong coupling. This enhances the strong coupling spectrum since, at equal temperature, the smaller number of degrees of freedom slows the expansion and leads to a smaller dilution of the signal.
To isolate the effect of the microscopic emission rate, in the lower panel of Fig. 4 we have scaled-out the number of degrees of freedom from both the frequency and the magnitude of the spectrum. By doing this scaling, the final spectra become a redshift of the thermal spectra of Fig. 1 and Fig. 2 together with a reduction of the amplitude common to all calculations. After this scaling, the extrapolated λ = 1.8 result is no longer a good approximation to the strong coupling result. However, after the degrees of freedom are factored out, we find a surprising semi-quantitative agreement between the strong coupling result and the perturbative extrapolation to λ = 0.5. Note that this value marks a limit of applicability of the perturbative analysis, since beyond this value the Debye mass is not separated from the T scale. This means that λ = 0.5 may be used as a good predictor of the strong coupling limit, at least in a static thermal system.

Discussion
The spectrum of thermally-emitted gravitational waves is an interesting observable which could potentially provide us with valuable information of the very hot stages of the Universe, which are hardly accessible to other probes. Nevertheless, this is also a very challenging observable, since the region of frequencies where the spectrum is maximal, in the vicinity of 100 GHz, is not within our current detection capabilities. Nevertheless, the interest in high-frequency gravitational waves is pushing the community to devise new detection technologies. Therefore, it is worth analyzing this mechanism of GW production in depth.
In this paper we have used N = 4 SYM as a tool with which to understand this process. Certainly, this theory is very different from the Standard Model and modelling the matter present in the universe with this theory is in no possible sense an approximation. However, this theory allows us to determine the gravitational-wave spectrum in very different dynamical regimes. In particular, and thanks to the gauge/gravity duality, we can access the plasma dynamics at strong coupling, when the plasma cannot be considered as a collection of interacting quasi-particles, in contrast to perturbation theory. This fact implies that calculations based on N = 4 may be used not only to get an intuition on how the emission rate depends on a particular parameter of the theory, but also whether this observable is sensitive to the presence or absence of quasi-particles in the plasma. There is little doubt that unbroken phase of the Standard Model is dominated by quasi-particles. However, thermal dynamics at very large scales or even in the dark matter sector, if it behaves thermally at some early stage of the evolution of the Universe, could have a strongly-coupled phase with no quasi-particles. The strategy to use the gauge/gravity duality to describe the dynamics at strong coupling has been very successful in understanding the properties of the QCD plasma close to the deconfinement transition (see [29] for an extended review and [30] for an analysis of thermal photon and dilepton emission, which are weakly-coupled probes of the QCD plasma, similar to the one described in our work). For other applications of the gauge/gravity duality to understand different processes of gravitational-wave emissions see [31][32][33][34][35][36][37].
In this work, we have focused on the spectrum of gravitational waves produced by a thermal system at momenta in the vicinity of the maximum of the emission rate. At very small momenta, the emission rate is solely controlled by the shear viscosity of the system [7]. Since viscosity is dominated by the least interacting degrees of freedom [38], we have not explored this region with our strong coupling computation. In addition, the rate of gravitational-wave emission is extremely small for those small frequencies, much smaller than other known sources of stochastic gravitational waves [7][8][9]. On the contrary, the maximum of the spectrum is more sensitive to interacting degrees of freedom. As our calculation shows, the qualitative features of the thermally-produced gravitationalwave spectrum are very similar for the two extreme microscopic pictures of the plasma constituents described by the λ → 0 and λ → ∞ limits of N = 4 SYM. Nevertheless, there are quantitative differences between those two extremes.
For very weakly-coupled plasmas, the emission rate of gravitational waves is suppressed by the smallness of the coupling, since this radiation cannot be produced in the absence of matter interaction. Naively, the perturbative spectrum is proportional to λ, but the complicated infrared dynamics of gauge theory plasmas alters this behaviour as a consequence of the emergence of the Debye mass scale. As the coupling increases, so does the emission rate. However, it is hard to assess whether this growth should be monotonous with the coupling constant. Nevertheless, the strong coupling computation provides a natural benchmark to estimate the maximum strength of the emission rate, at least for N = 4 SYM. In this paper, we have compared this limit with extrapolation of the leading-order perturbative calculation to intermediate values of the coupling. Even though next-to-leading order correction should be significant for the large values of the coupling we have considered, these results provide us with a second reference to assess the spectrum.
Remarkably, the extrapolation of the perturbative spectrum to intermediate coupling is not monotonous but rather possesses a maximum at λ ≃ 1.8. Because of this, leadingorder perturbative computations in this region show very little sensitivity to the particular value of the coupling. At fixed temperature, this maximum of the perturbative emission rate leads to larger output of gravitational waves than the strong coupling calculation at comparable momenta. However, if the fixed temperature computations are embedded into a cosmological expansion model, these two calculations yield a comparable spectrum of gravitational waves at present time. Nevertheless, this observation depends significantly on the change of the effective number of degrees of freedom from weak to strong coupling.
For close-to-conformal evolution, this dependence on the number of degrees of freedom can be factored out in such a way that the final spectrum becomes a shift of the static thermal spectrum. By comparing the weak and strong coupling calculations of the rescaled spectrum or, equivalently, the fixed temperature spectrum, we find that the weak coupling calculation λ = 0.5 result reproduces even quantitatively the strong coupling result. This is an interesting observation, given that λ = 0.5 marks a limit of applicability of the perturbative computation, since for this value the Debye mass equals the temperature. While this may just be a numerical coincidence, it is worth mentioning that at values of the coupling as large as λ = 0.5 the matter fields of the theory acquire a large thermal width, which obscures their interpretation as quasi-particle excitations. It would be interesting to test whether the strong coupling limit may be mimicked by leading-order results at some fixed value of the coupling constant for other gauge theories in which the strong coupling limit can be addressed via the gauge/gravity duality, like for example [39][40][41][42].