Soft gravitational radiation from ultra-relativistic collisions at sub- and sub-sub-leading order

Using soft-graviton theorems a well-known zero-frequency limit (ZFL) for the gravitational radiation flux $dE^{GW}/d \omega$ is re-derived and extended to order ${{\cal O}}(\omega)$ and ${{\cal O}}(\omega^2)$ for arbitrary massless multi-particle collisions. The (angle-integrated, unpolarized) ${{\cal O}}(\omega)$ correction to the flux turns out to be absent in the case of two-particle elastic collisions. The ${{\cal O}}(\omega^2)$ correction is instead non-vanishing and takes a simple general expression which is then applied to bremsstrahlung from two-particle elastic collisions. For a tree-level process the outcome is finite and consistent with expectations. Instead, if the tree-level form of the soft theorems is used at sub-sub-leading order even when the elastic amplitude needs an all-loop (eikonal) resummation, an unphysical infrared singularity occurs. Its origin can be traced to the infinite Coulomb phase of gravitational scattering in four dimensions. We briefly discuss how to get rid, in principle, of the unwanted divergences and indicate --without carrying out-- a possible procedure to find the proper correction to the naive soft theorems. Nevertheless, if a simple recipe recently proposed for handling these divergences is adopted, we find surprisingly good agreement with results obtained independently via the eikonal approach to transplanckian-energy scattering at large (small) impact parameter (deflection angle), where such Coulomb divergences explicitly cancel out.


Introduction
A well known result in gravitational bremsstrahlung is a general formula [1] for the (finite) zero-frequency limit (ZFL) of the spectrum of gravitational waves dE GW /dω emitted in a generic process as ω → 0. Such a formula can be derived either by classical considerations or by using the leading-order soft-graviton theorems going back to the classic paper of Weinberg's [2] and others before [3,4] and after [5][6][7][8] him.
At the same time, new techniques have been developed to compute graviton amplitudes through their connection to gauge theory amplitudes [35,36]. While the main thrust in JHEP05(2019)050 those papers has been towards understanding the ultraviolet structure of quantum gravity or supergravity some attention has also been devoted to the infrared (IR) divergences and their cancellation in sufficiently inclusive quantities.
Although at the time of Weinberg's paper these problems looked highly academic, the recent direct observations of GW from BH and NS mergers [37] can make the issue of gravitational radiation in highly energetic collisions relevant even outside the merger regimes (see e.g ref. [38]).
An obvious question to ask is whether the new developments in soft theorems can be used to extend the predictions about dE GW /dω away from the ZFL. This program has been recently undertaken by Laddha and Sen [39] with very encouraging results. Because of the IR problems with gravity in D = 4 most of the work in [39] has focussed on D > 4 but, very recently an extension to the physically relevant case has been attempted by Sahoo and Sen [40].
In an unrelated development the problem of high-energy gravitational scattering has been studied, particularly since the late eighties. In that case the original motivations were somewhat different: studying gravitational scattering at transplanckian energies [41][42][43][44] (see also [45][46][47]) can shed light on the information puzzle in regimes in which, at least classically, a black hole should form [48][49][50][51] and then seen evaporate via the Hawking process [52]. Another motivation, in the context of string theory, was to study new effects due to the finite size of strings such as tidal excitations [41,44,53] or possible short-distance modifications of gravity [54,55].
That program is still ongoing: in spite of some definite progress (see e.g. [56][57][58] and references therein) constructing a unitary S-matrix in the black hole formation regime, whereby the impact parameter b is smaller than the Schwarzschild radius R S = 2G √ s at the given CM energy E CM = √ s, has resisted so far all attempts. For this reason the attention shifted, at least momentarily, to the study of gravitational bremsstrahlung, an unavoidable phenomenon that should already occur well before the collapse regime is attained. Indeed, in the ultra-relativistic regime, gravitational radiation occurs already at next to leading order in the expansion of the eikonal phase δ(s, b) in Einstein's deflection angle θ E = 2R S /b. This is why the bremsstrahlung problem was first addressed at leading (i.e. θ 2 E ) order. Two very distinct methods have been followed: a classical one [59], using Huygens principle in the Fraunhofer approximation and a second one [60] (see also [61,62]) following the semiclassical approximations implicit in the Amati-Ciafaloni-Veneziano (ACV) approach. Amazingly, the two approaches gave the same result for dE GW /dω as long as ω √ s. The features of the spectrum are quite interesting: they are consistent with the ZFL for ω < 1/b, deviate only logarithmically from the ZFL constant for 1/b < ω < 1/R S , and finally drop above the "Hawking frequency" 1/R S . This intriguing inverse proportionality between the initial transplanckian energy √ s and the characteristic energy of the final gravitons /R s is also the leitmotif of other studies [63,64]. 1

JHEP05(2019)050
The basic motivation of this work is to build a bridge between the soft theorem (valid at sufficiently small ω but for arbitrary kinematics) and the results of [59][60][61]65] (valid at small-deflection angle but presumably in a wider frequency range) and check their mutual consistency in the overlap of the two regimes where both of them can be trusted.
To this end we will extend the calculation of the ZFL to the next two sub-leading terms by taking advantage of the known sub-leading corrections to the soft graviton theorems. Since in this paper we consider directly the physical case of 4-dimensional spacetime, we will have to face some infrared problems that would be absent in higher dimensions. Let us start therefore by outlining the assumptions of our approach and their possible limitations.
Consider a generic process i → f in which we have excluded both real and virtual softgraviton processes/corrections. 2 Let us denote by S (0) f i the bare (uncorrected) S-matrix for the process: (1.1) We will assume that, to the desired order in the low-energy expansion, the inclusion of virtual and real soft-gravitons dresses the bare S-matrix with a coherent-state operator In the last term of (1.2), which specifies the definition of the coherent state, a q , a † q are destruction and creation operators for a (soft) graviton of momentum q and maximal energy Λ. Λ E (with E the characteristic energy-scale of the process) defines what one means by "soft", λ is an infrared cutoff, and λ q is a process-dependent "function" of q to be discussed below.
If the initial state contains no soft gravitons the cross section for emitting any number of soft-gravitons (sgr) with a maximal total energy ∆E (basically identified with the energy resolution of the detector) will exhibit the well-known cancellation of virtual and real soft divergences in the form: with a characteristic dependence on ∆E. We have assumed that ∆E is negligible w.r.t. the total energy E of the process (otherwise there are correction terms, see [2]).
Let us now consider the expectation value of the energy carried by the soft-gravitons in the process at hand. It will be given by the expectation value (in the appropriate coherent JHEP05(2019)050 state) of the gravitational-energy operator d 3 q ωa † q a q : where ω ≡ |q|, the vacuum states refer to the soft-graviton Fock space, and we have used a well-known property of coherent states. Equation (1.4) is the basic equation to be exploited in order to find the connection between soft-graviton theorems and the spectrum of gravitational radiation from a given process. Indeed, factoring out the probability for i → f +soft gr., it follows from (1.4) that: where we have emphasized that λ q depends on the process i → f under consideration. 3 Expanding λ q in power of q will give the corresponding expansion of the GW spectrum. The leading term in |λ (i→f ) q | 2 goes like ω −2 which, given the phase space volume d 3 q, leads to a constant dE GW /dω, the above mentioned ZFL. Naively, the two next-to-leading corrections will correspond to terms in dE GW /dω going to zero as ω and ω 2 up to possible logarithmic enhancements to be discussed later.
There is, however, a subtlety concerning the sub-leading contributions to λ q which can already be seen at tree level. Unlike the leading term the non-leading ones contain differential operators acting on the bare amplitude. When we eventually take the square modulus of the graviton amplitude we have to be careful about which bare amplitude (S or S † ) these operators act on.
A more important issue is what happens to these differential operators when we replace the bare amplitude with the one dressed by the insertion of soft gravitons. It is well known that, in D = 4, the naive operators should suffer modifications because of IR divergences [22]. While the final answer should be free of such divergences one expects finite correction for the sub and sub-sub-leading terms.
There are actually two kinds of IR divergences that can induce finite modifications to the naive recipe. The first is directly related to the cancellation of real and virtual divergences we have already mentioned. It has been discussed, in particular, in [22]. Once the cancellation is achieved some new finite contributions may be left. To our knowledge, no systematic analysis of these finite corrections has been given so far. 4 In section 6 we 3 Our procedure should be justified if ∆E E but much larger than the typical energy ω of the gravitons whose spectrum we wish to compute. A more rigorous framework should be provided by the Kulish-Faddeev prescription [8]. Our procedure basically amounts to assuming that such a prescription extends to gravity up to the sub-sub-leading term. 4 One of us (GV) would like to thank Z. Bern for an interesting discussion about this point.

JHEP05(2019)050
will argue that, under the already mentioned condition ∆E ω, these corrections should be unimportant.
The second kind of divergences is related to the infinite Coulomb phase due to the longrange exchange of massless gravitons. This can also be seen as a divergent time delay in D = 4 which might be cured by properly defining the asymptotic states as in [8]. This case has been recently addressed in a series of interesting papers by Laddha and Sen [25,39] and by Sahoo and Sen [40] and some recipe about how to deal with them has been proposed.
In this case the infinity should be unobservable and simply disappears from any physical observable (such as the energy flux). We shall see, however, that the naive soft theorems at sub-leading level do not eliminate the divergence. An appropriate modification of the soft theorems will remove it but, as pointed out in [39] may also leave behind some finite (logarithmic) corrections. Their existence was recently confirmed in [65] within the eikonal approach.
For the time being we will proceed as if the naive prescription suffers no modification even at D = 4. The need and possible form of the modifications will be discussed in section 6.
The emission of soft gravitons at tree level is governed, at the first three leading orders, by the universal behavior [14,15,27] where κ 2 = 8πG, q µ and h µν denote the momentum and the polarization of the soft graviton, while J µν i = p µ i ∂/∂p i ν − p ν i ∂/∂p i µ denotes the angular momentum operator 5 of the ith 'hard' particle. In this paper, for simplicity, we consider the collision of spin-less particles for which an additional spin term in J µν i is absent. We will also take the massless limit, which is theoretically interesting, should be free of collinear divergences and, as a result, should apply to the case of ultra-relativistic collisions.
An important property of (1.6) is its gauge invariance i.e. the fact that it is invariant under the replacement: which can be easily checked to hold thanks to conservation of linear and angular momentum. It is straightforward to connect the "function" λ q of (1.2) to the soft operators S i of (1.6). This allows to compute the graviton spectrum (and thus the GW energy spectrum) to leading, next to leading and next to next to leading accuracy. Although we shall start with the fully differential distribution, in this paper we will consider the frequency spectrum after integration over the direction of the emitted gravitational wave. We will first review the leading order computation, and then consider the sub-leading corrections. In the rest 5 To compare with the literature note that we have omitted a factor i in the definition of J. All our momenta, including the soft graviton's, are taken to be incoming. We also use the − + ++ Minkowski metric. Finally, one should not confuse the soft operators Si in (1.6) with the S-matrix introduced in eqs. (1.1)-(1.2) and denoted by S.

JHEP05(2019)050
of the paper we shall assume that eq. (1.6) continues to be valid beyond tree level. Such an assumption (for the non-leading term) will have to be checked case by case.
The plan of the paper is as follows. In section 2, we will rederive the zero-frequency limit of the GW spectrum and point out the possible presence, in Weinberg's soft theorem, of sub-leading terms. In section 3, we present the general computation of the O(ω) correction to the GW spectrum. Sub-section 3.1 deals with the general case while in Subsections 3.2 and 3.3 we will show, by two different methods, that the O(ω) correction to the two particle elastic collision is absent. In section 4, we present the O(ω 2 ) corrections to the GW spectrum in the general case. In section 5.1 we specialize again to the case of 2-body scattering and then we apply it to two particular cases: in subsection 5.2 to a tree-level gravitational scattering and, in subsections 5.3, to the leading eikonal resummation of trans-planckian scattering at small deflection angle pointing out a possible infrared problem for the latter case. In section 6 we discuss the issue of IR divergences, their elimination, and the possible finite terms that originate from them. In section 7, we conclude and outline some directions for future investigation. Some details of the computations are relegated to appendices A and B. Appendix C contains an ab-initio calculation of graviton emission from the tree-level gravitational scattering of scalar particles in N = 8 supergravity, analogous to the process discussed in subsection 5.2, showing complete agreement between the two calculations.

The GW spectrum at leading order
To leading order, the emission of real soft gravitons (gravi-strahlung) typically contributes a factor (∆E/λ) B 0 to the inclusive cross section, where ∆E is the maximal energy allowed in soft gravitons and λ is an infrared cutoff.

JHEP05(2019)050
Integration over the soft (final) light-like momentum −q = (|q|, −q) = |q|(1, −n) produces the well-known infrared divergent result [2,64]: where Λ is a characteristic energy-scale of the process. Note that the mass scale µ can be chosen at will since log µ drops out thanks to momentum conservation to leading order in q.
The infrared divergence appearing in (2.4) corresponds to a bremsstrahlung spectrum for the number density dN 0 /dω given by and thus to an energy spectrum: which is nothing but the constant ZFL of the gravitational radiation associated with the specific process under consideration.
Noting that in the sum over i = j each pair i, j is counted twice, we can replace −(p i p j ) by −(p i + p j ) 2 i.e. by the Mandelstam variable associated with each distinct particle pair. In the special case of a 2 → 2 process an extra factor of two follows from the fact that two distinct pairs share the same Mandelstam variable (s, t, or u). The final formula in that case reads: In the small-t (deflection angle θ s ) limit this gives: in agreement with recent classical and quantum calculations [59][60][61].
In the next two sections we will extend the result (2.6) to the first two next to leading orders in ω. However, in order to have the complete result to sub-leading order one has to take into account possible sub-leading contributions already contained in the above calculation. Indeed, to subleading order the momenta p i appearing on the l.h.s. of (2.1) differ from those appearing on the r.h.s. In (2.4), (2.5), (2.6) we have also neglected higher order terms: in particular, if momentum conservation were written as i p i = −q, one would get, instead of (2.3): An even subtler point is that sub-leading corrections to eq. (2.5) depend on exactly which variables are kept fixed while one integrates over angles. Physically, we certainly want

JHEP05(2019)050
to keep the center-of-mass energy √ s 12 fixed and, by fixing ω, also s 34 is kept constant.
However, there is some ambiguity on the 3rd variable one wishes to keep fixed. 6 We will come back to such terms after discussing the other non-leading corrections stemming from (1.6). 3 The O(ω) correction to the GW spectrum The first sub-leading correction arises from the interference between S 0 and S 1 of eq. (1.6) After some tedious but straightforward algebra, the sum over helicities of the emitted graviton produces where Poincaré invariance has been taken into account to set i p i = i J i = 0 at intermediate steps. The resulting expression for B 1 reads where the arrows on J i , J j indicate whether the derivative acts on M N or on its complex conjugate. Clearly the basic integral to compute is: Unfortunately, for a given pair i, j this integral has collinear divergences (that cannot be cured by inserting a lower limit on the soft-graviton's energy). This can be easily seen in the fact that qp i goes like θ 2 qi when the angle θ qi between q and p i goes to zero, while phase space (in four dimensions) gives dθ qi θ qi . This problem disappears after summing over i and j but introducing a cut-off (meaning in this case a mass for the hard quanta) is quite cumbersome. A better way to proceed is to modify I µ ij in such a way as to make it collinear safe for each i, j and µ. We thus make the replacement The numerator now vanishes if q is parallel to either p i or p j , removing the collinear singularities of each term. Furthermore, one can check that the additional terms drop out after summing over i and j thanks to linear and angular momentum conservation.

JHEP05(2019)050
At this point the calculation can be performed in two ways: by splitting the threedimensional integral and carrying out explicitly the angular integration, or by introducing an extra Lorentz invariant δ-function to fix the GW frequency in a conveniently chosen Lorentz frame. In the first approach manifest Lorentz invariance is lost while in the second it is kept. We have checked that both approaches lead to the same final result. Hereafter we follow the second, more elegant method.

Covariant calculation of B 1
In order to arrive at a Lorentz-invariant frequency spectrum let us define the following Lorentz covariant four-vector integral: where P is, a priori, an arbitrary four vector and Λ is a constant with dimension of energy. Clearly K µ ij has nice transformation properties under Lorentz. A physically interesting choice, adopted hereafter, consists in identifying P with the total momentum of the n incoming (or m outgoing) particles in a generic n → m process. In this case Λ takes the meaning of √ s ω 0 where s = −P 2 is the Mandelstam variable of the corresponding channel and ω 0 is the center-of-mass frequency at which we wish to compute the spectrum. The Lorentz-invariant (graviton number) spectrum dB 1 dω 0 is then given in terms of an integral of the type (3.6), i.e.: from which one can compute the GW energy spectrum upon multiplication by ω 0 . Lorentz covariance allows us to expand K µ ij in the form: where K P , K i , K j are functions of the four Lorentz-invariants P 2 = −s, (P p i ), (P p j ) and (p i p j ). By its definition (3.6) K µ ij is orthogonal to both p µ i and p µ j and therefore must be of the form: Contracting K µ ij with P µ and using the already known integral appearing in (2.4): as well as

JHEP05(2019)050
determines the scalar K to be: , (3.12) where, also for future use, we have introduced the convenient quantity: Note the absence of singularities in (3.12) whens ij vanishes. Inserting this result in (3.7) we get (after renaming ω 0 as ω): (3.14) Let us elaborate further the result (3.14) in the general case. Since the vector Q µ ij is orthogonal to both p i and p j we can replace ∂p j ) µ and obtain, after multiplying by ω, We see no reason for why this expression should vanish for a generic n → m process and explicit calculations seem to confirm this. We will see however that (3.15) leads to a vanishing result in the case of a two-body collision. We will show this in two different ways: by working in a special Lorentz frame (the so-called Breit frame), and by using an appropriate Lorentz-covariant definition of the derivatives.

Breit-frame argument for B 1 = 0 in two-body scattering
In order to satisfy all kinematic constraints one can choose the Breit frame (BF) where In this frame the Mandelstam invariants read and obviously satisfy s + t + u = 0. The only subtlety that we have to deal with is how to define d/dp iz in order to preserve p · k = 0, where p = (p, 0, 0) is the 'longitudinal' 3-momentum and k = (0, k, 0) is the 'transverse' 3-momentum. For B 1 this will turn out to be irrelevant since d/dp i,z combine with p i,z = 0 in the BF.
Neglecting an overall factor, it is convenient to rewrite dB
The relevant combinations are Combining the various terms one finds (reg indicates regularisation by a small mass) An interesting feature of the above result is that only combinations of derivatives of the Since these combinations vanish when applied to the constraint E 2 − p 2 − k 2 = 0 their action does not depend on which independent variables we use to express the amplitude.
We now note that dB by simply exchanging the sense of the arrows which amounts to simply changing an overall sign. Therefore dB  In the case of an elastic 2 → 2 process things simplify. Recalling that now P = p 1 + p 2 = −p 3 − p 4 , let us first notice that the contributions with i = j vanish trivially. Also the contributions with i, j = 1, 2 and i, j = 3, 4 vanish since the vectors Q µ 12 , Q µ 34 vanish identically.
We are thus led to consider the remaining pairs: i, j = 1, 3, i, j = 2, 4 and i, j = 2, 3, i, j = 1, 4. These need special attention since it is not a priori obvious how the partial derivatives are defined since the N (= 4) momenta are constrained by momentum conservation p i = 0 and by the mass-shell conditions p 2 i = 0. In terms of the Mandelstam variables s, t, u of the four-point amplitude the constraint s + t + u = 0 holds.
By studying carefully how the soft theorems work at the level of the five-point function (see e.g. the one in appendix C) with a soft, but finite momentum graviton, one can argue that the correct way to define the derivatives is to first replace the Mandelstam variables as follows: and by then letting the derivatives act on the modified Mandelstam variables as if all the momenta were independent (in a sense they are, because in the five-point function . With these rules the basic derivatives become: As a consequence one can easily check that: as one would naively expect. The above rules are also consistent with angular momentum conservation (which we have used already). Indeed: and, after adding to (3.26) J 2 , J 3 and J 4 , one easily finds: These partial-derivative rules will be taken up again when discussing the O(ω 2 ) correction.

JHEP05(2019)050
For instance, consider the contribution from (i, j) = (1, 3) to which we must add, of course, the one from (i, j) = (3, 1). As a result, the sum − → ∂ 1 + − → ∂ 3 appears. Because of the relative signs with which p 1 and p 3 appear in (3.24) only the term ∆ µ u ∂ u will survive. However, Q 13 ·∆ u = 0 from the transversality of Q 13 . Very similar arguments show that all contributions to B 1 vanish identically 7 in agreement with the Breit-frame conclusion (3.22).
We see however no reason for the same trivial result to hold for N ≥ 5. Also, as we shall see in the following section, a non-vanishing result will emerge at the next order in ω. 4 The O(ω 2 ) correction for a generic massless process At tree level, the sub-sub-leading terms come from the sum of three contributions: where s stands for the graviton's polarization and As in the previous section we have introduced arrows to indicate which amplitude (S 0 or S * 0 ) the operators J i , J j act on. However, unlike in the sub-leading case, we are now facing the problem of how to order the double derivatives appearing in S 2 . Hereafter we will use the prescription 8 that all derivatives act on S 0 or S * 0 before any possible multiplying factor. We will come back to this possible ambiguity when discussing specific processes. This being said, in the following we shall usually omit to write S 0 , S * 0 in the formulae. Summing over s = ±2, one gets where: and Π µν,ρσ was already defined in (2.2). Thanks to gauge invariance, momentum conservation and Lorentz invariance, we expect that all terms depending onq vanish after summing over i and j. This is easily seen if one observes thatq always appears in combination with q in Π. Since qJ i q = 0 = qJ j q, q will have to contract either with p i or with p j (or with both). In any case the factor qp i (or qp j ) will cancel the pole and the sum over i (or j) will vanish thanks to momentum or 7 At least if we neglect log corrections that appear to sub-leading order beyond tree-level. On the other hand, recent studies [40,65] that keep those effects into account appear to confirm this conclusion. 8 We acknowledge useful correspondence with Matin Mojaza and Paolo Di Vecchia about this important issue.
where W ij µν is given in (4.5) and The contraction W ij µν q µ q ν gives the following terms (from B 02 + B 20 and B 11 , respectively): . (4.10) As it was already the case for B 1 , the individual integrals for fixed i and j are affected by 'collinear' divergences, when n = q/ω is parallel to either v i or v j . Although these divergences cancel after summing over i and j, it is convenient to shift the integration variable q in such a way that integrals be finite for fixed i and j. The obvious choice for the shift is the one already used for B 1 .
That shift is ill defined for i = j, but there were no such contributions in B 1 . Terms with i = j are present in B 2 and have to be treated separately. Instead, for i = j, we apply the replacement In the following, to simplify notation, we will often suppress the indices ij inq ij . Notice that while q 2 = 0,q 2 = −2qp i qp j /p i p j Also at variance with what happened for B 1 , shifting q does not leave individual integrals appearing in B 2 unchanged. The explicit calculation is reported in appendix A. The final result of such calculation is:

JHEP05(2019)050
where D i = p i ∂ i (no sum). We now proceed as in the case of B 1 by introducing a Lorentz invariant constraint such as δ((qP/Λ 2 ) + 1), with P taken to be the total incoming momentum, P = p 1 + p 2 , and ω 0 = Λ 2 / √ s (where s = −P 2 = 4E 2 = E 2 CM ) the center-ofmass frequency at which we wish to compute the gravitational wave spectrum.
The last integral can be easily computed in the CM frame where P = (2E, 0) and q = −ω(1, n) (emitted radiation, η = −1) In order to evaluate the remaining 'shifted' integrals: one observes that M µν ij = M νµ ij (symmetric) and p i,µ M µν ij = 0 = p j,µ M µν ij (bi-transverse). As a consequence M µν ij can be written in the form After 'tracing' with η µν and contracting with P µ P ν one gets where we encounter again: (4.16) The first 'scalar' integral η µν M µν ij is easily computed while the second 'scalar' integral P µ P ν M µν ij requires more work and gives the result (see appendix B, eq. (B.10)):

Contracting with (B)H and (A)P yields
Including the coefficient functions A and B, defined in (4.19) and (4.20), and inserting the result in eqs. (4.12), (4.14), we arrive at our final expression for B 2 (after renaming again ω 0 as ω): where we recall that Q µ ij ≡ P µ − P p j p i p j p µ i − P p i p i p j p µ j and that, by convention, all derivatives act on the amplitude before any possible multiplier.

General covariant result for 2 → 2 scattering
As in the case of B 1 a big simplification occurs for a four-point amplitude. We will continue using the recipe for the derivatives that led to a vanishing result for B 1 .
In the case of C 1 things are very simple. Using eq. (3.25) we get:

JHEP05(2019)050
Let us now consider C ij 2 : it can be re-expressed as follows: where: Let us consider, as an example, the case of L 12 2 + L 21 2 and use our rules (3.24) for taking derivatives. Then: We have also used the orthogonality between H 12 and p 1 , p 2 to get rid of the derivatives w.r.t. s. The same result is obtained for L 34 2 + L 43 2 and, mutatis mutandis, for the other pairs. Summing up the different pairs, each with its own K ij 2 factor gives: Consider finally C 3 . In this case there is no contribution from i, j = 1, 2 and i, j = 3, 4 since, as already noticed for B 1 , Q 12 = Q 34 = 0. Like for C ij 2 we can decompose C ij 3 as: where: (5.10) Consider, for instance, L 13 3 . Using the rules (3.24) and taking into account the orthogonality between Π 13 and p 1 , p 3 , we easily find:

JHEP05(2019)050
and Adding up finally all the contributions we obtain: where we recall that: and Note that the operators D and ← → ∆ st , . . . are unambiguous in the sense that they give the same result when acting on A(s, t), The previous expression can be further simplified by using the easily proven identity: allowing us to re-express B 2 as follows: A few comments on (5.17) are in order: • The result is symmetric, as it should, between t and u. It is not, instead, with respect to s since we are working in the s-channel center of mass; • The (1, 1) contribution to (5.17), being an absolute square, should be positive. This can be checked to be the case by isolating the terms that contain one derivative acting on S and one acting on S † . Since the combination of operators appearing in ← → ∆ 2 st and ← → ∆ 2 su are negative-definite, their pre-factors must be negative as well. It is straightforward to prove that this is indeed the case in the s-channel physical region.

JHEP05(2019)050
• A non trivial check of (5.17) consists in considering the t → 0 (or equivalently, given its symmetry, the u → 0) limit. If we go back to eq. (1.6) for a 2 → 2 process, and consider the forward limit p 2 + p 3 , p 1 + p 4 → 0, we can easily check that both S 0 and S 2 (at least naively) vanish in that limit while S 1 does not. This is because, in the limit, the contributions of i = 2 and i = 3 have opposite denominators and equal (opposite) numerators for S 0 and S 2 (S 1 ). The same happens of course for the i = 1 and i = 4 contributions. This is confirmed by the vanishing of the leading term S * 0 S 0 in the forward (t = 0) or backward (u = 0) direction, see eq. (2.7). Naively, one should expect the same to be the case for the S * 0 S 2 interference term. Such terms are easily identified in (5.17) and give: where one should stress that the ordering of the derivatives is important for the cancellation.
• Extra terms can originate from a careful evaluation of the leading contribution at sub-leading level. An example is the constant −4 appearing in (2.9). Also, as already mentioned, the spectrum dE GW /dω depends on which variables are held fixed while the graviton's angular variables are integrated over. Two of them are obviously the total center-of-mass energy and the graviton's frequency (or, equivalently, (p 1 p 2 ) and (p 3 p 4 )). The angle-integrated spectrum, however, must depend on a third variable somewhat related to the scattering angle of the underlying 4-point function. One symmetric choice would be to keep (p 1 p 3 ) + (p 2 p 4 ) − (p 1 p 4 ) − (p 2 p 3 ) fixed, but it's by no means unique. In the next section we will see an explicit example in which these extra sub-leading terms are crucial for carrying out an important check.
We end up this section by applying the general result (5.17) to two examples of gravitational bremsstrahlung from two-particle collisions. As we will see in the first example the vanishing of the S * 0 S 2 term at t = 0 does not occur. This is because contributions of the type t 2 − → ∂ 2 t are very sensitive to the amplitude they are acting on. They do not vanish, for instance, if the amplitude has a pole at t = 0 as in the example presented below. In fact, already the vanishing of S 2 in eq. (1.6) for p 2 + p 3 → 0 crucially depends on which amplitude the J i operators act on. In any case there are (tree level) amplitudes for which the above check must hold and they represent a very good test of our final result.

Tree-level gravitational scattering
Let us first consider the tree-level gravitational elastic scattering of two different scalars. In this case, the amplitude receives only a t channel contribution giving (up to an irrelevant numerical factor)

JHEP05(2019)050
We will then apply eq. (5.17) to eq. (5.19) by defining: where, as already explained, we define the derivatives appearing in B 2 to act on M or M † first, i.e. before any possible multiplier. We also recall that the differential operators appearing in (5.17) can be taken to act on any one of the expressions in (5.19) according to convenience. We first note that the operators ← − D 2 and − → D 2 appearing in (5.17) can be replaced, after the above mentioned ordering of the derivatives, by , respectively. But since D gives just the s, t, u dimensionality of M, which is 1, these operators give simply a vanishing result.
Let us now turn to the other two differential operators appearing in eq. (5.17). For the operator ← → ∆ 2 su it is most convenient to use the expression for M that does not contain u so that we only need to consider the derivatives w.r.t. s. These are readily computed to give, after some simple algebra: Similarly, for the operator ← → ∆ 2 st it is more convenient to work with the s-independent form of A. This gives: Inserting now (5.21) and (5.22) in (5.17) we find the rather elegant result: This can then be converted into the following result for the sub-sub-leading correction to the (unpolarized) flux: (5.25) and, in spite of appearance, there are no singularities at either t = 0 or u = 0. Re-expressing the Mandelstam variables in terms of the cosine of the scattering angle x ≡ cos θ s , we can rewrite eq. (5.25) as The function f (s, t, u) is displayed in figure 1. It is positive, with a maximum of 3/2 at t = 0 (reached with an infinite slope) a value of 1/2 at u = 0 and a minimum of about 0.25 at x ∼ −0.3. We stress however again that only the sum of the leading and non-leading contributions have an non-ambiguous meaning. If, for instance, the leading term is defined through eq. (2.3), the O(ω 2 ) correction is given by because of the −4 appearing in eq. (2.9). As a very non trivial check of our procedure (in particular of our recipe for defining partial derivatives) we present, in appendix C, an ab-initio tree level calculation of the single-graviton emission amplitude from two-body massless-scalar scattering in N = 8 supergravity. Quite remarkably, the soft expansion of the result confirms both the vanishing of the O(ω) correction to the unpolarized flux and its precise functional form (5.25) or better (5.28) at O(ω 2 ). The full agreement requires the already mentioned careful identification of sub-leading terms implicitly contained in the leading terms of the two calculations. They are the same up to the constant −4 appearing in eq. (2.9). The latter is thus an important ingredient for the success of the check.
To summarize, we have found that in this case the leading correction to the ZFL of (2.6) is of relative order ( ω) 2 /Q 2 where Q ∼ √ −t is the momentum transfer in the process. The correction looks like a quantum effect if we define the classical limit as → 0 at fixed √ s and ω (which is the usual way one deals with gravitational wave fluxes in General Relativity). However, if we use the uncertainty principle to replace Q by /b, it becomes a classical-looking correction of relative order (ωb) 2 which is in agreement with

JHEP05(2019)050
expectations about soft-graviton theorems. The positive sign may look surprising since, by energy conservation, the constant ZFL value has to leave the way to a decreasing spectrum at higher frequency; however, a maximum at ω ∼ 1/b has been found to occur [65] in the process discussed in the next subsection.

Gravitational scattering in the leading eikonal approximation
As a second example we consider the same gravitational scattering in the very high (i.e. transplanckian) energy regime where an all-loop resummation is needed and has been carried out in [44] (see also [42,43]) in the small deflection angle limit. In that case the elastic S-matrix takes a very simple form in impact-parameter space (we recall that the impact parameter b is related to the orbital angular momentum by b = 2J/ √ s and is related by a Fourier transform to the momentum transfer Q): where L is some infrared cut off screening the (unobservable) infinite Coulomb phase already well-known in electromagnetic scattering. Note that this expression has the form of a semiclassical approximation (valid at Gs 1) with Gs playing the role of a classical action. Indeed (5.29) leads to known classical features of gravitational scattering: • The derivative of the exponent w.r.t. b provides the correct gravitational deflection angle (the generalization of Einstein's deflection angle to the case of massless-particle collisions).
• The derivative of the exponent w.r.t. E = √ s gives the (Shapiro) time delay as a function of the impact parameter. Note that while the deflection angle is independent of L, the time delay is not. Fortunately what matters are normally time-delay differences for which L again drops out.
Let us now apply B 1 and B 2 to (5.29) in order to compute the low-energy gravitational radiation accompanying the (otherwise elastic) collision. The first qualitative remark to be made is that, to leading order in , we can keep only the action of the differential operators on the exponent (just like in the usual WKB approximation). As a result each derivative brings down a −1 factor that precisely compensate for the explicit positive powers of appearing in dE GW /dω. In other words we obtain spectra which have a smooth classical limit unlike in the case discussed previously.
At a more quantitative level we find that, also in this case, B 1 of eq. (3.15) gives a vanishing result, when acting on (5.29). The "surprise" comes from B 2 . If we use our final expression (5.17) for B 2 and we let the derivatives act just on the exponent we find: √ s. Thus we recover, as expected, a classical correction to the ZFL.

JHEP05(2019)050
Unfortunately, (5.30) is IR divergent. The divergence comes from the infinite Coulomb phase which is typical of four-dimensional physics. Such an infinity should disappear from any physical observable meaning that the initial recipe for computing the sub-leading terms has to be modified. In recent papers [39,40] the authors have proposed to replace the infrared cutoff L by ω −1 . This recipe has been supported by known classical results [66] and also, at O(ω) level, by the independent IR-singularity free method developed in [65].
If we apply such a recipe to (5.30) we get the O(ω 2 ) correction to the (unpolarized) energy flux in the form: This result agrees with the one obtained in [65] up to the overall (positive) constant that cannot be fixed in the small-angle approximation used in [65]. As pointed out there, since (5.31) represents a positive leading correction to the ZFL, it implies a maximum of dE GW 2 /dω displaced from ω = 0 by an O(1/b) amount. It also provides another independent check of the validity of the recipe proposed in [39,40]. In the next section we will offer some ideas on how one can try to justify that recipe through a suitable modification of the soft theorems themselves.

Elimination of (and finite terms from) infrared divergences
In order to understand the origin of our (hopefully spurious) IR singularity in (5.30) it is better to step back and consider the sub-leading correction S 1 before computing the cross section. Let us also restrict ourselves to n = 4 If we use for M 4 (p i ) the (Fourier transform of) the eikonal expression (5.29) and let the derivatives present in J i act on the exponent we find: where the restriction of the sum over the i, j pairs to be both incoming (η i = η j = +1) or outgoing (η i = η j = −1) comes from the way the Coulomb phase originates (see last section of [2]). We note that essentially the same term came out from the approach of [25,39] where it was argued that the IR cutoff should be replaced by a log(ω −1 ). We see clearly here the origin of the Coulomb divergence. This propagates to the S * 0 S 1 contribution to B 1 but cancels with the S * 1 S 0 contribution thanks to its over all imaginary phase. We will now argue that such a divergence is spurious and should be cancelled by contributions not included in the naive definition of S 1 .
Consider indeed two classes of sub-leading contributions according to whether the final soft graviton is emitted from an outgoing or an incoming leg. In the former case the initial JHEP05(2019)050 state participating in the 2 → 2 sub-process gives the usual contribution to the Coulomb phase, but the two hard final particles fail to do so since their total momentum is the initial one minus the soft graviton's. The mismatch would produce exactly terms like in (6.2) with η i = η j = −1. But, obviously, the mismatch is compensated by the fact that also the final soft graviton contributes to the total Coulomb phase by its own rescattering on the hard final particles. For the latter kind of contributions (emission from initial particles) the situation is similar for the final particle's contribution to the phase. For the initial particles contribution one can argue that the Coulomb-divergence comes from the asymptotic (initial and final) states and therefore the initial state to be considered in not the one entering the blob but the actual initial state carrying the whole energy.
In view of the above it looks that the naive recipe for S 1 is wrong on two accounts. The contribution to (6.2) with η i = η j = +1 should be simply omitted while the one with η i = η j = −1 should be supplement with a new term due to the rescattering of the soft graviton on the other final particles. We note that a similar contribution also appears in the treatment of [60] and is crucial in order to recover agreement with [59].
We may ask whether other modifications of the naive recipe for S 1 are needed because of the standard IR divergences. It was argued in [22] that the correct replacement in that case is: where 1/ plays the role of our log L. This suggests that the true IR divergences are cured as usual by the interplay between real and virtual contributions. Furthermore, we will argue that, provided ∆E ω there should be no important finite leftover contribution to dE GW /dω.
The argument goes as follows. Consider a generic process in which there is a certain number of hard external particles of characteristic energy E and a number of soft external particles of energy around ω E, the energy/frequency at which we wish to compute the spectrum of GWs. The virtual corrections can be split into those with momenta larger or smaller than ω. The former contribute only when the gravitons are exchanged between the hard particles, the latter contribute to all graviton exchanges. Consider now the real emission corrections to that process distinguishing those with emitted gravitons of energy less than ω and those with higher energy up to the ∆E upper bound. The latter will cancel the IR divergences of the corresponding virtual contributions leaving behind the usual logarithmic dependence on E/∆E. The crucial point is that, instead, the real gravitons softer than ω exactly cancel their virtual counterpart (since for both the real and the virtual gravitons the upper limit is ω). Therefore, provided ∆E ω, the net result is nothing but the usual infrared factor for the no-emission amplitude which cancels out in dE GW /dω. Obviously, some dependence on ∆E/E will remain as long as that quantity is finite.

JHEP05(2019)050
In conclusion, while we expect important corrections to the naive recipe to come from the IR divergence related to the Coulomb phase (which has no real counterpart), the usual cancellation mechanism between real and virtual IR divergences should leave room to no substantial correction as long as ω ∆E E. A complete and general calculation of the finite corrections coming from the Coulomb phase goes beyond the scope of this paper.

Summary and outlook
In this paper we have combined the low-energy theorems for graviton emission, together with some reasonable assumption about how multi soft graviton emission exponentiates to give a coherent state, in order to compute the spectrum of gravitational radiation emitted to leading, sub-leading, and sub-sub-leading order in the frequency. For simplicity we have only considered the case of massless (equivalently highly relativistic) spin-less particles, and we have summed over the gravitational wave (GW) polarizations and integrated over the angular distribution of the radiation. The results have been expressed therefore in terms of the GW energy flux dE GW /dω.
At leading order we could easily reproduce the known [1] constant zero-frequency limit (ZFL). At sub-leading order we obtained a general result which happens to vanish for GW emission from a two-body collision. For the case of the unpolarized and/or angle-integrated flux this is in agreement with other results [25,65].
At sub-sub-leading order (i.e. to O(ω 2 )) we derived a nice and fairly compact expression for the general case which further simplifies in the case of a two-particle collision process. The expression passed a number of non trivial checks including one against an ab-initio calculation in N = 8 supergravity.
An important conclusion of our method is that the soft theorems, while presumably valid at tree-level or in D > 4, must be amended at higher orders in D = 4 because of IR problems. Actually there are two kinds of IR problems that can invalidate, in principle, the soft-graviton theorems.
The first is the usual IR catastrophe which is solved by taking into account both real and virtual IR divergences and their cancellation for physically measurable, finite resolution, observables. We have argued, in section 6, that these divergences should be also harmless as long as dE GW /dω is concerned and the energy-resolution is parametrically larger than ω.
The second kind of divergences comes from the familiar infinite Coulomb phase, also known in QED. Usually, such a phase is unobservable and cancels out when computing physical quantities. However, when the differential (angular-momentum) operators, appearing in the soft theorems at non-leading order, act of such a divergent phase they generate infinities that do not cancel out in dE GW /dω. 9 This IR problem has already been noticed in recent papers [23][24][25], whose authors came up with a recipe for dealing with the problem at sub-leading order. Such a recipe has been confirmed by checks against known results and against calculations based on the eikonal method [65] which are free of such problems.

JHEP05(2019)050
Using the recipe to the sub-sub-leading order we find that the leading correction to the ZFL has exactly the structure found independently in [65] which implies a bump in the GW spectrum at ωb ∼ 1.
For the future it would be clearly important to develop a more rigorous approach to dealing with IR problems in D = 4 and, by reversing the Cachazo-Strominger argument [27], to understand the outcome in terms of an anomaly in the extended BMS symmetry responsible for the soft-graviton theorems at sub-leading level (see [67] for a recent review). The other direction of research would be to extend the formalism to massive and/or spinning particles in the initial and final states. Such calculations would nicely complement those presently under way (see e.g. [68]) for computing the conservative part of the effective potential at the third post-Minkowskian (3PM) level by providing its dissipative (radiation) counterpart.

JHEP05(2019)050
which, by rewriting the i =j as a sum all over i, j and subtracting to it the diagonal terms, can be rewritten as The (−3/2) factor in eq. (A.1) comes from three −1/2 factor for each of the three initial terms in (11); the factor (+5/2) is the sum of a coefficient 1 from the first, a coefficient 0 from the second and a coefficient +3/2 from the third.
• For B 20 and B 02 one has to be careful and take into account that, for i = j, We find: where (↔) means reversing the arrow-orientation. The last term of eq. (A.3) is the sum of two shift terms: from i = j (giving a numerical factor 1/2) and from i = j (providing a 3/2 factor).
Finally, combining (02) + (20) + (11), we find: where the last term is obtained by combining the last term of eq. (A.2) with the last term of eq. (A.3). This is the result we inserted in (4.12) modulo the fact that, according to our ordering of the derivatives, the second term in (A.4) has to be ordered accordingly.
B Non covariant calculation of P µ P νM µν Relying on Lorentz invariance, one can compute the relevant integral in the CM frame whereby

JHEP05(2019)050
The angular integral can be rewritten as The first two integrals were already computed, for B 0,1 in the 'covariant' approach in section 3. In the mass-less limit under consideration, after introducing a small mass regulator, one finds where (in the mass-less limit) Thanks to the shift, as noticed in section 3, the Lorentz invariant final result is free of divergences in the massless limit m i → 0. This requires cancellations between L i , L j and the log terms in M 0 ij . The M ab ij integral matrix can be decomposed as In order to determine the 'scalar' integrals α, β, γ, one can project M ab ij along the three independent components:

JHEP05(2019)050
where, choosing reference frame so that v j = |v j |(0, 0, 1), v i = |v i |(sin θ i,j , 0, cos θ i,j ) and n = (sin θ n,j cos φ n,j , sin θ n,j sin φ n,j , cos θ n,j ), In order to compute this integral matrix, we need to calculate After some long and tedious algebra one gets Combining with the other two terms one finally gets As a check notice that the result vanishes, as expected, , as expected, one gets a constant 16π. Using the explicit expressions for M 0 ij , L i , L j one can check Lorentz invariance and finiteness (as m i → 0) for generic values of v i v j . Indeed the result can be written as In this appendix, we study the soft limit of the 4-scalar + 1 graviton amplitude in order to compare the exact results with the ones obtain from soft theorems. We will show the perfect agreement between the two approaches provided one takes carefully into account all sub-leading terms, including those somehow hidden below the leading one. In addition to the fermions, the N = 8 SUGRA multiplet contains the graviton h µν , which is a singlet of the SU(8) R-symmetry, 28 gravi-photons A
The expression for V can be further simplified to By the same token, one finds U qp 1 qp 2 qp 3 qp 4 = −4 C 1 + s 2 t 2 C s/t + u 2 t 2 C u/t + s 2 u 2 C s/u + u 2 s 2 C u/s + t 2 s 2 C t/s + t 2 u 2 C t/u , (C. 29) where (recall Q i = qp i /2) We would like to determine the soft graviton spectrum at this order (ω 2 ). To this end one should integrate expressions of the form Integrating over the graviton phase space, from the U term, one gets:

JHEP05(2019)050
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.