Graviton particle statistics and coherent states from classical scattering amplitudes

In the two-body scattering problem in general relativity, we study the final graviton particle distribution using a perturbative approach. We compute the mean, the variance and the factorial moments of the distribution from the expectation value of the graviton number operator in the KMOC formalism. For minimally coupled scalar particles, the leading deviation from the Poissonian distribution is given by the unitarity cut involving the six-point tree amplitude with the emission of two gravitons. We compute this amplitude in two independent ways. First, we use an extension of the Cheung-Remmen parametrization that includes minimally coupled scalars. We then repeat the calculation using on-shell BCFW-like techniques, finding complete agreement. In the classical limit, this amplitude gives a purely quantum contribution, proving that we can describe the final semiclassical radiation state as a coherent state at least up to order $\mathcal{O}(G^4)$ for classical radiative observables. Finally, we give general arguments about why we expect this to hold also at higher orders in perturbation theory.

1 Introduction 1 2 Graviton particle statistics from on-shell amplitudes 4 2.1 Graviton emission probabilities in the KMOC formalism 4 2.2 Mean, variance and factorial moments of the graviton particle distribution 7 3 Tree amplitudes from Feynman diagrams

Introduction
The two-body problem in general relativity has been receiving increased attention since the first detection of gravitational waves. Quantum field theory methods have recently proven to be very useful to understand the perturbative long-distance regime of the twobody dynamics, offering a new perspective for understanding the inspiral phase in the post-Minkowskian (PM) expansion in the spirit of effective field theories (EFT) [1,2].
On-shell scattering amplitudes techniques, powered by locality, unitarity and double copy, have been used to get compact analytic expressions for the state-of-the-art binary dynamics for spinless pointlike bodies at 3 PM order and partially at 4 PM order [3][4][5][6][7]. A handful of alternative and complementary approaches have also been developed in recent years. The relativistic eikonal expansion [8][9][10][11][12][13][14] and semiclassical worldline tools [2,[15][16][17][18] offer many insights on the binary problem, both at the conceptual and at the practical computational level. Moreover, the formalism can be extended to include both spinning bodies [19,20] and finite size effects [21,22] in terms of additional higher-dimensional operators. All these approaches share the need for careful analysis of the classical limit, as done in seminal work by Kosower, Maybee and O'Connell (KMOC) [23]. In the conservative case, a dictionary has been found [24,25] which enables the analytic continuation of observables from hyperbolic-like scattering orbits to bound orbits, which ultimately are of direct relevance to LIGO.
The dynamics of the binary in the presence of radiation is much less understood compared to the conservative case. This is very important, for example to establish a direct connection with the waveforms [26][27][28][29][30]. Unitarity dictates that, even at the classical level, observables are IR-finite only when we include both real and virtual radiation, as stressed in [31,32]. This is crucial to obtain a well-behaved scattering angle at high energies [33], as was proven by a direct calculation of radiation reaction effects [11,12,[34][35][36]. A similar principle holds for less inclusive observables like gravitational energy event shapes [37].
We would like to understand the exact structure of the final semiclassical state, including classical radiation. Many recent insights, coming both from a pure worldline description [16,[38][39][40] and from a different parametrization of the kinematics in the classical limit [34,35,41,42], suggest that we should expect an (eikonal) exponentiation at all orders in the impact parameter space. The situation is less clear when we allow particle production. Since we expect the description of a classical wave for a pure state to be possible in terms of a single coherent state [27,43], a naive crossing-symmetry argument suggests that it should also be possible to describe the final radiation in terms of coherent states. A lot of attention has been devoted so far to the soft expansion where coherent states arise naturally from classical currents [44][45][46][47][48][49], but the dynamics of how these states are generated by the scattering process is much less clear. 1 In this work, we compute the expectation value of the graviton number operator using the KMOC formalism, and we show how this is connected to unitarity cuts involving amplitudes with gravitons. Similar ideas in a purely off-shell Schwinger-Keldysh formalism have been developed in [51,52]. Since coherent states correspond to Poissonian distributions at the level of the particle emission, deviations from such structure imply that we cannot represent the final graviton state as a coherent state. We will show that the leading contribution is given by a unitary cut involving the 6-point tree amplitude A (0) A complementary perspective is provided by the study of the factorization of radiative observables in the classical limit, as discussed in [53].
There are various approaches which can be taken to compute the amplitude in ques-tion. Applying a traditional Feynman-diagram method to gravitational theories is notoriously difficult, due to the multiplicity of gauge-dependent vertex rules contributing to the amplitudes. Inspired by the simplicity of on-shell amplitudes methods, Cheung and Remmen [54] rewrote the Einstein-Hilbert action in a simpler form through the introduction of an auxiliary field, in the spirit of the first-order formalism developed by Deser [55]. The result is a set of Feynman rules for pure gravity which can be used to compute graviton amplitudes in a very efficient way [56]. Here, we will extend this construction to matter by adding minimally coupled scalar fields. 2 An alternative approach is to use on-shell BCFW recursion [58], which has a number of benefits over off-shell approaches coming from Lagrangians. The first is that the only objects needed are exceptionally simple seed amplitudes which form the base point of the recursion. This approach can in principle be used to reproduce all tree amplitudes of a variety of massless QFTs, including Einstein gravity [59,60]. While the introduction of massive matter does not necessarily obstruct the BCFW construction of higher-point amplitudes [61], the method generally relies on having massless external particles whose null momenta are used to construct a linear momentum shift, and further generally requires the absence (or good behavior) of boundary terms in a corresponding complex integral. Recent works have explored the prospect of applying shifts to massive legs, as well as combining the shift with a soft recursion relation to construct higher-point amplitudes in general massive theories [62][63][64][65]. Here, we introduce a new shift capable of reproducing the hard collinear factorizations in a mixed gravity-scalar theory.
Using this "equal-mass shift" we will show that it is possible to obtain a compact form of A (0) , whilst maintaining gauge invariance at every stage in the computation. We then use the standard BCFW shift to compute the leading classical scaling of A (0) , which is one of the main results of the paper.
A summary of the paper is as follows. In section 2, we study the graviton number operator expectation value in the KMOC formalism for the two-body problem, which is given in terms of unitarity cuts involving on-shell amplitudes with graviton emissions. The deviation from the Poissonian statistics is equivalent to a deviation from coherence in the final semiclassical state with radiation. In section 3, we lay out the Feynman-diagram computation, establishing the Feynman rules in the extended Cheung-Remmen formalism with minimally coupled massive scalar particles. In section 4, we repeat the computation of the five-and six-point amplitudes using an on-shell equal-mass shift and BCFW recursion. In section 5, we study the scaling of these amplitudes in the limit → 0, and we prove that the six-point tree amplitude does not contribute to the total energy emitted with classical gravitational waves. Moreover, assuming coherence, we also establish new relations in the classical limit between unitarity cuts of amplitudes involving an emission of gravitons in the final state. Section 6 contains our concluding comments. Appendix A discuss the derivation of the KMOC formalism based on the Schwinger-Keldysh approach, and appendix B summarizes the connection between Poissonian statistics and coherent states.
Conventions We work in the mostly plus signature, for consistency with the gravity action presented in [54].

Graviton particle statistics from on-shell amplitudes
In this section, we study the particle statistics distribution of the gravitons emitted in the scattering of a pair of massive point particles of mass m A and m B in general relativity, using methods of perturbative QFT. In particular, we relate the expectation value of the graviton number operator to a sum of unitarity cuts involving scattering amplitudes with external gravitons.

Graviton emission probabilities in the KMOC formalism
LetP n be the probability of emitting n gravitons in the scattering of a pair of massive particles as described above. Unitarity implies that ∞ n=0P n = 1. In quantum field theory, this statement is equivalent to a completeness relation in the Hilbert space, where |k σ 1 1 ...k σn n k σ 1 1 ...k σn n | is the n-graviton particle projector on states with definite momenta k 1 , ..., k n and helicities σ 1 , ..., σ n , whose values are indicated by the single + and − signs.
We denote the scattering matrix operator by S, the momenta of the incoming (resp. outgoing) massive scalar particles by p 1 , p 2 (resp. p 3 , p 4 ), and the outgoing graviton momenta by {k i } i=1,...,n . It is clear that the probabilityP n is given by taking the expectation value of the n-graviton particle projector, As it is written, (2.2) is formally divergent, as is known from the study of infrared divergences in quantum field theory (see [66]) because of the contribution of zero-energy gravitons. We will therefore work with a finite-resolution detector λ > 0, which implies that we will study only the probabilities of gravitons emitted with an energy E k > λ. Correspondingly, we will replacē As we will see later, we will not be interested in the single probability but in a particular infrared-safe combination of probabilities. Therefore λ will be used only as an intermediate regulator, and in the end we will send λ → 0.
We would like to scatter classically two massive point particles with classical momenta m A v A and m B v B with an impact parameter b µ . Since the main purpose of this paper is to take the classical limit from a quantum field theory calculation, we use the KMOC formalism [23] and take instead as our incoming state The wavefunctions ψ A (p 1 ), ψ B (p 2 ) are defined as where N 1 , N 2 are normalization factors, c,j = /m j is the Compton wavelength, and w,j is related to the intrinsic spread of the wavefunction for the j-th massive particle (j = A, B). We will also require the "Goldilocks conditions" which ensure that wavefunctions such as those in (2.6) effectively localize the massive particles on their classical trajectories as → 0. We expand the S-matrix in terms of the scattering matrix T , For the expectation value of the graviton projector operator, only the amplitudes with at least one graviton emitted are going to contribute. We can read off from (2.2) the probability of emitting n gravitons with energies E k i > λ, (2.9) We now expand |ψ in from (2.4) in terms of the wavefunction ψ A (p 1 )ψ B (p 2 ) and ψ in | in terms of ψ A (p 1 )ψ B (p 2 ), and we introduce the momentum transfers [23], Writing the momentum integrals explicitly and inserting the momentum constraints and amplitudes, we then arrive at We now conveniently define a set of symmetrized variables for the external momenta [35] which has the nice property of enforcing exactly the condition p A · q = p B · q = 0. In terms of these new variables [20], we have where we use the double bracket notation · introduced in [23], which contains the implicit phase space integral over p A , p B and the appropriate wavefunctions where Note that (2.11) is expressed in terms of unitarity cuts involving n gravitons and the two massive particles in the intermediate state. The same result can be obtained by applying the LSZ reduction with the appropriate KMOC wavefunctions from the in-in formalism, as shown in appendix A.
2.2 Mean, variance and factorial moments of the graviton particle distribution In classical physics, we are interested in knowing whether the final graviton particle distribution is exactly Poissonian or super-Poissonian (the most general case). We refer the reader to appendix B for a brief review of the two cases. Poissonian statistics are known to be equivalent to having a single coherent state representing the quantum state for the classical radiation field. Here we give a short argument [27] for why we expect a single coherent state, based on the fact that the we expect the incoming state to be a pure state in the classical limit and on the unitarity of the S-matrix. The work of Glauber in 1963 [67,68] shows that every quantum state of radiation (i.e. every density matrix) can be written as a superposition of coherent states, where P σ out (α k ) is a well-defined probability density (P σ out (α k ) ≥ 0) in the coherent state space in the classical limit, and |α σ k represents a coherent state of a graviton excitation ("harmonic oscillator") of momentum k and definite helicity σ, which we can write generically as where a † σ (k) and a σ (k) are the creation and annihilation operators of a graviton of helicity σ. This representation is known as the Glauber-Sudarshan P-representation [67,69], and it is widely used in the quantum optics literature. In quantum field theory, we need to consider an infinite superposition of harmonic oscillators for all momenta k ∈ R 1,3 , and therefore we will promote (2.16) to 3 where now Since we are dealing with scattering boundary conditions and our incoming KMOC state |ψ in is a pure state, the unitarity of the S-matrix SS † = 1 implies that |ψ in is mapped to outgoing pure states. Therefore, the outgoing radiation state must be a superposition of pure states, But thanks to a crucial theorem of Hillery [43], we know that every such superposition of pure states is trivial in the classical limit → 0, We therefore expect, on general grounds, to be able to describe the final radiation state for a scattering process involving point particles with a single coherent state. From the pure amplitude perspective, the same question is hard to answer unless we work strictly in the soft approximation [37,47,50]. But in general, we can address this question perturbatively by studying the mean, the variance and the factorial moments of the particle distribution. A similar approach has been taken by F. Gelis and R. Venugopalan [51,52,70] in the standard in-in formalism, which we try to specialize here from a fully on-shell perspective and in the classical limit.
The graviton number operator is defined aŝ Having defined the expectation value of the number operator in the final state gives the mean of the distribution, which can be expressed in terms of unitarity cuts in a similar fashion to the derivation of equation (2.11), as depicted in Figure 1. The mean of the distribution is defined as where |r 1 r 2 X denotes the state with n X gravitons and two massive particles of momenta r 1 and r 2 , and λ dΦ(X) stands for the phase space integration for the gravitons. We define the variance of the distribution as If the variance is equal to the mean, i.e. if then the distribution is consistent with a Poissonian distribution. This means that the deviation from the Poissonian distribution, characterizes the deviation from the coherent state description. We claim here that the difference between the mean µ λ and the variance Σ λ is an infrared-safe quantity in perturbative quantum gravity. While the probability of emission of n gravitons is generally ill-defined because of infrared divergences, there is a non-trivial cancellation which happens for ∆ out . Indeed, the contribution of zero-energy gravitons to the final state, which give rise to the infrared divergent contributions, is known to be exactly represented by a coherent state. This can be proved either from a Faddeev-Kulish approach [37,50] or from a path integral perspective [38,40,53]. Let us denote the mean and the variance of this coherent state for zero-energy gravitons by µ E k ∼0 out and Σ E k ∼0 out respectively. In Appendix B, we show that for such a coherent state of soft gravitons we have 4 This is the reason why the cutoff λ was removed in (2.27). 5 We can easily check by induction on the number of loops and legs that 6 where we have explicitly extracted the scaling in the gravitational coupling G of the product of an L 1 -loop amplitude with an L 2 -loop amplitude with n gravitons. The lowest order contribution to ∆ out is of order O(G 4 ), which corresponds to This leading term is the unitarity cut involving the 6-pt tree amplitude A (0) . It is important also to understand the higher order terms in ∆ out , since they will give non-trivial amplitude relations if we assume coherence at all orders. From the definition (2.27), we have Let us examine the first several terms appearing explicitly in the expansion of (2.31), It is not necessary to specify α σ E k ∼0 (k) for the argument to work. The interested reader can find additional details in [50]. 5 This argument does not apply directly to non-abelian theories because of the presence of collinear divergences, which for perturbative gravity are known to cancel exactly [71]. It would be interesting to develop this idea further, along the lines of [72,73]. 6 To avoid cluttering the notation, we keep the λ dependence implicit in P where we have organized each different line according to the expected behavior of the terms in the classical limit. We expect that the first three lines of (2.32) are related to "quantum" contributions and are therefore irrelevant in the classical limit. The last line of (2.32), instead, contains a combination of unitarity cuts which will give non-trivial quadratic relations between "classical" loop amplitudes with a higher number of emitted gravitons of the form P (L 1 ,L 2 ) n with n ≥ 2 and L 1 + L 2 ≥ 1, and 5-point amplitude contributions involving P (L 1 ,L 2 ) 1 . We will discuss this interpretation in more detail in section 5, where we will also emphasize the relevance of the 5-pt amplitude for the calculation of classical radiative observables.
It is important to consider also higher moments of the statistical distribution for the graviton number production. We can define a generating functional from which all higher moments can be derived, Therefore, the knowledge of all graviton emission probabilities P λ n is enough to completely determine the distribution of the particles above the energy cutoff. In practice, we can rely on perturbation theory and therefore computing the first few moments is enough to accurately determine the particle distribution. We can also defined connected moments (or "cumulants"), like the variance and its higher order generalizations. Having defined a generating functional for a Poissonian distribution we would expect, given a certain waveshape α σ (k), that because all the cumulants should be equal. In particular, the variance is a special case for m = 2, i.e. Σ (2),λ = Σ λ . For our purposes it is more convenient to consider factorial moments Γ (m) , which correspond to a linear combination of the connected moments discussed above. We define the factorial moments For a Poissonian distribution it is possible to prove that and therefore we can also consider in perturbation theory other infrared-safe combinations of probabilities like where for m = 2 one can check that we recover the difference between the mean and the variance in (2.31). By expanding (2.39) we get immediately It is interesting to consider the first terms in this expansion of ∆ where we have organized the terms similarly to what was done in (2.32). We will explore the deep consequences of assuming coherence at all orders, i.e. ∆ (m) out = 0, in section 5. In [53], it is shown how coherence properties are linked to the factorization of radiative observables in the KMOC formalism. 7 In classical physics, we expect only the 1-point function to play a role for any observable of interest. Such an observable is essentially uniquely determined by the classical equations of motion and the retarded boundary conditions at t → −∞: all two-point and higher-point functions then have to factorize as → 0. There the following relation was established,

⇐⇒
Zero-variance property in the Glauber-Sudarshan coherent state basis 7 Similar statements about the classical factorization have been made in [37] for infrared divergences and in [27] for the classical expansion.
which implies that Poissonian distributions in the number operator basis correspond to a degenerate distribution (∝ δ 2 (α σ − α σ )) in the Glauber-Sudarshan space.

Tree amplitudes from Feynman diagrams
In this section, we extend the parametrization of the pure lagrangian used by Cheung and Remmen [54] to the case of real scalar fields minimally coupled with gravity. This will make use of an auxiliary field, the connection, whose job is to effectively resum higher order graviton pure contact vertices in the same spirit as the first order Palatini formulation developed by Deser [55,74]. We can then compute in a straightforward way all the tree level amplitudes we need for this work.
Let us consider the lagrangian of two real scalars minimally coupled with gravity in D = 4 dimensions, where we work in the mostly-plus signature and have used the following conventions: , We introduce the auxiliary field A a bc , which allows us to rewrite the pure gravity lagrangian as Before setting up the perturbation theory in the new variables, it is useful to unmix the graviton and auxiliary field by doing the shift and adding the gauge fixing term Using we get explicitly up to O(h 4 ) a lagrangian of the form The quadratic terms in the lagrangian are given by 9) and the interaction terms are 8 (3.10) In the massless limit m A , m B → 0, the interaction terms become purely trivalent. In that case, it is possible to set up the standard Berends-Giele recursion relations. But even with the mass terms, the final expressions are more compact than in the standard perturbative expansion of gravity: the gravity pure self-interactions are nicely resummed by the auxiliary field, which makes it possible to avoid the cumbersome expressions for higher point vertices (at least at tree level, where ghosts are absent). The Feynman rules for the propagators are then 11) and the rules for the interaction vertices are where all momenta are chosen to be ingoing. At this point one can implement these Feynman rules in the xAct package [75], which we use extensively in the following calculations.
For the purposes of simplifying computations, we adopt the following conventions for the momenta of our amplitude: and we define the momentum invariants s ij = −(p i + p j ) 2 , with Mandelstam invariants defined as s = s 12 and t = s 13 in the particular case of four-point kinematics.

Four-point and five-point tree amplitude
We have only one diagram in the 4-pt case, given in Fig. 2. The Feynman rules give the well-known result 9 (3.14) For the 5-pt amplitude, we have explicitly computed the 7 diagrams pictured in Fig. 3. Notice that the first 6 diagrams are in one-to-one correspondence with the analogous calculation in scalar QED [77], while the last one is related to the graviton self-interaction.   9 See for example eq. (3.1) of [76], with D = 4 and κ = √ 2κ4.

Six-point tree amplitude
We have computed the 68 diagrams in Fig. 4 for the 6-point tree amplitude. In order from the top left of the picture in Fig. 4, the first 42 of these diagrams can be compared with the analogous calculation in scalar QED done in [53], which in particular involve the 3-point and the 4-point vertices with one matter line and one or two gravitons. The remaining 26 diagrams are classified into the following three types: • 21 diagrams involving the graviton self-interaction; • 3 diagrams with the auxiliary field propagator; • 2 diagrams with a 5-point contact vertex with 3 gravitons and one matter line.
We have highlighted in red the contribution of the auxiliary field, which is crucial to obtain the correct result. The calculation of these tree level amplitudes agrees exactly with an independent onshell BCFW calculation presented in the next section.

Tree amplitudes from on-shell recursion relations
In this section, we compute the necessary tree-level amplitudes for the theory defined in equation (3.1) by using an on-shell diagrammar 10 to recursively construct all the amplitudes in the theory. A diagrammar requires basic amplitudes to serve as the atoms of the computation, and the on-shell recursive framework of BCFW [58]. In massless theories there are straightforward arguments to construct three-point amplitudes from little-group scaling [59,60]. The simplicity comes from the on-shellness of the momenta, which is maintained throughout the computation and simplifies the expressions needed as input.
We begin with a brief review of BCFW recursion, in preparation for the new shift that we will introduce to compute the 5-point tree amplitude A and to set the stage for its application to the 6-point tree amplitude A (0)

Review of BCFW
The basic mechanism of BCFW recursion is understood through elementary complex analysis. The derivation begins by introducing a complex variable z and considering a linear shift in (a subset of) the momenta p i in the (yet-to-be-determined) n-point tree-level amplitude: where the shifted momenta are defined aŝ The choice of r i corresponds to a choice of shift.
As tree amplitudes are rational functions, we can consider A (0) n ({p i }) as a meromorphic function of z which we denote as A (0) n (z). We then evaluate the contour integral where the z I are the poles in the complex plane, and the integration contour γ ∞ := lim R→∞ γ R , where γ R is a circular contour around the origin with radius R. The choice of the vectors r i will to some extent determine the large-z behavior, but importantly must also satisfy [80]: • For all i, j, we have r i · r j = 0, which ensures linearity of deformed inverse propagators in z; • On-shellness of the shifted momenta:p 2 i = −m 2 i , which implies r i · p i = 0; • Conservation of momentum is maintained on the shift, i.e. i r i = 0 .
With an appropriate choice of shift, and for generic kinematics, and the non-trivial residues on the right-hand side are thus encoded by the kinematic poles of the amplitude. In particular, the first condition implies that the poles in A (0) n (z) are simple poles. The residues are defined by the product of lower-point on-shell amplitudes in the same theory and the scalar propagator, where L and R stand for the "left" and "right" amplitude in the factorization, and The momentum channels which contribute a residue are those which contain at least one shifted external momentum in both {p L } and {p R }, and the poles corresponding to each channel are the solutions of the linear equationŝ Note also that each pole contributes only a single residue, so partitioning into {p L } and {p R } should take into account global momentum conservation to avoid overcounting. A "good" shift on A (0) is defined as any shift for which the left-hand side of (4.3) vanishes, behavior which corresponds to the vanishing of the residue at infinity, also known as the "boundary term", For amplitudes in massless theories, it is understood what constitutes a good shift for various helicity configurations in various theories [59,[81][82][83][84]. Then, by combining (4.4) with (4.3), we get the recursive formula Later in this section we introduce a new kind of shift which is applicable to massive legs as well. In particular, it will be only the first item, the on-shellness of the momenta, that needs modification to accommodate this case.
In the following section we apply BCF shifts [85] exclusively to massless legs: they are labelled as [i, j and they modify the external legs as follows: which implies that (4.10) We now proceed to apply these shifts to graviton-scalar amplitudes.
We use the spinor-helicity formalism throughout, adopting the shorthand of [80] whereby Feynman-slashed four-momentum is replaced by the momentum labels with products denoted by simply concatenating momentum labels Differences of momenta are similarly denoted, whilst sums of momenta are combined into an upper-case P : 4.2 Building blocks of the amplitude diagrammar Figure 5. Generic on-shell diagrams that are the atoms of the diagrammar, corresponding to the amplitudes defined in (4.13) (resp. (4.14)) for (a) (resp. (b)).
We begin by looking at amplitudes with a single flavor of massive scalar, which we pick as flavor A without loss of generality. To construct these amplitudes we require pure gravity amplitudes as well as minimally coupled graviton-scalar amplitudes. The diagrams for the three-point amplitudes needed are depicted in figure 5.
The massless three-point graviton amplitudes are (4.14) where we have introduced a reference spinor χ. Although it may appear as though the amplitudes (4.14) depend on the choice of χ, this is not the case, as long as the denominators do not vanish. Using the amplitudes (4.13)-(4.14) we can apply BCFW recursion to construct fourpoint amplitudes. Up to helicity conjugation and permutation (crossing) invariance, there are two independent configurations: We can apply a [3, 4 shift to construct both, 11 In massless theories the validity of such a shift follows directly from the scaling of the twopoint propagator and polarization tensors [82,84], and this analysis appears to hold for the massive case too, as it results in the correct amplitudes (see for example [86,87]). The momentum shift involves a subset of the physical poles of the theory, and thus the amplitudes can be reproduced by the diagrams in figure 6. At four points, some simple algebra reproduces a compact form of the amplitude from the factorizations (4.20) The spurious double pole cancels upon summation with the symmetric term, and the tech- 11 For simplicity, we will suppress the a,ḃ indices from here on. nique also gives the correct result for the mixed-helicity configuration, 12 .

(4.22)
Finally, we consider the four-point two-flavor amplitude computed from a single Feynman diagram and given in equation (3.14), that is equivalent to There are well-established on-shell constraints on the classical contribution of this amplitude to eikonal scattering; it consists of a single residue in the form of a product of three-point amplitudes subject to a shift prescription which defines the residue in s 13 [88]. The full QFT amplitude requires further information to fully reproduce equation (4.23). Because of the simplicity of the Feynman diagram calculation, we treat it as a fundamental amplitude in our diagrammar, and it joins the basic building blocks in equations (4.13)-(4.14).

The equal-mass shift
The results discussed in section 4.2 relied upon the presence of massless particles in the processes in question, but here we are interested in the amplitude with two massive particles with different flavors and just a single massless graviton, as depicted in figure 7. This raises the question of whether we can construct this amplitude with any kind of shift. In fact this is possible, but first we need to consider what actually makes on-shell recursion effective. The principal advantage of the BCFW method is that it allows us to construct higherpoint amplitudes from on-shell expressions. When we are dealing with massless theories/particles, this also implies that the on-shell condition for a particle is also satisfied: p 2 i = 0. These two statements are not completely equivalent when considering theories with equally-massive particles (particles 1 and 3): an on-shell expression need not be in terms of momenta and masses which satisfy the on-shell conditionŝ but can be loosened such that the mass is shifted, but by the same value for both particles: The mass m A is thus treated like a kinematic variable rather than an invariant defining "on-shellness". Crucially, the equal-mass expressions now used in the recursion remain equal-mass expressions, and the diagrammar can be used to build amplitudes in the theory just like the massless case. This approach still requires at least one massless external particle, which we label particle 5 and assume to have positive helicity, without loss of generality. The three-line shift that satisfies the requirements of on-shell recursion iŝ where one can easily verify that Thus the condition (4.25) is satisfied, and equal-mass amplitudes can be used in the recursion. Similarly to the BCFW shift, shifting the anti-holomorphic spinor |5] produces a boundary term in A (0) 5 (z), i.e. it is a "bad" shift. From comparison with the extended-Cheung-Remmen Feynman diagram computation of section 3, we confirm that the holomorphic shift is a good shift for the five-point tree amplitude.

Five-point tree amplitude
We now apply the equal-mass shift to the tree-level amplitude with two flavors of pairs of minimally-coupled massive particles and one graviton.

The equal-mass shift we use is
(4.28) The factorization on the equal-mass poles are depicted in figure 8, and the shift yields a total of five terms, where the factorizations correspond to residues at the following poles: . (4.30) It is convenient to organize the calculation in terms of the variables 13 which are antisymmetric under the exchange of the corresponding pair of momenta. Each residue in (4.29) yields an expression containing spurious poles, which are not present in the full amplitude. For example the P 52 factorization gives with the spurious poles x ij|kl proportional to denominator factors evaluated at other residues Through algebraic manipulations the spurious poles in the full expression can be cleared, and the amplitude can be symmetrized in K A and K B . The final expression 13 Note that because the momenta are all incoming, the Ki are not momentum transfers. is and the negative helicity case is obtained simply by switching the square brackets for angle brackets. We find perfect agreement with the Feynman diagram calculation from section 3 when tested on rational kinematic points.

Six-point tree amplitude
The six-point tree amplitude can be computed using a standard BCFW shift [5, 6 , 14 where we consider which generates 10 factorization diagrams. All of these are of the general types of factorizations are shown in figure 9. We make use of the permutation invariance of the scalar particle by defining {I 1 , with the complement set labelled as J i . There are four factorizations for each of the left and middle diagrams and two for the last, giving a total of ten residue contributions to the amplitude. Schematic representation (up to crossing symmetry) of the three independent factorizations which are relevant for the construction of the six-point tree amplitude A (0) We confirm numerically the vanishing of the boundary (large-z) terms from the Feynmandiagram expression. 15 Moreover, we have verified that the reproduction of the amplitude, as the Feynman-diagram and on-shell calculations produce the same result on all (rational) numerical points tested.

The classical limit of scattering amplitudes with radiation: graviton interference is a quantum effect
In this section, we use the explicit calculation of the six-point tree amplitude A (0) 6 of the previous section to prove the coherence of the emitted semiclassical radiation field up to order O(G 4 ) for radiative observables. Moreover, assuming coherence to all orders as suggested by the arguments of section 2, we derive an infinite set of non-trivial relations between unitarity cuts in the classical limit. Those are relevant for the calculation of physical radiative observables, such as the waveform or the total linear and angular momentum emitted by the gravitons, because they suggest that only the 5-pt amplitude is required for the classical calculation and all the higher multiplicity amplitudes are not explicitly needed.
In order to take the classical limit, we follow the rules established in [23]. We express the massless momenta in terms of their wavenumbers and the momentum transfers of (2.10), q j = q j , w j = w j for j = 1, 2; (5.1) and we use the parametrization of the massive momenta from (2.12), which define the classical trajectory. They are therefore associated to classical velocities v A and v B , Note that in section 4 we used notation which was more compact for the purposes of computing the amplitudes. We can translate to the notation introduced earlier in equation (2.11) by noticing that Crucially, we also need to restore the powers of in the coupling as We use these equivalences to infer the scaling of the amplitudes. We begin by extracting the leading classical scaling of the five-point and six-point amplitude, and we then discuss the consequences of coherence for classical radiative observables.

Classical limit of the five-point tree amplitude
We begin by computing the classical limit of the five-point tree amplitude, which was given previously in [16,77] by an equivalent large mass expansion. An interesting alternative derivation can be made in supergravity theory by using the Kaluza-Klein compactification of amplitudes of massless particles in five dimensions, by taking advantage of a straightforward application of double copy [89]. The manifestly gauge invariant expression for A (0) 5 given in equation (4.34) can easily be written in terms of the polarization tensor for the graviton through the identification ∼ , (5.5) and the following scalings also hold: Moreover, we can safely neglect the quantum shift in the massesm j →0 = m j . Using equation (5.6), we can simply apply power counting to each of the terms in expression (4.34). We deduce that, upon including the contribution from κ, the terms which contribute to leading behavior as → 0 are We can make the following replacements in order to match the notation in [77] at leading order in the classical expansion, 16 This implies that the leading order behaviour of the five-point tree amplitude is of order −7/2 . As we will see later, this will imply that the amplitude contributes to the total classical energy emitted in gravitational waves. In particular, we get where f µν f ρσ is proportional to the linearized Riemann tensor and can be expressed in terms of the polarization tensor ε µν Upon substituting the relation the amplitude can thus be expressed which matches the result in [77] analytically.

Classical limit of the six-point tree amplitude
To compute the leading terms of the classical expansion of A (0) 6 , we directly extract the scaling of the BCFW residues in equations (4.36)-(4.38). In the following, we will use explicitly the rules extracted in equations (5.1), (5.4), and (5.6). First we consider the terms which originate from the factorizations of the general type (4.36), For the scaling of the three-point amplitude A 3 (I 1 ,P I ,5 ± ), we first note that a shift in momenta does not modify the scaling, which can be seen from the fact that z 5I 1 takes the form (5.15) so that it scales in the same way as p 5 . We can thus rearrange the amplitude to extract the scaling: wheref 5 ≡ f as defined in equation (5.5), but in terms of shifted momenta and polarization vectors. The shifted five-point amplitude inherits the scaling of equation (5.9), Thus upon including the contribution from the pole, each term of the form (5.13) has the leading scaling behavior We now show how taking only the leading-classical term trivializes the kinematics. Using equation (5.14) we havê and we can make the statement This is not the only simplification in the leading classical limit. We also observe that from p I 1 = −p I 2 + O( ) we have has the opposite sign under the I 1 ↔ I 2 switch, so these contributions cancel pairwise, giving A (0) An identical argument for the terms of type (4.37) also gives So the permutation invariance naturally leads to a drop in inverse-scaling. Finally, describing the scaling of terms of the type (4.38), requires the scaling from the four-point single-flavor amplitude. From counting we find And as the massless poles contributes −2 , the massless factorizations manifestly scale as Thus we conclude that We expect similar arguments to hold at higher points, which would imply that the general scaling of the (n + 4)-point amplitude is

Coherence of the final radiative state
Using the classical scaling discussed in equations (5.1)-(5.4), we can rewrite the graviton emission probability in our problem as The leading order contribution in the classically relevant region is It is the scaling of the energy of the emitted radiation that determines if the amplitude contribution is classical or quantum, and in the following we take this as a guiding principle. The expectation value of the energy operator is given by the same unitarity cuts appearing in the mean of the graviton particle distribution, but weighted in the phase space integration by an energy factor E j := ω j for each of the emitted gravitons. The scaling in the classical limit has to be such that the total energy carried by the emitted gravitons, i.e. by the classical gravitational wave, is finite in the classical limit. While each separate probability of the emission of n gravitons (5.31) is infrared divergent when λ → 0, in this paper we are interested only in the deviation from a Poissonian distribution in the → 0 limit, As we have shown in section 2, this is an infrared-safe quantity. A naive power counting in from Feynman diagrams for the five-point and six-point tree gives a series expansion starting with the following types of terms, but as we have seen in the preceding subsections, it turns out that some of the lower-order terms are zero, The cancellation of the leading term in the expansion was shown already in [77] for A 5 , but the new result here is the double cancellation of two leading terms in the expansion for A 6 . 17 This has physical consequences, as we have seen: we find that where for simplicity we have kept the powers of coming from the coupling in (5.4) implicit inside the probabilities. 18 This will be assumed for all the rest of our arguments in this section. Therefore, while the 5-point tree-level amplitude gives a classical contribution to classical radiative observables, the 6-point tree-level amplitude gives a "quantum" contribution which proves that we can describe the final semiclassical radiation state as a coherent state at least up to order O(G 4 ) for classical radiative observables.

Classical relations for unitarity cuts from all-order coherence
Assuming coherence to all orders in perturbation theory implies a set of (integral) relations between loop and tree amplitudes with emission of gravitons. 17 A similar result has been obtained in scalar QED in [53]. 18 Alternatively, we should have written We have decided to avoid this cumbersome notation here.
For example, we expect that unitarity cuts involving tree-level amplitudes with two or more gravitons emitted, and their conjugates, would give vanishing contributions in the classical limit. The reason is that having a coherent state as an exact final semiclassical state for the radiation would imply that all the gravitons emitted are uncorrelated. Indeed, our conjectural classical scaling for tree-level amplitudes in (5.29), This follows directly from our main result (2.40). While the expansion of Γ (n),λ out starts at order G 2+n , the lowest order contribution to (µ λ out ) n is of order G 2n+n : clearly then for n ≥ 2 the equation (5.39) must hold, as a simple consequence of coherence.
In order to make definite statements about the probabilities at higher orders, we need to combine them at a given loop order, so let us define we get, after imposing all the constraints (5.39) and (5.42) in the expansion in the coupling, ) .

(5.46)
We see now that the contributions in the first line manifestly involve the six-point tree amplitude and six-point loop amplitudes. We expect from (5.42), and also from the uncertainty principle [53], that the amplitudes A (L) n with n ≥ L + 2 always give a vanishing contribution to observables in the classical limit because they do not contribute to the classical field. In particular, this applies to the six-point tree-level amplitude A (0) 2 and higher-point tree-level amplitudes. But we cannot prove this directly from the coherence property, so we therefore assume that this is the case. We can then conjecture that which is equivalent to saying that the leading classical term in the expansion of the L-loop (4 + n)-point amplitudes will not conspire with the quantum scaling of the (4 + n)-point tree amplitude to give a classical contribution. It would be nice to have a direct check of (5.47) and its higher order generalizations. A first consequence of (5.45) and (5.47) is 48) which is equivalent to the statement that the seven-point one-loop amplitude is classically suppressed. More generally, from the equations (5.41) and (5.47) a very interesting set of relations follow directly for successive values of the loop order and the order in the coupling G, beginning with the following: The similarity between the top and bottom line emerges from the combinatorics of the power matching. This type of relation has a straightforward generalisation to families relating n-graviton amplitudes at L 1 + L 2 = 2n − 2 + k for k = 0, 1, 2, . . . , which can be generated from the matching of orders of G in the factorial moments. Those relations have the common feature that they relate particular combinations of unitarity cuts involving more than one graviton emitted at higher loops to other unitarity cuts involving the 5point amplitude at a lower loop level. We have represented the simplest of these relations, involving the one-loop amplitude with two gravitons emitted and the tree-level amplitude A (0) Fig. 10. Figure 10. If the final graviton particle distribution is Poissonian, this implies non-trivial classical relations between cut contributions of classical amplitudes with more than one graviton emitted to the cuts of the 5-point amplitude The outcome of this section is that we have strong evidence that the fundamental data to describe the final semiclassical state are encoded in the 5-point amplitude at all orders in the coupling constant, providing that (5.47) and its higher order generalizations hold. All the higher-multiplicity amplitudes are either suppressed in the classical regime, or related to the 5-point amplitude by a classical relation. This suggests, purely from the S-matrix perspective, that we can describe the radiation in the two-body problem entirely with a coherent state where the 5-point amplitude plays an essential role, as suggested in [53].

Conclusion
In an effective field theory approach to the scattering of compact bodies in GR, we can reduce the problem to considering a pair of minimally coupled massive scalar particles interacting with gravitons in perturbation theory. The KMOC formalism provides a rigorous framework to take the classical limit of quantum scattering amplitudes for massive particles [23], and this was recently extended to the scattering of waves by using coherent states [27]. Essentially, this is an on-shell classical limit of the in-in formalism at zero temperature which is built into the standard framework of quantum field theory. For the two-body problem in general relativity, the incoming KMOC state for two massive particles is a pure state. The unitarity of the S-matrix then dictates that ingoing pure states are mapped to outgoing pure states, and classical pure states for the radiation field are known to be described exactly by one coherent state [43].
In this paper, we have found evidence of this fact by studying scattering amplitudes with external gravitons. In particular, we have studied the properties of the final graviton particle distribution using the Glauber-Sudarshan representation [67][68][69]. We have considered the mean, the variance and higher-order factorial moments of the distribution by taking the appropriate expectation values of the graviton number operator in the KMOC formalism. Since coherent states are characterized by exact Poissonian statistics, the deviation from a coherent state structure is conveniently parametrized by the difference ∆ (m) between the factorial moments and the expected value for a Poisson distribution. Given that zero-energy gravitons in our problem obeys exactly a Poissonian distribution, ∆ is infrared finite. In the perturbative expansion, we proved that the leading contribution is related to the unitarity cut involving the six-point tree amplitude A (0) 6 (φ A φ B → φ A φ B h 1 h 2 ) and its conjugate. This is expected, since the deviation from coherence has to come from the correlation between graviton emissions. 19 The crucial problem is therefore to compute this tree-level amplitude and and its classical scaling. To do that, we developed two new approaches. First, we extend the Cheung-Remmen parametrization of the pure Einstein-Hilbert action to include minimally coupled scalars. The obtained Feynman rules are very compact, and we were able to compute analytically the full amplitude with 68 diagrams. Second, we constructed on-shell recursion relations for the case of tree-level amplitudes with two different massive particles flavors coupled to gravity: a new "equal-mass shift" is used to construct the 5-point amplitude A and the standard BCFW shift was then used to compute the 6-point amplitude A . While the large z−scaling behavior is nontrivial, a direct calculation shows that the boundary terms vanish, justifying our approach. We found perfect agreement between the two approaches, and we also agree with known results in the literature for the 5-point amplitude [77].
Regarding the classical limit, we label as "classical" the amplitudes which give a contribution to the total energy emitted in the classical limit in the KMOC formalism. The unitarity cuts appearing in such an expectation value are the same as for the probability of graviton emission, and therefore the scaling of the amplitudes appearing in those cuts determines whether we get a classical or a quantum contribution. It is known that a naive power counting does not give the correct answer, as this was already pointed out in [77] for the 5-point tree by doing an equivalent large mass expansion. Here we showed this in a manifestly gauge-invariant way in the spinorial formalism, by defining suitable kinematic variables which have a well-defined expansion. We confirmed the classical result for A (0) 5 obtained in [77], and we found that the six-point amplitude A (0) 6 gives a purely quantum contribution. A BCFW-like argument suggests that A (0) n →0 ∼ −3− n 2 , which would mean that tree-level amplitudes with an higher number of emitted gravitons should also give a quantum contribution. This result also resonates with some conjectural classical relations that we found between unitarity cuts of scalar graviton amplitudes, which point towards a characterization of the coherent state only in terms of the (all order) 5-point amplitude data. While this is often implicitly assumed in some wordline descriptions [16,28,29], our result provides a direct justification from the S-matrix perspective. Further developments along these lines have been pursued in [53]. 20 Our work has further connections in several other directions. For example, in the case of quantum field theory with external sources, unitarity cuts involving vacuum diagrams have been related to the Abramovsky-Gribov-Kancheli (AGK) cancellation in the context of context of reggeon field theory models in [51,52,94]. There it was shown that a Poissonian distribution of the cut reggeons naturally explain the AGK cancellation, and this actually inspired part of the ideas developed in this work. Furthermore, the set of infinite amplitude relations we found must have some overlap with the ones related to soft theorems [37,95], which in general are also valid beyond the classical regime. It would be nice to make this connection more precise. Finally, we have only discussed the classical perturbative longdistance regime of the scattering, but quantum and classical non-perturbative effects can make radiation incoherent and will introduce correlations between the waveform detected at different locations. It would be interesting to explore this further.
We conclude with some open questions. First, it is known that the classical description breaks down at sufficiently high energies because of quantum radiation reaction effects, which ultimately make the emitted gravitons interfere with each other [31,32,46,96]. This is actually important to have a consistent resummation of radiation reaction effects, and perhaps a simpler setup where analytic calculations are possible at very high orders -like working in a fixed background -can give us some useful lessons in this direction [97][98][99][100][101][102][103][104]. Second, we are still lacking a rigorous proof of the general validity of on-shell recursion techniques in the case of a pair of massive particles minimally coupled with gravity, which would be helpful in establishing rigorously the all multiplicity tree-level classical scaling discussed in this work. Finally, we have restricted our attention to point particles. But spin and tidal effects (and possibly other higher-dimensional operators) can also be relevant, and it is not clear if coherence will persist once those operators are added to the lagrangian. 20 Some of the ideas in this direction have been presented at the Paris-Saclay AstroParticle Symposium 2021 by P. Di Vecchia, C. Heissenberg, R. Russo and G. Veneziano in a series of seminars [91][92][93]. Figure 11. The two branches of the Schwinger-Keldysh contour C run from above (+) to below (-) the real time axis.
Using the interaction representation for the quantum fields, we can write 23 } is a set of two copies of the interaction lagrangian in the pure gravity theory where all the fields belong the same branch of the contour C. At this point we can rewrite the initial expression as where the ordering P corresponds to x 0 ∈ C + , y 0 ∈ C − .

(A.5)
We have therefore unified the treatment of the two time orderings in a compact way and have arrived at a simple path integral representation for the general in-in expectation value. Indeed, one can write a generating functional The generic SK graviton propagator in the (+)/(−) basis can then be written as where η, η can take values ±1, and P µναβ := − 1 2 (η µα η νβ + η µβ η να − η µν η αβ ) in de Donder gauge. It is manifest that we can choose any basis for the SK formulation, for example the time-ordered/anti time-ordered (also called (+)/(−)) basis as in the previous calculations or the retarded/advanced basis, and the result will be independent of that choice.
The direct connection with the standard Feynman integral perturbative expansion can be seen directly at the level of the generating functional. We can express the SK generating functional in terms of the Feynman generating functional and its conjugate Thanks to this equation, one can compute diagrams in the SK formalism by stitching together ordinary Feynman diagrams and their complex conjugates.
To make the connection with the KMOC formalism, we need to add matter coupled with gravity and to consider as our initial state |ψ in . Essentially all the previous arguments can be generalized to extend the discussion for a correlator of a set of massive scalar and graviton fields. Then we have ψ in | a † σ (k)a σ (k) |ψ in = ε µν σ (k)ε αβ σ (k) d 4 x d 4 y e ik·(x−y)/ x y ψ in | h µν (x)h αβ (y) |ψ in , we must take the LSZ reduction for the massive external states with the appropriate KMOC wavefunction ψ A (p 1 ) and ψ B (p 2 ) as defined in (2.6), ..φ(x 1 )φ(x 2 )... , (A.12) which in the limit → 0 will effectively localize the massive particles on their classical trajectories characterized by a 4-velocity v A and v B and by an impact parameter b µ . The in-in formalism is a set of off-shell techniques in QFT which in principle can be used to compute the expectation value of any quantum field or polynomial thereof, including for example the stress tensor and its conserved charges. Here we have shown that taking an appropriate LSZ reduction on the external legs and using appropriate wavefunctions for the massive particles, we naturally obtain the KMOC formalism. Under LSZ reduction, the contraction arising from time-ordered (+) or anti-time ordered (−) correlators of fields {h µν , φ} in the Schwinger-Keldysh formalism maps to S-matrix elements (with the +i prescription) and their conjugates (with the −i prescription). Moreover, the contraction of fields belonging to different branches of the contour ((+) and (−) or vice versa) gives the unitarity cut contributions. See Fig. 12 for a pictorial representation of these different contributions. We hope that this will help to address some concerns raised by Damour in [106,107] on getting classical observables from scattering amplitudes with a definite i prescription. Figure 12. The contributions of the type (a) (resp. (b)) arise from purely time-ordered fields (resp. anti-time ordered) and correspond, under LSZ reduction for the external legs, to on-shell contributions which are linear in the amplitudes. On the other hand, terms of the type (c) mix fields on different branches of the Schwinger-Keldysh contour, which corresponds in unitarity cut contributions between one amplitude and its conjugate in the on-shell formalism.
This derivation gives some insight into the relation between the SK formalism and the KMOC formalism relevant to fully on-shell calculations, like the radiated energy, angular momentum, or more localized observables like the waveform and gravitational event shapes (essentially by considering only the on-shell radiative contribution of the fields arising in the large r limit). But it also extends beyond this. In particular, it explains some recent derivations of off-shell metric configurations from "amplitudes" with one off-shell graviton leg [108]. In that case one avoids taking the LSZ reduction of the graviton field whose expectation value is taken. A simple example is given by the (linearized) metric generated by on-shell matter particles coupled to gravity. For example, this justifies the results obtained in [108] for the derivation of gravitational shock wave configurations from the 3-point function with two massless on-shell scalars and one off-shell graviton. The same argument can be repeated for any other on-shell matter configuration coupled to one off-shell graviton, essentially making use of the (linearized) stress tensor [109,110]. Alternatively, one can work fully on-shell but in (2, 2) signature, as shown in [49,111,112].

B Poissonian distributions and coherent states
The graviton coherent states introduced in the main text can be expanded in the Fock space basis of a definite number of gravitons, [dΦ(k i )α σ (k i )] |k σ 1 ...k σ n , (B.1) and a direct calculation of the probability of detecting n gravitons with helicity σ gives P σ n :=δ σσ exp − dΦ(k)|α σ (k)| 2 1 n! dΦ(k)|α σ (k)| 2 n , (B.2) which corresponds exactly to Poissonian statistics. A straightforward calculation of the mean and the variance in a coherent state gives Poissonian statistics are equivalent to having a coherent state, as can be seen by computing P σ n for a generic probability distribution in the Glauber-Sudarshan representation, Tr P σ n ρ radiation,GS = σ=± D 2 α σ P σ (α) P σ n , (B.4) which requires P σ (α) = δ 2 (α σ − α σ ) to match the Poissonian distribution. In classical physics, however, we can have more general statistics for the classical radiation field. In particular, the variance of the distribution can be greater than the mean, 24 µ ρ < Σ ρ , (B.5) which defines the so-called super-Poissonian statistics. This applies, for example, to thermal classical distributions. In our case, as discussed in the main text, the fact that we are working with pure states that are evolved with a unitary map suggests that all the classical states will have to obey the minimum uncertainty principle [23].