One-loop Gravitational Bremsstrahlung and Waveforms from a Heavy-Mass Effective Field Theory

Using a heavy-mass effective field theory (HEFT), we study gravitational-wave emission in the scattering of two spinless black holes or neutron stars of arbitrary masses at next-to-leading order in the Post-Minkowskian expansion. We compute the contributions to the one-loop scattering amplitude with four scalars and one graviton which are relevant to the calculation of the waveforms, also presenting expressions of classical tree-level amplitudes with four scalars and up to two radiated gravitons. The latter are obtained using a novel on-shell recursion relation for classical amplitudes with four scalars and an arbitrary number of gravitons. Our one-loop five-point amplitude is expressed in terms of a single family of master integrals with the principal value prescription for linearised massive propagators, which we evaluate using differential equations. In our HEFT approach all hyper-classical iterations and quantum corrections to the amplitude are dropped at the diagrammatic level, thereby computing directly contributions to classical physics. Our result exhibits the expected factorisation of infrared divergences, the correct soft limits, and highly nontrivial cancellations of spurious poles. Finally, using our amplitude result we compute numerically the corresponding next-to-leading corrections to the spectral waveforms and the far-field time-domain waveforms using the Newman-Penrose scalar $\Psi_4$.


Introduction
The extraordinary observation of gravitational waves by the LIGO and Virgo collaborations [1][2][3][4][5], 100 years after Einstein's prediction, has initiated a new era of exploration of our universe, with the promise of major discoveries in fundamental areas from black holes to particle physics. With the increasing precision and scope of current and future experiments, there is a pressing need for accurate theoretical templates for the gravitational-wave signal. A similar demand for ever more precise theoretical predictions has boosted the development, over the last few decades, of highly efficient methods to compute scattering amplitudes of elementary particles to high perturbative orders. 1 It is then remarkable that amplitudes, and modern methods devised for their computation, have now been put to use in tackling problems in classical gravity.
The connection with amplitudes was revealed more than 50 years ago in [7,8], where corrections to the Newtonian potential were computed from one-loop Feynman diagrams. Those papers also appreciated that loop diagrams contribute to classical physics, a point vigorously strengthened in [9,10]. An amplitude-based approach was applied at the second Post-Minkowskian (PM) order in [11][12][13][14] to compute corrections to the Newtonian potential using modern amplitude techniques [15,16]. More recently, several works have pursued this approach to compute the conservative part of the potential at 3PM [17][18][19][20][21][22][23] and 4PM [24][25][26], also including radiation [27][28][29][30][31], in the presence of classical spin  and in theories where Einstein gravity is modified by higher-derivative interactions [60][61][62][63][64][65][66]. The fact that gravity is a non-renormalisable theory does not prevent one from making predictions: treating gravity as an effective theory [67], non-local/non-analytic effects arising from the low-energy theory can be reliably calculated and disentangled from a yet-unknown ultraviolet completion. Observables that can be computed in this way include the deflection angle between two heavy objects (black holes or neutron stars), the Shapiro time delay, and waveforms.
Amplitude techniques have thus emerged as powerful alternatives to a variety of other approaches, such as the effective one-body formulation [68][69][70][71][72][73], and the worldline approach started in [74,75] and further developed in a relativistic setting in [76][77][78][79][80][81][82][83][84][85][86][87][88]. In these works one performs an expansion in Newton's constant G while keeping the dependence on the velocities exact (the PM expansion), which is natural from the quantum field theory viewpoint, as opposed to the Post-Newtonian (PN) expansion , also studied for spinning objects [116][117][118][119][120][121][122]. In order to have a sensible perturbative expansion one must require that GM/b ≪ 1, where M is the typical mass of a heavy object and b is the impact parameter. Quantum effects can be discarded since the characteristic Schwarzschild radius of the objects involved is much larger than their Compton wavelengths, GM ≫ /M, which combined with the previous relation requires one to work in a regime where /M ≪ GM ≪ b.
An ideal amplitude-based method tailored to compute classical observables should possess two features: first, it should easily disentangle quantum corrections from the classical contribution; and second, it should avoid the subtraction of the so-called hyperclassical (sometimes called super-classical) terms from complete amplitudes. Indeed, in the approximation we are considering GM 2 is not small, and a resummation of higher perturbative orders is mandatory. Remarkably, this is achieved in impact parameter space (IPS) [123][124][125][126][127][128], where amplitudes are believed to exponentiate, and it is precisely the IPS amplitudes that turn out to be relevant both for computing classical observables such as scattering angles and waveforms. The extraction of this "eikonal" exponent at a certain perturbative order requires a delicate subtraction of terms that reconstruct the exponentiation at lower perturbative orders; while computing such terms provides a consistency check of the result, an efficient method should preferably avoid this.
An important step in this direction was taken in [24], where the conservative part of the potential at 4PM was computed from the radial action. Motivated by the WKB formalism, in [129] an exponential representation of the S-matrix alternative to the eikonal was proposed, with a clear procedure to compute the matrix elements of the hermitian operator N, defined through S := e iN . Alas, such methods still require the computation of complete amplitudes and a subsequent subtraction.
The approach we follow in the present paper, initiated in [23,130], is based on a Heavy-mass Effective Field Theory (HEFT). Specifically, it was proposed in [23] that classical observables can be computed entirely avoiding the subtraction of iterating terms. This framework was tested in the conservative sector by deriving the scattering angle for two heavy spinless objects at 3PM, and in this paper we will show how the HEFT approach can be used to incorporate radiation emission. The key finding of [23], which makes the HEFT approach ideally suited for computing classical quantities in gravity, is that only a particular subset of diagrams contributes to the classical observables, namely those that are two massive particle irreducible (2MPI). Conversely, diagrams that are two massive particle reducible compute hyper-classical terms, and in the HEFT approach one simply drops them from the get go.
The relevance of a heavy-mass expansion arises from the fact that the momenta exchanged by the heavy particles are much smaller than their masses, thus it is natural to consider an expansion in the heavy masses. This is precisely the situation one encounters in heavy-quark effective theory [131][132][133][134]. Such a set-up in gravity was first considered in [38,47], and in [130] some of the present authors were able to combine the heavy-mass expansion with the colour kinematic/duality [135][136][137], producing compact expressions for amplitudes where the BCJ numerators are manifestly gauge invariant. All-multiplicity expressions for D-dimensional amplitudes with two heavy scalars and an arbitrary number of gluons or gravitons were then presented in [138], where the underlying BCJ kinematic algebra was also related to a quasi-shuffle algebra, further studied in [139][140][141]. Curiously, the number of terms in a numerator with n−2 massless particles is the Fubini number F n−3 , which counts the number of ordered partitions of n−3 elements (or, more mundanely, F n is the number of possible outcomes of an n-horse race, including ties). The HEFT amplitudes enjoy several important properties which make them particularly convenient as building blocks of loop integrands: in addition to the already mentioned gauge invariance, the BCJ numerators which build these amplitudes are local with respect to the massless particles, have poles corresponding to the propagation of the heavy particles, and factorise into products of lower-point ones on the massive poles.
We now come to discuss the main observable quantity of this paper, that is the gravitational waveforms produced in the scattering process between two spinless heavy objects at one loop in the PM expansion. At leading order, the waveforms for spinless objects were computed in [142][143][144] and reproduced recently in [79] (and in [80] as a onedimensional integral), while in [83,86] this was generalised to include spin. Waveforms in the frequency domain were obtained recently in the zero-frequency limit [145] and for generic frequencies [22] using amplitudes-based techniques. 2 A precise definition of waveforms in terms of amplitudes, which we will employ in this paper, was proposed in [151] and further studied in [152], based on the KMOC approach [153]. Remarkably, it turns out that the only input required to compute the waveform is the one-loop HEFT amplitude obtained from 2MPI diagrams, which is precisely what our HEFT computes efficiently. While this is perhaps not surprising, since the HEFT allows one to compute directly classical physics, this is certainly a welcome finding of our investigation.
The first goal of this work is then the computation of the one-loop HEFT, or classical, amplitude. 3 The key ingredients in this computation, which enter the unitarity cuts, are the HEFT amplitudes with two scalars mentioned earlier, but in addition we find that amplitudes with four scalars and one and two gravitons are also needed in the evaluation of a particular class of snail-like cut diagrams. In order to compute these four-scalar amplitudes we devise a novel incarnation of the BCFW recursion relation, where we shift the momentum transfers. Such shifts have the advantage of leaving unmodified the linear propagators corresponding to massive particles within the HEFT amplitudes (with two scalars), and lead to the large-z behaviour that is necessary to avoid boundary terms in the recursion.
Combining all cuts we first obtain the one-loop HEFT amplitude integrand, which is reduced using LiteRed2 [155,156]. There is an additional layer of simplicity introduced by re-parameterising the heavy momenta in the HEFT using what we will refer to as p-orm-variables. Indeed, we can reduce to just one basis of master integrals with a single iε prescription, i.e. any linear propagator will appear with a principal value prescription without the need to distinguish between ±iε prescriptions. The expression for our integrand is obtained using D-dimensional amplitudes and thus is valid in D dimensions. We then proceed to evaluate all relevant integrals around D=4 using the method of differential equations [157][158][159][160] in canonical form [161], adapted to the study of classical gravitational dynamics [19]. The integrated result for our amplitude is infrared divergent, and the divergence is in agreement with Weinberg's universal formula [162]. A number of further highly non-trivial checks on our result are also presented, including subtle cancellations of several spurious singularities.
Armed with the result for the integrated amplitude, we then proceed to compute the spectral waveform (or waveform in the frequency domain) which is, schematically, W ∼ d 4 q 1 d 4 q 2 δ (4) (q 1 +q 2 −k)δ(p 1 ·q 1 )δ(p 2 ·q 2 ) e iq 1 ·b M (1) 5,HEFT (q 1 , q 2 ; h) , (1.1) where M 5,HEFT (q 1 , q 2 ; h) is the one-loop HEFT amplitude with the emission of a graviton of momentum k and helicity h, once again obtained from only 2MPI diagrams, similarly to the conservative case. Interestingly, the spectral waveform defined above is infrared divergent, as is the one-loop amplitude; this was already noted in [150], where it was observed that in the time domain this divergence can be absorbed by a redefinition of the time variable and is thus, reassuringly, unobservable. We can then restrict to the finite part of the one-loop amplitude, and move on to evaluate numerically the two-dimensional integral that gives the spectral waveform. We do this for several values of the mass ratios of the heavy objects and as a function of the frequency ω of the emitted gravitational wave. Finally, we Fourier transform the spectral waveforms to obtain the time-domain waveforms in the far-field region. A convenient quantity to compute is the Newman-Penrose scalar Ψ 4 [163], which represents the sec-ond time derivative of the gravitational strain in the far-field region. The results of our numerical evaluations are presented in a number of plots in the frequency and time domains for various mass ratios of the heavy objects. The interested reader can find Mathematica notebooks with expressions for the HEFT amplitudes at tree level with one and two emitted gravitons, and at one loop with one emitted graviton in our Gravity Observables from Amplitudes GitHub repository.
The rest of the paper is organised as follows. In Section 2 we briefly review the fivepoint kinematics of the process at hand, introducing the parameterisation employed in subsequent calculations. In Section 3 we provide a self-contained introduction to the HEFT expansion we use, including a few illustrations thereof: the HEFT expansion of Weinberg's soft factor, and that of the gravitational Compton amplitude. Section 4 discusses our new recursion relations, providing expressions of the classical gravitational amplitudes with four scalar and up to two gravitons required in later sections. In Section 5 we perform the key computation of the paper: that of the one-loop HEFT amplitude with four scalars and one graviton using unitarity. All cuts are then merged into a single integrand, which is reduced to master integrals using LiteRed2 [155,156]. The final result for the integrand is shown in (5.40) and (5.41). The relevant family of master integrals is presented in Section 6. Section 7 discusses the final integrated result for the one-loop amplitude, shown in (7.1), along with several consistency checks of our computation. In Section 8 we move on to discuss the waveforms. We begin by reviewing the KMOC approach to such quantities, and show how to compute waveforms from the HEFT. We then calculate waveforms numerically for several mass ratios, in the frequency and time domains, illustrating our results in several plots. Section 9 summarises our conclusions and prospects for future work. A few appendices complete the paper: in Appendix A we present a detailed evaluation of the integrals listed in Section 6 using the method of differential equations; in Appendix B we review Weinberg's classic result for the infrared-divergent part of one-loop amplitudes in gravity and extract its classical part, which we use in the main text as one of the checks on our results; in Appendix C we show in an example that the exponentiation in impact parameter space of two massive particle reducible diagrams occurs also in the presence of radiation; and finally, in Appendix D we give details of one of the unitarity cuts used in our calculation.
Note added: While this paper was in preparation, we became aware of [164] and [165], which appear concurrently with our work and partly overlap with it. The key results of these papers are in agreement. We thank the authors for communication and for sharing copies of their drafts prior to publication.

Kinematics of the five-point scattering process
We are interested in computing the scattering amplitudes of two heavy scalars of masses m 1 and m 2 , accompanied by the emission of a graviton of momentum k: Here we have introduced the convenient "barred" variables [19,166], defined as (2. 2) The advantage of this parameterisation is that, using the on-shell conditions, one can show thatp i.e. the momentum transfers q 1,2 are exactly orthogonal to the barred momentap 1,2 of the heavy scalars. It is also useful to introduce barred masses, where i runs from one to the number of heavy particles (two in this case), andp i :=m ivi . As we shall see in Section 3, the HEFT perturbative expansion is organised in powers of them i . We also mention that it will sometimes be useful to write the q i in terms of the radiated momentum k and a single, average momentum transfer q as To describe our five-point scattering process, we need to specify five independent Lorentz-invariant products, which we choose as where as usual we define the four-velocities using To show that w i ≥ 0 and q 2 i ≤ 0, we simply go to a frame where v i = (1, 0) (for fixed i) and k = (ω, 0, 0, ω). Then v i ·k = ω ≥ 0, and from momentum conservation one finds that q 2 i ≤ 0 in the heavy-mass limit (this can be checked by going again to the rest frame of particle i). Furthermore, y is the relativistic factor 1 √ 1−˙ x 2 , where˙ x is the relative velocity of one of the two massive bodies in the rest frame of the other. For example, we can choose the rest frame of particle 1, where v µ 1 = (1, 0, 0, 0), and then v µ 2 = y(1,˙ x). Hence, y ≥ 1, where y = 1 corresponds to the static limit. We will also regularly use barred versions of these invariants, namelyw i :=v i ·k andȳ :=v 1 ·v 2 .
In the following, we will denote the HEFT amplitudes with two or four heavy particles (plus any number of gravitons), as A and M respectively. Complete amplitudes will be denoted as A and M. Finally, all of the amplitudes in this paper will be matrix elements of i T where we write the S-matrix as S=1 + i T .

Basics of the HEFT perturbative expansion
In this section, we give a short introduction and review of the salient features of the HEFT perturbative expansion. As an invitation to the subject, we illustrate this expansion by applying it to several examples: first, to the three-point and four-point gravity amplitudes with two heavy particles, and then to the Weinberg soft factor for the emission of a soft graviton in two-to-two scattering -the process which is the main focus of this paper. We then outline the computational strategy used to compute loop corrections to classical quantities in general relativity, following [23]. A final related application will be discussed in Appendix B in connection with the structure of infrared divergences of gravitational amplitudes at loop level.

The three-point amplitude
Our first simple example is the gravitational three-point tree-level amplitude The two massive scalars carry momenta p 1 and p 1 ′ , with p 2 1 =p 2 1 ′ =m 2 while the graviton has momentum q and polarisation tensor ε µ q ε ν q . If all of the momenta are real in Minkowski signature then momentum conservation and the on-shell conditions imply that q=0. However, this amplitude is still non-zero since the polarisation vector ε q is well-defined in the limit of zero graviton energy. Additionally, we will frequently use amplitudes like the above in unitary cuts and BCFW diagrams, where it is necessary to make the momentum q complex and non-zero.
We now wish to perform the HEFT expansion of this amplitude and therefore we use the barred variables [19,166] introduced in the previous section which satisfyp·q = 0. Furthermore, we define the barred mass and velocity as Expanding the three-point amplitude (3.2) for largem, while keeping q fixed, we find that there is only one term, which is of orderm 2 We define this O(m 2 ) term as the three-point HEFT amplitude, and we always label such amplitudes in diagrams with the letter "H".

The gravity Compton amplitude and its HEFT expansion
We now move on to the tree-level gravitational Compton amplitude, which was derived e.g. in [12].
As before the massive scalars carry momenta p 1 and p 1 ′ , with p 2 1 =p 2 1 ′ =m 2 , and the two gravitons have momenta ℓ 1,2 , with ℓ 2 1,2 = 0. Momentum conservation relates these as p 1 =p 1 ′ +ℓ 1 +ℓ 2 , and for later convenience we also introduce q:=ℓ 1 +ℓ 2 . The four-point Compton amplitude in the full theory can be written in a way that makes its doublecopy structure manifest: where the denominators and numerators are (3.8) Here (D 21 , N 21 ) = (D 12 , N 12 ) ℓ 1 ↔ℓ 2 and N [12] = N 12 − N 21 . It is also useful to note that Remarkably, one can combine the terms in (3.7) into the compact expression where the numerator is the square of the corresponding Compton amplitude in Yang-Mills, and we introduced the linearised field strength tensors F µν In order to perform the HEFT expansion of the Compton amplitude it is essential to make the Feynman iε prescription explicit, as we will see below, and to once again use the barred variables (3.10) As previously stated these satisfyp·q = 0 which avoids inconvenient feed-down terms in the expansion arising from dot products of the form p 1 ·q = q 2 /2. Using the barred variables, we can rewrite the denominators as (3.11) Again, we writep :=mv and expand the massive propagators for largem keeping the momenta of massless particles fixed: (3.13) One can then use this to expand the quantity (p 1 ·F 1 ·F 2 ·p 1 ) 2 in (3.9). The delta-function supported term gives (3.15) Combining (3.14) and (3.15), also reinstating the coupling constant dependence, we arrive at The normalisation of (3.17) and (3.18) is consistent with the three-point HEFT amplitude we found in the previous section. We now define the HEFT amplitude involving two massive scalars as the term in the HEFT expansion which is homogeneous inm and of O(m 2 ). Hence, the amplitude A 4 in (3.18) is the HEFT four-point Compton amplitude. This can also be obtained via the double copy [135,136] from its Yang-Mills counterpart, as shown in [23]. The two terms in the expansion of the Compton amplitude in (3.16) can be expressed diagrammatically as follows, where the line cut in red corresponds to the delta function πδ(p·ℓ 1 ) in (3.17). Note that the ordering of the two three-point HEFT amplitudes on either side of the red cut does not matter.
We can now make a few observations on the general structure of the expansion we have just seen in this example: two three-point amplitudes joined by a "cut propagator". This term is of O(m 3 ), and we will refer to it as the "hyper-classical term". Note thatm = m 2 − q 2 4 , hence this parameter does not have a fixed scaling, which is however recovered in the large-m limit. 5 In this terminology, the HEFT amplitude is then the classical amplitude.
2. Propagators of massless particles are untouched by the HEFT expansion, and hence are treated with the standard Feynman iε prescription. However, our HEFT amplitude contains squared linearised propagators, and it is clear from the preceding derivation (see e.g. (3.12)) that such propagators appear with the derivative of the principal value prescription 6 (3.20) For instance, in (3.13), by 1/(p·ℓ 1 ) 2 one really means the combination (3.21) We will see in Section 3.4 how these features extend to generic amplitudes. We also comment that when we reduce our integrands using integration by parts identities (IBP), we are left with a basis of master integrals which contain only a simple power of the linearised propagators. These can then be treated with the standard principal value prescription.

Weinberg soft factor and its HEFT expansion
An interesting limit of the five-particle process introduced in Section 2 is the soft limit where the graviton momentum k → 0. In [162] it was shown that in this limit the amplitude factorises into the elastic four-point amplitude (without the graviton) multiplied by the universal Weinberg soft factor. Note that this statement for the leading soft singularity holds at all loop orders in gravity, and will be used in Section 7 as a consistency check of our one-loop result.
We now want to apply our HEFT expansion to this soft factor, which will decompose into a delta-function supported term and a HEFT term. The kinematics of our 5 Further comments on the subtle distinction between the 1/m and the expansions can be found in Section 3.5. 6 This is known as the Hadamard's partie finie regularisation.
scattering of two heavy bodies with the emission of one graviton has been described earlier in Section 2. In terms of the variables introduced there, Weinberg's factor for the emission of a soft graviton has the form [162] where we have kept the Feynman iε. Next, we rewrite it using the barred variables introduced in (2.2). One then expands the denominators using and where we retain only terms up to O(m −2 i ). We also set q 2 = −q 1 := q where appropriate. Doing so one obtains delta-function supported terms along with the HEFT part of the soft factor (effectively derived by setting all the iε to zero): (3.25) A few comments are in order here.
1. First, we observe that the convenience of the barred variable stems from the fact that S W is neatly decomposes into the sum of the two terms (3.24) and (3.25). Using unbarred variables, the result would be given by the sum of the unbarred versions of (3.24) and (3.25) and in addition we would also have the "feed-down" term 3. Finally, we note that S HEFT W can be recast in the interesting forms 7 [167] S HEFT where a·F (k)·b := (a·k)(ε·b) − (b·k)(a·ε) and ε µν (k) := ε µ (k)ε ν (k) as usual. Note that S HEFT W is O(k −1 ) and linear in q.

Diagrammatics of the HEFT expansion with one heavy source
We now summarise some of the key aspects of the HEFT expansion of tree amplitudes with two massive scalars of mass m and an arbitrary number of gravitons, following the discussion of [23].
In gravity we must sum over all possible orderings of gravitons scattering off a heavy particle: for instance, we have contributions coming from the following two schematic Feynman diagrams, where the scalar propagators in the diagrams above are Switching to barred variables (3.10), as required by the HEFT expansion, with q = Q L + Q R , we can rewrite the propagators as where we usedp·q = 0. Combining the two diagrams in (3.29) using (3.32) leads to a factor of 8 (3.33) The first contribution arising from (3.29) has the form and we denote it diagrammatically as The cut line in the middle stands for the delta function in (3.34), and the two blobs represent complete amplitudes, which can then be HEFT-expanded by applying this procedure recursively. Among the terms arising from (3.29), the maximally connected one is that where the blobs are themselves HEFT amplitudes, denoted as A i+2 (1 L . . . i L ,p) and A j+2 (1 R . . . j R ,p). This contribution is of O(m 3 ), where a factor ofm −1 comes from the delta function and each of the two HEFT amplitudes is of O(m 2 ). The second term in the second line of (3.33) will instead give a new contribution to the HEFT amplitude with i L + j R gravitons.
In conclusion, a generic amplitude with two massive lines and n−2 gravitons can be expanded as where P(n−2, h) denotes the partitions of n−2 gravitons into h non-empty subsets, and the summation is over all partitions with h=1, . . . , n−2. Finally, P j denotes the j th subset of graviton indices of a given partition P with length i j and total momentum ℓ P j .
A few comments are in order before concluding this section.
1. The term with h=1 has no delta function, and is the HEFT amplitude. As discussed earlier, it is of O(m 2 ). A practical way to obtain HEFT amplitudes is to set ε=0 in all massive propagators in a generic tree-level amplitude, and then expand it, retaining only terms of O(m 2 ).

2.
A term in the HEFT expansion involving a product of h HEFT amplitudes (connected by h−1 cut propagators) is precisely of order O(m h+1 ). This is a key feature of the HEFT expansion, which allows for a clear organisation of the calculation in powers of thems.
3. The integrands constructed using HEFT amplitudes contain linearised propagators raised to powers, such as that in the last term of the second line (3.33). As was the case for four-points in Section 3.2, these are regularised using the derivative of the principal value prescription. Note that no modifications are done to the massless propagators, which follow the usual Feynman iε prescription.
4. The dots in (3.36) stand for terms that are subleading in the HEFT expansion and which are not relevant for classical physics.

5.
Finally, similar combinations of diagrams such as in (3.29) that give rise to cut propagators have appeared in several other works, e.g. [21,23,127,[168][169][170]. Importantly, in [23] and in this work, by using the 1/m expansion we have achieved a fully systematic and general HEFT expansion, free from unpleasant feed-down terms (such as (3.26)), up to and including terms that are relevant for classical physics.

vs HEFT (or 1/m) expansion
The HEFT expansion is closely related, but subtly distinct from the expansion, and here we briefly outline some of the differences, also highlighting the advantages of the HEFT expansion.
In the expansion around small values of , one scales Newton's constant as G→G/ and the graviton momenta (or sums thereof) as k → k while keeping their wavenumbers fixed [153]. Equivalently, one can take the heavy-mass limit m→∞: the graviton momenta then scale as O(m 0 ) and hence are subleading compared to any massive momentum p=mv which scales as O(m). The expansion in inverse powers of the masses is therefore equivalent to the expansion in small . This prescription is conceptually clean but has some unpleasant features. The on-shell conditions for the incoming, p i , and outgoing, . Hence, contracting massless momenta with massive ones leads to expressions that in general do not have a homogeneous degree in . As a consequence, the hyper-classical part of an amplitude will generate terms that feed down to the classical part of the same amplitude. Going to increasingly higher orders in the loop expansion, one will be faced with a proliferation of such feed-down terms from hyper-classical contributions. 9 As we have seen in the preceding sections, the natural expansion parameters in the HEFT are the inverse barred mass variables 1/m i , and each term of the expansion has a fixed degree inm i . An advantage of the HEFT expansion is that in terms of the barred variables introduced in Section 2, momentum conservation implies the exact statementp i ·q i =0. As a consequence, the HEFT expansion is free of feed-down terms from hyper-classical terms such as those appearing in the conventional expansion. Of course, the HEFT expansion still contains hyper-classical terms, for example those in (3.19). However, these are subtracted in the classical observables we wish to compute, as we shall see explicitly in Section 8.
As can be seen from their definition (2.4), the barred masses do not themselves have a fixed scaling in terms of , therefore each term in the HEFT expansion does not have a fixed scaling in . Of course an expansion in can be obtained from the HEFT expansion instantaneously by using (2.4). In other words, the HEFT expansion is a reorganisation of the expansion.
4 Tree-level amplitudes with four heavy scalars from a BCFW recursion

Diagrammatics of the HEFT expansion with four heavy scalars
To compute the classical five-point one-loop amplitude we will also require HEFT amplitudes involving four external scalars, in addition to the two-scalar HEFT amplitudes discussed above. The HEFT expansion for these amplitudes is a natural extension of the two scalar case. We start with an amplitude involving four scalars and n−4 gravitons, denoted M n , To perform the HEFT expansion of the amplitude above amplitude we expand in both the barred massesm 1 andm 2 associated with the two massive lines. The analysis of the massive propagators is essentially identical to the two scalar case in Section 3.4, except that there are now two such massive lines in every diagram.
The simplest example is the elastic process, for which we only need to keep the leading order in the large-m i expansion where + · · · are subleading terms which are irrelevant for classical physics. This leadingorder contribution scales likem 2 1m 2 2 and is what we define as the HEFT amplitude with four heavy scalars: M 4 . We will see how to calculate this amplitude in the next section explicitly, but for now let us continue with another example of the four-scalar HEFT expansion.
The tree-level five-point amplitude, with four scalars and one graviton, when expanded in the large-m i limit, contains iteration terms. These appear since the amplitude contains massive propagators which, when expanded, give delta functions exactly as in (3.33). As usual, we can write the expansion of this amplitude diagrammatically as follows where the three-point amplitudes are the HEFT amplitudes we found in (3.5). The red cut line denotes the same delta function as in the two-scalar case: πδ(p i ·k) where i = 1, 2 if the first/second scalar line is cut. The new object here is the four-scalar onegraviton HEFT amplitude, denoted M 5 , which scales asm 2 1m 2 2 and will be calculated in the next section.
For amplitudes like those in (4.1) with more radiated gravitons, the HEFT expansion proceeds analogously. There are iteration terms containing delta functions, and a term containing no delta functions which is O(m 2 term is what we define as the HEFT amplitude with four scalars and n−4 gravitons, and we denote it as M n . As before, in diagrams these are labelled with the letter "H".

The HEFT BCFW recursion relation
Here we present a novel and highly efficient method to construct HEFT amplitudes with two pairs of scalars and any number of gravitons valid in D dimensions. In order to do so we will use a D-dimensional version of BCFW on-shell recursion relations [171,172] with a carefully chosen shift that leaves unmodified the linearised propagators of massive particles, and only invokes factorisation channels that involve gravitons. The only required inputs are the HEFT amplitudes with a single pair of massive scalars available in any dimension and for any multiplicity [130,138,140], and the well-established factorisation on poles corresponding to massless propagators.
We now describe the shift using the kinematic setup and conventions introduced in (2.1) for one radiated graviton. If there are several radiated gravitons we simply replace where the k i are the graviton momenta and the corresponding polarisation vectors are denoted by ε i .
A convenient choice for the D-dimensional shifts turns out to be one where, unlike in the usual BCFW recursion [171,172], one shifts the internal momenta q 1 and q 2 : where z ∈ C and r is a null vector obeying for all gravitons i = 1, . . . , n. The first condition makes the shifted propagators scale as z −1 for large z; the second and third condition are important to guarantee that the shifted amplitude M(z) → 0 as z → ∞ as we show below. Naively, (4.6) seems to impose too many constraints on r to allow for a non-trivial solution. The way out is to demand that r lives in a space whose dimension is larger than the spacetime dimension D. For this to be useful our knowledge of HEFT amplitudes in any dimension is crucial. The solution for the shifts involves q i ·r, which need to be nonvanishing, hence also the q i must live in this larger spacetime, while the null vector r lives only in the extra dimensions. Note that our shifts differ from those employed in [173], which also addressed the computation of amplitudes with massive scalars and up to two gravitons, however before taking the classical limit.
Importantly, it is easy to see that the HEFT amplitudes appearing in the on-shell diagrams are completely unaffected by the shifts because of their structure and the conditions (4.6). This makes these diagrams particularly efficient to compute HEFT amplitudes as we will demonstrate for n=0, 1, 2 gravitons in Sections 4.4, 4.5 and 4.6.

Proof of large-z behaviour
In the following we want to show that HEFT or classical amplitudes with two pairs of scalars vanish as z → ∞ and hence there is no problematic boundary term.
We will infer the large-z behaviour from general properties of the Feynman rules and the special properties of the shift vector r introduced above. First, these amplitudes scale asm 2 1m 2 2 and subleading powers inm i will always be dropped. 10 A generic Feynman diagram contributing to a four-scalar multi-graviton process involves two two-scalar-m-graviton vertices, with massesm 1 andm 2 , respectively, and up to n multi-graviton vertices connected by graviton propagators. In order to find the large-z behaviour we only need to trace the momentum flow of the shifted momentaq i through a given diagram, where i labels the various propagators that can appear.
Let us consider a Feynman diagram with s−1 pure multi-graviton vertices and s graviton propagators with shifted momenta, connected to two vertices with pairs of massive scalars. Each propagator will contribute a factor of 1/q 2 i which scales as 1/z, and hence the propagators produce a total factor of z −s . Next there are s−1 multi-graviton vertices which are quadratic in the momenta of internal or external gravitons. Each vertex potentially contains two factors of shifted momentaq i . These either contract with k i , ε i orp i , which removes the z-dependent term in anyq i because of (4.6); or, the two shifted momenta contract with each other, which gives a factor linear in z. Hence, in the worst case one gets an overall scaling of z −s z s−1 = z −1 for a diagram with s propagators. Finally, the two-scalar-multi-graviton vertices scale asm 2 i , from two powers ofp i . Even when the shift modifies thep i , the associated z-dependent corrections are subleading in the 1/m i expansion and cannot contribute to the HEFT amplitude. In conclusion, the HEFT amplitudes have favourable large-z behaviour under the BCFW shifts introduced above for any multiplicity. Therefore we can bootstrap any classical amplitude with two pairs of massive scalars from a BCFW-like recursion of the multi-graviton Compton classical amplitudes.

Four-point amplitude: elastic scattering
This corresponds to the classical 2 → 2 amplitude without radiation. In this case q 2 = −q 1 := −q andq = q + zr. The on-shell conditionq 2 = 0 is solved by z = −q 2 /(2q · r), however this shift only appears in the polarisation vectors in the BCFW subamplitudes due to the judicious choice of shift vector r.
There is a single on-shell diagram in the q 2 -channel and the ingredients are the three-point amplitudes (3.5) and this diagram can be evaluated instantly as 11 where in the last step the sum over internal graviton polarisations was performed using 12 where d φ has the following values for the cases of pure gravity or N = 0 supergravity: where by f | v , f | η we denote replacing ε μ q ε ν q inside f by v µ v ν and η µν , respectively. 11 We remind the reader that in this paper we denote as A and M the amplitudes containing two or four massive scalars, respectively. 12 For an interesting discussion of completeness relations see [174].

Five-point amplitude: process with one radiated graviton
(4.12) In Figure (4.12) we have drawn the two recursive diagrams that contribute to the five-point HEFT amplitude, M 5 . The on-shell conditions giveq 2 1 = q 2 1 + 2z 1 q 1 ·r = 0 andq 2 2 = q 2 2 − 2z 2 q 2 ·r = 0 with the solutions z 1 = −q 2 1 /(2q 1 ·r) and z 2 = q 2 2 /(2q 2 ·r). The two diagrams give contributions of the form where the numerators N i are obtained from appropriate products of a three-point and a four-point HEFT amplitude and performing the relevant state sums. As was noted earlier in the elastic case, the three-point amplitudes are not modified by the BCFW shifts of the q i . This is also true for the new ingredient we need in this case, namely the four-point HEFT amplitude. In the first diagram of Figure 4.12 this amplitude has the form, referring to (3.6), where we are allowed to replaceq 1 by q 1 because r·k = r·ε k = r·p i = 0. However, note that the polarisation vector in F −q 1 remains εq 1 . We can also set k·q 1 = k·q, where q = (q 1 − q 2 )/2 is the average momentum transfer defined in (2.5).
With these preliminaries, we can now compute N 1 : where we have used (4.11) to perform the state sum. Similar manipulations give in terms of which the five-point amplitude with one radiated graviton is This result matches the form [175] of the classical five-point tree-level amplitude, first computed in [27]. Note that both N 1 and N 2 contain the spurious pole 2k·q which cancels when we sum the contributions from both BCFW diagrams.

Six-point amplitude: process with two radiated gravitons
In the six-point case, four distinct recursive diagrams contribute and they are given by (4.18) Once again we solve the on-shell conditions for the deformed momentaq 2 1 =q 2 2 =t 2 1 =t 2 2 = 0 to get the value of the z-pole for each BCFW diagram. We then calculate each diagram by gluing together the appropriate subamplitudes via a state sum. These amplitudes include yet another new ingredient: the five-point single heavy source HEFT amplitude, again derived in [23] This five-point HEFT amplitude can be written compactly in terms of a set of BCJ numerators [23,130,138] as follows (4.21) The first of these numerators is given by and the rest are related by permuting the massless legs ℓ 1 , ℓ 2 , ℓ 3 . Hence the six-point tree-level HEFT amplitude with two heavy sources and two radiated gravitons is where we have defined numerators N i for each BCFW diagram in the same manner as before. As was the case for the five-point amplitude, the BCFW shift inq 1 simply drops out of the amplitude for the same reasons as before.
5 One-loop five-point amplitude via unitarity

Strategy of the calculation
In this section we construct the one-loop integrand via unitarity cuts, by gluing treelevel HEFT amplitudes. The classical amplitude is obtained from two massive particle irreducible (2MPI) diagrams, which are of O(m 3 2 ). These two terms are simply related by swapping 1 ↔ 2, hence we focus here on the former. We will confirm in Section 8 that these are precisely the terms needed for the waveforms.
We also mention in passing that the O (m 3   1m  3 2 ) hyper-classical diagrams, corresponding to two massive particle reducible HEFT diagrams, factorise when Fourier transformed to impact parameter space. This was seen in the conservative case in [23], and also happens in the presence of radiation, as we show in Appendix C.
The HEFT tree amplitudes that enter the unitarity cuts have either two or four massive scalars plus several gravitons, and have been described in Section 3.4 and Section 4, respectively. They are all manifestly gauge invariant since the dependence on the graviton polarisations occurs only through the corresponding linearised field strength tensors. The three-point HEFT amplitude (3.5) is an exception which depends directly on the polarisation tensor, but it is nevertheless gauge invariant.
The cut diagrams contributing to the classical amplitude at O(m 3 1m 2 2 ) are and the following cuts which subsume the above cuts There are also HEFT diagrams where the graviton is emitted from an incoming leg, for example the following swapped version of diagram C 1 (and similarly diagram C 3 ) However this gives exactly the same contributions as C 1 , which can be seen explicitly using the loop momentum reparameterisation ℓ 1 → ℓ 3 = −ℓ 1 − q 1 and the following property of the HEFT delta function where we have used (2.3). Thus we just need to multiply the contributions of diagrams C 1 and C 3 by a factor of 2. Note that whenever we cut two gravitons as in C 1 and C 2 we also include a symmetry factor of S = 1 2! for the two identical particles crossing the cut. As noted earlier the contribution ofm 2 1m 3 2 is obtained by swapping q 1 ↔ q 2 and p 1 ↔ p 2 .
We now have to combine carefully the information from the various cuts to construct the complete integrand. Naively summing the cut integrands from all the diagrams leads to an over-counting since there are terms in the full integrand detected by more than one of the cuts above. The correct procedure, called cut merging, ensures that terms detected in several cuts are only counted once.
An example of this issue is the particular contribution to the full integrand contained in the overlap of all cut diagrams. It is easy to see that it corresponds to the following triple cut with all massless propagators present: There is also a mirror version of this diagram where the emitted graviton appears on the left of the diagram, but this makes an identical contribution.
We denote the operation of cut merging as a union and we find that the contributions in the overlap of diagrams 2 * C 1 and C 2 are also exactly identical to the triple cut diagram (5.5). Similarly, we found explicitly that the contributions detected by cut C 1 and C 2 are exactly contained as a subset of the contributions from cuts C 3 and C 4 . The exact identification and matching of overlap terms is facilitated by the use of a minimal set of independent scalar products, as we will be explain in more detail in the next section.
The integrand found in this process is given by a linear combination of tensor integrals. As we show in the next sections this bare integrand can be reduced further in a two-step process: first, we will convert the tensor integrals into a sum of loop momentum independent coefficients times scalar integrals, second, we will reduce the scalar integrals to a family of master integrals using integration by parts relations (IBP).
The relevant (master) integrals are scalar one-loop Feyman integrals of the form with the propagator structure coming from the pentagon master topology where the top line corresponds to a linearised massive propagator taken to some integer power 13 , while the bottom line is a HEFT delta function which is present in all diagrams. In the following we will refer to the propagators by the labels D i as defined in Table 1 D Table 1: Propagator basis 13 In the family of master integrals the linearised massive propagator D 5 corresponding to the top line always appears with power a 5 = 1 and is regulated with the principal value prescription. If a 5 > 1, it is regulated according to (3.21) and its generalisation to higher powers, however such integrals can always be reduced to master integrals with a 5 = 1. and in this notation the master topology is given by By using a minimal basis of scalar products and using all possible identities, we can merge integrands in different cuts before performing IBP reductions.
We will study the various cuts in more detail momentarily, but already at this stage we can identify the topologies that occur in the various overlaps of cut diagrams. The overlap terms between cut diagrams C 1 and C 2 are precisely those with the pentagon master topology (5.9), where all the massless propagators are present (possibly with higher powers of some propagators) and also the box topology where we collapse the massive propagator D 5 . The overlap terms between cut diagrams C 3 and C 4 include in addition the box topology where we collapse just propagator D 1 .
The final step is the reduction to a basis of master integrals using IBP relations. We will describe our basis in full in Section 6, where we also present explicit results for each integral.

Cut one
For the first diagram (5.1), which we have denoted by C 1 , the integrand is (5.10) In the above and for the remainder of this section we suppress the explicit Feynman iε and principal value prescriptions as they can be reinstated unambigously. The three and four-point amplitudes with two scalars and one or two radiated gravitons are given in (4.7) and (3.18).
First we perform the sum over intermediate states where for some reference momentum q, and with d φ defined in (4.10). Note that if the dilaton is included in the state sum, then the polarisation tensors ε αβ = ε α ε β are no longer traceless. Since the HEFT amplitudes are manifestly gauge-invariant, we are entitled to make the simplification P αβ → η αβ [174] in (5.11). This corresponds to the state sum shown and used already earlier in (4.9). Using the full projector (5.12) with reference momentum q, however, also allows an intermediate check since the reference momentum can be seen to drop out of the result explicitly. In addition, some care is needed when dealing with diagrams involving a three-point graviton amplitude, which is not written in terms of field strengths -an example of this situation is the diagram (5.5), and in such cases we need to use the full projector P αβ given in (5.12).
Once we have performed the state sum, the integrand is a function of where A, B can be any of the vectors ℓ 1 , q 1 , q 2 ,v 1 ,v 2 . The scalar products involving field-strength tensors F k are not all independent due to the identities [176] A·F k ·B k·C + C·F k ·A k·B + B·F k ·C A·k = 0 , where A, B, C can be any vector. The first of these relations is simply the Bianchi identity in momentum space, and the last two follow trivially from the antisymmetry of F k and the fact that k is on shell. Using these relations we can write the integrand in terms of the independent tensor structures which involve products of a pair of field strengths taken from the following list, where we have fixed the second vector contracted into F k to always bev 2 , In all of our calculations the internal cut lines are in D dimensions, while external momenta are kept in four dimensions. Four-dimensional external kinematics allows us to rewrite the tensor structure q 2 ·F k ·ℓ 1 by expanding F k in terms of a basis formed by taking anti-symmetric products of the vectorsv 1 ,v 2 , q 1 , q 2 : and then solving the linear system for the coefficients a, b, c, d, e, f . The coefficients are then written in terms of the traces and hence we will be left with an integrand made of products of these structures.
The fact that the external kinematics is restricted to four dimensions implies even more identities due to the vanishing of the Gram determinant G(v 1 , v 2 , q 1 , q 2 , ε k ). This is equivalent to the fact that that any fully anti-symmetrised tensor with more than four Lorentz indices must vanish in four dimensions. This gives for example the following and contracting relations like this with the F k ,v 1 ,v 2 , q 1 , q 2 gives us an additional relation between tensor structures that we use to simplify the integrand further. The additional relation we generate is quadratic in the field strengths and allows us to reduce the possible combinations of field strengths that appear in the integrand to two such traces, Now one can rewrite all scalar products in the numerator depending on the loop momentum ℓ 1 in terms of inverse powers of propagators (assuming the cut conditions). This gives a fully tensor-reduced integrand which can be written in terms of loop momentum independent coefficients c i and scalar integrals which are sub-topologies of the master topology (5.9), possibly with higher powers of propagators.
We then perform IBP reduction using LiteRed2 [155,156] and assuming the cut conditions of cut diagram C 1 . We find the following four master integrals, 20) and therefore the contributions are given by a sum of these integrals multiplied by their coefficients, which we write as where the overall factor of 1/2 was introduced for convenience since this contribution is doubled up when we include the swapped graph 5.3.

Cut two
For the second diagram, the integrand C 2 is given by In addition to the three and four-point HEFT amplitudes required up until now, we now also need the five-point tree-level HEFT amplitude which is given in (4.21).
The analysis and simplification of this diagram follows the same steps as for cut C 1 , however now there seem to emerge new graph topologies which are distinct to those contained in (5.8). Explicitly, propagator structures like the following can appear where ℓ 4 = ℓ 1 +k and all massive propagators are linearised. In fact, in them expansion of the HEFT, both of these topologies can be rewritten into the form of (5.8) as we now shall demonstrate.
The origin of topology T 1 can be seen by expanding the five-point HEFT amplitude in terms of BCJ numerators using (4.21) (5.25) The BCJ numerators N themselves never contain any massless propagators, hence all terms in the topology T 1 must come from the denominator 1/l 2 4 q 2 2 associated with the second BCJ numerator above. Thus, we can eliminate the topology T 2 by reparameterising the loop momentum for this particular (second) term as ℓ 1 → −ℓ 1 − q 1 , which at the level of the integrand is equivalent to the replacements ℓ 1 ↔ ℓ 3 , ℓ 4 ↔ ℓ 2 . This leaves us with the expression We would like to note that the rewriting above relies on the HEFT specific condition v 1 ·q 1 =0 imposed by the delta function and the principal value prescription for the linearised massive propagators in the HEFT. The delta function (regularised linear propagators) are even (odd) when the sign of the momentum is flipped, and this allows manipulations which otherwise would change the iε prescription.
Next, we consider the terms in cut C 2 belonging to the topology T 2 which can appear in either of the two remaining BCJ numerators where g(ℓ 1 ) is shorthand for the rest of the integrand. The topology T 2 can contain two powers of linearised massive propagators, for example, 1/(v 2 ·ℓ 1 ) 2 . However, in what follows we will assume for simplicity only single powers as the analysis in either case is the same. The first step is to perform partial fractions on the two linearised propagators to yield Next we again re-parameterise loop momentum by ℓ 1 → −ℓ 1 − q 1 , but only in the second term above, which is a subtopology of (5.8). Graphically this process can be written as and, as has been explained above, such manipulations are possible because we have principal valued and/or delta function cut linear massive propagators.
Hence, as was the case for cut C 1 , the integrand from cut C 2 contains the same basis of propagators D i given in Table 1 and corresponding to the master topology (5.9). We then follow the same method: first, tensor reduce to scalar integrals in this master topology, and second, perform IBP reduction assuming the cut conditions of diagram C 2 to get a set of master integrals. In this case these we get the box I 5 , the pentagon I 6 previously found in cut C 1 , and two additional topologies To summarise, the contributions from cut C 2 are where the coefficients c 5 and c 6 exactly match those found from cut C 1 .

Cut three -the first "snail" diagram
In the above cuts, which gave contributions C 1 and C 2 , we always cut the massless propagator D 1 with momenta ℓ 1 . There are, however, other unitarity cut diagrams relevant for classical physics that are O(m 3 1m 2 2 ) in the HEFT expansion. These are the diagrams C 3 and C 4 which only involve a single cut massless propagator and a single HEFT cut of the massive scalar on the bottom line of (5.33). These new diagrams probe all the master integrals found in cuts one and two but also involve four new master integrals with the collapsed massless propagator D 1 . We find these new contributions to be crucial for capturing the full infrared divergence of the classical five-point oneloop amplitude. That is, with these contributions the infrared-divergent part of our amplitude is given by a five-point tree HEFT amplitude multiplied by a infrared phase, which exactly matches Weinberg's prediction [162]. We explain this matching in detail in Appendix B, and now move on to computing the contributions from these new diagrams.
The first new cut we consider is C 3 in (5.2) which we replicate below for convenience which features the gluing of a four-point tree-level HEFT amplitude with two scalars into a new ingredient, the five-point tree level HEFT amplitude with four scalars. This five-point amplitude was derived in Section 4 using BCFW recursion relations and is given in (4.15), (4.16) and (4.17). The contribution from cut C 3 is then given by (5.34) Next we perform the state sum, tensor reduce the result and finally perform IBP reduction to a set of MIs. In addition to the master integrals found in cut diagram C 1 , we also have to include the following master integrals I 1 := j 01010 = , I 2 := j 01011 = , (5.35) The contributions from this cut are then given by which include the contributions from cut C 1 as expected. Once again these contributions will be doubled up when we include also the swapped graphs in (5.3).

Cut four -the second "snail" diagram
In addition to cut C 3 we also have contributions coming from cut diagram C 4 which involves gluing in the six-point tree-level HEFT. The cut diagram has the form The six-point tree-level amplitude with four scalars and two radiated gravitons was derived in Section 4. Computing this contribution is more involved than the previous cut diagrams, hence we relegate the details to Appendix D.
Cut diagram C 4 contains those master integrals probed by cut C 2 but in addition the box I 4 previously found from cut C 3 and the following new master integral The contributions from cut diagram C 4 are then given by and contain those contributions coming from cut C 2 as expected and where c 4 matches the same coefficient appearing in cut C 3 .

Final result before integration
We can now merge the contributions from all the cuts according to (5.6) and present the one-loop five-point amplitude at orderm 3 1m 2 2 in the HEFT expansion in terms of the master integrals and their coefficients: In the next section, we outline our strategy used to evaluate the master integral topologies (6.3) and present complete analytic results in dimensional regularisation (around D=4). In Section 7 we will then discuss the full, integrated result of the one-loop five-point amplitude.
14 Them 3 1m 3 6 The one-loop integrals from differential equations

The structure of the integrals
The integrals we are considering have the form We compute the master integral (MI) basis using LiteRed2 [155,156]. The full list of MIs from the merger of the contributions from the cuts in Section 5 are where we dubbed with a tilde the integrals which have pinched the massless propagator D 1 in (5.9). These are new types of integrals that appear in the bremsstrahlung process, and their contribution is fundamental for the infrared behaviour of the amplitude, as will be discussed in Appendix B.
In general, these integrals depend on the five kinematic variables (ȳ,w 1 ,w 2 , −q 2 1 , −q 2 2 ) introduced in (2.6). The strategy to compute the integrals is the following: 1. The delta function and the principal value make most of the integrals finite. Then it is easy to evaluate the integrals by direct integration of Feynman parameters, except for I 4 , I 5 and I 6 .
2. We consider the two remaining box MIs, I 4 , I 5 , and their sub-topologies separately setting up the differential equations for a single kinematic variable for each of them (ȳ andw 1 , respectively).

3.
We find the ǫ-canonical form [161] for each of these two linear systems of differential equations and solve them (details can be found in Appendix A).

5.
Finally, since we are interested in the scattering amplitude to order O(ǫ 0 ), we can write the pentagon I 6 as a linear combination of the four boxes [158], as we show later in this section.

The analytic form of the master integrals
In the following, we present the analytic expression of the MIs in dimensional regularisation up to O(ǫ): where, for convenience, we definedǫ = ǫ e (γ E −log π)ǫ , γ E is the Euler-Mascheroni constant and µ IR is the infrared scale introduced in dimensional regularisation. We have expanded the I 1 integral up to and including O(ǫ) since its coefficient arising from IBP reduction is infrared divergent. The final integral to be computed is the pentagon I 6 , which can be decomposed in terms of the boxes I 3,4,5 and I 4 [180,181]: I 5 . (6.5) The infrared-divergent part of I 6 and its imaginary part take a particularly simple form:

Final result after integration and checks
Our final result before integration was written in (5.40) and (5.41). The expression for the one-loop amplitude in (5.40) is expanded in our basis of master integrals, with the ten coefficients c i andc i potentially having spurious Gram determinant singularities at ǫ µνρσv µ 1v ν 2 q ρ 1 q σ 2 = 0. Instead, we chose to present the result in terms of the functions that appear in the integrals, which makes the analytic structures more transparent: where the coefficients d IR , R, i i , l i are rational functions of the kinematics variables and are homogeneous in the linearised field strength F µν k , and µ IR is an infrared scale. In particular, each of the coefficients is a linear combination of the various c i s andc i s and we notice these combinations are always manifestly free of spurious Gram determinant singularities. Here we are considering the amplitude in dimensional regularisation up to terms O(ǫ 0 ) and the second subscript of the coefficients specify the power of ǫ in their Laurent expansion.
A first non-trivial check of our amplitude is to confirm that the coefficients of: 2. the logarithms of the infrared scale −(lw 1 + lw 2 ), 3. the associated imaginary part i 1 are all the same, as dictated by unitarity and the Callan-Symanzik equation (for a recent discussion, see [182]): The fact that the infrared-divergent part is proportional to the tree-level amplitude can be also seen as a direct consequence of Weinberg's exponentiation of infrared divergences, which is reviewed in Appendix B: 3) The second imaginary contribution i 2 is also proportional to the tree-level amplitude: but is not accompanied by corresponding infrared-divergent or scale dependent terms. To understand the difference compared to the previous case, we need to take into account them expansion. Indeed, we have shown in Section 3 that this expansion implies the use of the principal value prescription for the linear propagators, which in turn has the effect to make the otherwise divergent integrals I 2 , I 3 , I 3 and I 4 finite, without altering their imaginary part. Then, i 2 in (7.4) should be thought of as the imaginary part corresponding to infrared divergent/scale dependent contributions that have been moved to iteration terms 15 .
Moreover, the amplitude (7.1) has spurious singularities in the physical region which have to cancel when we combine different logarithmic contributions near the poles. First of all, let us emphasise that the argument of the logarithm and the square roots appearing in I 1 and I 2 are positive in the physical region, as discussed in Section 2.
At leading order in the soft limit, i.e. O(ω −1 ), only the two-massless triangles 15 This can be checked explicitly by expanding Weinberg's soft phase in terms of the p i variables, instead ofp i to order m 3 1 m 2 2 or m 2 1 m 3 2 .
contribute, and they reproduce exactly Weinberg's soft factor in the HEFT (3.28): where q ≃ −q 1 ≃ q 2 as in Section 3.3 and M is the classical four-point one-loop amplitude: Finally, we notice that: 1. lw 2 and lȳ have a double pole atŷ 1 =w 3. c 1 , c 2 and l q have a pole of order four at 16ŵ 4. the rational terms R have single poles in bothŷ 1 andŷ 2 and a pole of order three inŵ 1 .
For example, nearȳ ≃ŷ 1 , the logarithms simplify as log ȳ 2 − 1 +ȳ (7.7) and one can show that the double poles and the logarithms cancel when we combine lw 2 and lȳ. Nevertheless, the leftover term has still a simple pole inȳ ≃ŷ 1 . This pole is only cancelled once we take into account the rational contribution R, which is needed to restore locality and cancel spurious poles [183,184]. The spurious singularities in y ≃ŷ 2 andw 1 ≃ŵ 1 share the same fate, even though showing it explicitly is more involved because the terms to be combined are more complicated.
As a final check, the present authors and those of [164] have performed independent comparisons of their two results for the one-loop amplitude, finding perfect agreement. 16 The combination −4q 2 1w 2 1 − q 2 1 − q 2 2 2 is non-negative. Indeed, if we choose the frame wherē v 1 = (1, 0, 0, 0) and k = ω(1, 0, 0, 1), then this becomes 4 ω 2 (q 2 1,x + q 2 1,y ). This means that the poles sit in the physical configuration for which k, q 1 and q 2 are taken to be aligned (modulo boosts), which is expected to be smooth. Likewise, one can show thatŷ 1 corresponds to k being orthogonal to q 2 in the rest frame of particle 1.

Blitz review of the KMOC approach
In this section we review the connection between waveforms and amplitudes following the KMOC approach [151][152][153]. The two heavy objects are taken to be in an initial state represented as where the wavefunctions φ(p 1 ) and φ(p 2 ) are peaked around the classical values of the momenta. We use the same conventions as [151], for the massive objects. Similarly, for gravitons we choose where h denotes the helicity. The waveform we are interested in is related to the expectation value of the Riemann tensor [151][152][153], where S = 1 + iT and R µνρλ (x) is the Riemann tensor, evaluated at the position x of the observer, in the far future of the event. Expanding R µνρλ (x) as 5) one finds that in ψ|R µνρλ (x)|ψ in = 0, and [151][152][153] in ψ|R µνρλ (x) T |ψ in where we introduced barred variables in the delta functions as in (2.2). In the last equality we have also changed integration variables from Next we consider the term in ψ|T † R µνρλ (x)T |ψ in in (8.4). It can be rewritten in a similar form to the previous one noting that Following identical manipulations as before and combining the two non-vanishing contributions from (8.4), we arrive at the following expression for the expectation value of the Riemann tensor: Of course one can also follow the same procedure for h µν (x). With the free-field expansion Note that R out µνρλ (x) ψ or h out µν (x) ψ are effectively computed using where from now on we will drop the subscript "in" in the state |ψ .

From KMOC to HEFT
Our next task is to compute the quantity that appears in (8.9),

13)
17 In our conventions the linearised Riemann tensor is at one loop in the PM expansion. The first term is the complete amplitude, where we note that in our conventions where we have also indicated the dependence on the polarisation h of the graviton. Next we insert a sum over intermediate states in the second term in (8.13), The first intermediate state contributing at one loop is |n = |r 1 r 2 with r 1 , r 2 being two scalars, where the superscripts in (8.16) denote the loop order. This contribution is (after an expansion inm i ) nothing but a two massive particle reducible diagram, shown below in (8.17): This contribution is of O(m 3 1m 3 2 ) and is thus hyper-classical compared to the classical one-loop amplitude computed in Section 5. Therefore subtracting the expression in (8.15) from (8.13) has the effect to peel the hyper-classical contribution off the complete one-loop five-point matrix element. In our HEFT approach, the subtraction of such terms is achieved by simply dropping all two massive particle reducible diagrams, very similarly to what was done in [130] in the elastic case; in other words, the HEFT directly computes the classical part of (8.13). This is clearly one of the strengths of the HEFT.
One could also consider the case where |n = |r 1 r 2kh , with an additional inter-mediate graviton but these all give a vanishing contribution because of the kinematics. Finally, there are also hyper-classical iteration diagrams such as those in (4.3) involving a three-point HEFT amplitude with an external graviton. These terms are also not 2MPI and hence are also subtracted by pieces in (8.16), this time involving a disconnected five-point amplitude. However, even if these terms were not subtracted, they only involve zero-energy gravitons and hence do not contribute to R out µνρλ (x) ψ , but would to h out µν (x) ψ . Two brief comments are in order here. First, we note that hyper-classical contributions such as (8.16) exponentiate in impact parameter space, which we show explicitly in Appendix C. We also note that a similar one-loop cancellation in the expression of the waveforms at that order was advocated in [152], which studied generalisations of the eikonal in the presence of radiation [145,185,186].
Summarising, we have found that the quantity of interest is directly the one-loop matrix element computed in the HEFT from 2MPI diagrams, which we calculated in Sections 5 and 6, and summarised in Section 7: where M 5,HEFT is given in (5.41). We can then use this result to evaluate the earlier expressions (8.9) and (8.11) for the expectation value of the Riemann tensor and the gravitational field.

From HEFT to waveforms
Combining (8.9) and (8.11) with (8.18), we can write the one-loop expectation value of the Riemann tensor or the metric in terms of the 2MPI five-point HEFT amplitude: and (8.21) At one loop, this expression reduces to and we have defined As usual, k h denotes a graviton with helicity h=±. Note that having eliminated any hyper-classical terms from the waveform, we are now free to express the HEFT amplitude in terms of unbarred variables since any feed-down terms will be quantum.
is directly related to the waveforms, but before making this connection more precise we would like to pause and make a few comments:

Dependence on b
First, the dependence of (8.22) on b can be made more explicit: changing variables from (q 1 , q 2 ) to the variables (q, k) introduced in (2.5), we can rewrite (8.22) as We will henceforth set b 2 =0 and b 1 =b and drop the overall phase in (8.24).

Infrared (in)finiteness of the gravitational waveform
The one-loop HEFT amplitude M 5,HEFT in (8.22) contains infrared divergences, which in gravity give a non-vanishing contribution to the waveform. This is in agreement with earlier computations performed in the PN expansion [75]. While such divergent phases drop out of quantities such as cross sections, they are still present in the waveform, which is linear in the classical (or HEFT) amplitude. The question then arises as to what is their fate. The answer was suggested in [150], where the waveform in the timedomain was considered and it was noted that the divergent phase simply shifts time (in the exponential e −ik·x in (8.9)) by an amount proportional to 1/ǫ. Using the classical limit of Weinberg's formula computed in (B.10), we see that this shift has the form with k µ =ω n µ . Physically this time shift is not relevant as ultimately we only deal with time differences -we measure time with respect to when an experimenter begins tracking the wave signal. As a consequence, we can safely drop the infrared divergences in the one-loop amplitude.
It is also interesting to note that, by contrast, the electromagnetic analogue of the classical waveform is free of infrared divergences. Indeed, the infrared-divergent Coulomb phase in QED [162], evaluated in the HEFT expansion is e i 4πǫ e 1 e 2 y √ y 2 −1 +O(m −2 i ) , and only has hyper-classical contributions. The absence of classical infrared divergences in QED has a very clear physical interpretation: the photon does not interact electromagnetically with the system at large distances, while the graviton is influenced by the total (ADM) mass of the binary system. For example, this difference is crucial when computing the waveshape of [152] at next-to-leading order, which is finite in electrodynamics but infrared-divergent in gravity.
In conclusion, we can just focus on the simpler four-dimensional integral 5,HEFT,fin (q 1 , q 2 ; h) , (8.27) where M (1) 5,HEFT,fin is the infrared-finite part of the expression given in (5.41), and dµ (4) is the measure introduced in (8.23) evaluated for D=4. We will then safely use W in (8.27) within (8.19) and (8.20) instead of W .

Waveforms and Newman-Penrose scalar from the HEFT
Several quantities can now be introduced to describe the waveforms. One that is commonly used in the study of gravitational waves is the Newman-Penrose scalar [163] Here W is the Weyl tensor, in our case equal to the Riemann tensor, and where ζ is a reference vector chosen such that ζ·ε (±) =0, and ζ·n=1. Ψ 4 (x) is often used to illustrate the gravitational waveform [187][188][189], and is also the quantity considered in the open-source numerical relativity code GRCHombo [190,191].
Following standard steps, we can now compute Ψ 4 (x). First, using (8.5) and (8.29) the dotting gives where as usual the ε (±) µ vectors satisfy Note that the matrix elements that appear in (8.30) are of course the same as in (8.9), which (as summarised by (8.18)) correspond to the 2MPI one-loop five-point amplitude.
Next, one introduces retarded time u:=t−r, where r:=| x| is the observer's distance, and rewrites the plane waves in (8.30) using k·x=ω(t−rx·n)=ωu+ωr(1 − cos θ). At large r the exponentials oscillate very fast, and a well-known stationary phase approx-imation argument [192] gives where f ( k)=f (ω, ωn) is a function of the graviton momentum k; note that after the minimisation, the directionn of k is aligned to that of x. We then introduce Ψ 0 4 (x) from 33) to find, in the far-field domain, where we remind the reader that u is the retarded time. This is the result we will use to compute waveforms in the time domain.
Three final comments are in order here.
1. We recall that W was defined in (8.27). It is constructed out of the finite part of the 2MPI HEFT amplitude, which contains only classical physics.
2. Furthermore, we observe that at tree level which follows from the form of the tree-level five-point amplitude and the definition of W in (8.27). 19 From now on we drop the integrations 2 j=1 dΦ(p j ) |φ(p 1 )| 2 |φ(p 2 )| 2 as we are assuming that the wavefunctions φ(p i ) are peaked around the classical value of the momenta of the heavy objects, and are furthermore simply spectators in the evaluations.
3. Finally we comment that a quantity widely used to characterise the waveforms is the gravitational strain h, of which Ψ 4 is the second derivative with respect to the retarded time u, Ψ 4 =d 2 h/du 2 . This can also be obtained from our previous formulae: (1,x) .
In the next section we will perform numerical integrations and will present various plots of W (b, k ± ) and ω 2 W (b, k ± ) which will illustrate the waveforms in the frequency domain; we will then move on to show the waveforms in the time domain as obtained from (8.35). Note that from now on we will refer to W (b, k ± ) simply as the spectral waveform.

Set-up of the integration for waveforms
In this section we address the computation of the one-loop waveform introduced in (8.27), which enters the Newman-Penrose scalar Ψ 0 4 (8.34). A convenient way to perform the integrations in (8.27) was discussed in [151]. After integrating out q 2 using the delta function, one can parameterise the remaining integration variable q 1 as (renaming q 1 → q for notational simplicity), where and Choosing b to be the asymptotic impact parameter, we also have that b·v 1 = b·v 2 = 0. The Jacobian is then d 4 q = y 2 − 1 a=1,2,v,b dz a , so that The delta functions set We also note that

Waveform for binary scattering
We now move on to evaluate the waveforms numerically using our one-loop result for the HEFT amplitude. For completeness we will first briefly review the tree-level waveforms, before considering the one-loop case. In the following we parameterise the kinematic data as therefore working in the rest frame of the first heavy object. We also choose the impact parameter as b= √ −b 2 (0, 0, 1, 0). The polarisation vectors ε (±) are related to the graviton polarisation tensors corresponding to positive/negative helicity as where we have also written the relation of the polarisation tensors of ± helicity to the two standard transverse-traceless (TT) polarisation tensors plus (+) and cross (×). 20 The z b integration can, in principle, be performed analytically by closing the integration contour in the lower-half plane. Therefore, this integration can be rewritten as a sum over residues on poles and integrals over discontinuities on branch cuts using Cauchy's theorem. To compute the waveform in the time domain, it is convenient to perform the ω integration first, which evaluates to (derivatives of) a delta function and a PV of u − z b √ −b 2 for the real and imaginary part of the amplitude, respectively.

Tree level
At tree level the waveform is given by 2 term in the HEFT expansion which was computed in Section 4.5. W ± is shorthand for W (b, h ± ), and the subscript m 1 m 2 indicates the mass dependence of the corresponding term. We have also introduced which parameterises the relative mass ratio of the two massive objects. The hyperclassical terms in (4.3) are subtracted in the calculation of the waveform, as explained in Section 8.2. As mentioned previously, this allows us to work with m i and v i instead of their barred versions, since any feed-down terms we generate will be quantum.
At tree level the waveforms depend on the mass ratio χ only through the prefactor χ(1 − χ) in (8.47). Thus, they are maximised when both masses are equal, for a given total mass. As such we have only plotted the equal-mass case in Figures 1 and 2.
Finally, we mention that tree-level waveforms for non-spinning objects were derived in [79,144] in the time domain, see also [22] for a derivation in the frequency domain and [80] for a one-parameter integral representation of the time-domain waveform.

One loop
In the present section we evaluate numerically the following quantity, where the subscript "fin" means that we are dropping all infrared divergences in the corresponding amplitudes, as discussed near (8.27). This corresponds to isolating the new contributions of the one-loop amplitude from the lower order in the PM expansion and the tails [146,150]: As before, W ± is shorthand for W (b, h ± ) and the subscripts m 2 1 m 2 , m 1 m 2 2 indicate the mass dependence of the corresponding terms. In the plots displayed below we will show results for several values of χ.
We begin by plotting the quantity ω 2 W ± which is the spectral version of the Newman-Penrose scalar (8.35), for various choices of χ. For the positive-helicity waveform at θ → π 4 , φ → π 2 , and y → 2 this is shown in Figure 3, where we stripped off a dimensionful factor of κ 5 (m 1 +m 2 ) 3 (4π) 4 (−b 2 ) 3/2 . The corresponding negative-helicity waveform is shown in Figure 4. In the frequency domain, the most interesting part of the spectrum is contained in the region 20]. Beyond that, the amplitude is very small and tends to zero as ω → ∞. At one loop, the dependence on the mass ratio χ follows the same pattern as at tree 21 From the field theory viewpoint such exponentiation is natural. Indeed, we know from [162] that the infrared divergences exponentiate as per (7.3) and we expect them to be accompanied by an infrared-running logarithm, i.e. schematically 1 22 Recall that dµ (4) is proportional to (m 1 m 2 ) −1 .
-60 - level, in that the equal-mass case has the largest waveform due to the prefactor χ(1−χ) in (8.52). However, due to the two terms in (8.52) the χ-dependence of the waveform is not as simple, as Figures 3 and 4 show. This is a new feature at one loop. The time-domain waveform is obtained using (8.35) by performing a numerical Fourier transform. For example, in the equal-mass case, corresponding to χ = 0.5, the result is shown in Figure 5 (up to an overall factor of κ 6 (m 1 +m 2 ) 3 (4π) 6 (−b 2 ) 2 ). The small oscillations are due to the use of a finite frequency domain in the numerical computations and vanish when the range of frequencies is enlarged.

Conclusions
The HEFT approach provides a powerful method to compute classical effects both in the conservative sector [23] and in the presence of radiation. It has the clear advantage of computing directly classical quantities, leading to integrals with linearised propagators and well-defined iε prescription that can either be computed directly or using differential equations. It holds the promise to be efficiently applicable to higher PM orders and to other problems.
Several issues are left to investigate. Having computed the graviton emission amplitude in a scattering process, it is important to determine the corresponding one for bound states. It would be remarkable if an analytic continuation similar to that of [193][194][195] or a generalisation of the Bethe-Salpeter equation approach of [196] could be applicable also to waveforms.
A number of concrete problems can also be tackled with our method and results: one could determine the one-loop waveforms fully analytically at one loop, or study the radiated energy, power and angular momentum. Going in a different direction, there are intriguing differences and complementarities between the HEFT approach initiated in [23] and pursued in this paper, and the eikonal approach [22, 123-127, 145, 152, 185, 186, 197, 198]. In the former, which appears to be intimately related to the N operator discussed in [129], experience so far indicates that classical contributions can be computed directly without pollution either from quantum or hyper-classical terms, which can be discarded at the diagrammatic level. Clarifying the relationship between these approaches would be highly desirable. We hope to come back to some of these questions in the near future.

A Integrals from differential equations
A.1 The differential equation for j a 1 ,1,a 3 ,a 4 ,0 with respect to w 1 The subset of MIs which appear in our basis as sub-topologies of j 1,1,1,1,0 are and the differential equation with respect to w 1 in D = 4 − 2ǫ looks like If we normalise the integrals in the basis through the transformation the system of differential equations takes the canonical form (A.8) In order to make the singularities of the integrals manifest, we need to rationalise the system of differential equations and write it in dlog forms. We define with α ≥ 1. Then, the system of differential equation can be written in term of forms as where we omitted any dependence on the kinematic variables other than α, and The solution to the differential equations are given by where j 1 (α 0 , ǫ) is a chosen boundary value of the Feynman integral which needs to be fixed, Pe is a path-ordered exponential.
Finally, we are left with the evaluation of the boundary value of the I 5 integral and we choose to compute its asymptotic behaviour near the singular point α ∼ 1 (w 1 ∼ 0), using the geometric approach to the method of regions [178] implemented in the Mathematica package asy2.1.m [179]: where in the second step we performed the x 2 integration, redefined ( keeping only the leading term in the α ∼ 1 expansion and used the Cheng-Wu theorem [199] to make the integration on x 4 trivial. The ǫ-expansion in D = 4 − 2ǫ gives us the boundary value of the integral we are after. A.2 The DEs for j 0,1,a 3 ,a 4 ,a 5 with respect to y The subsets of MIs which appear in our basis as sub-topologies of j 0,1,1,1,1 respectively are 14) and the differential equations with respect to y in D = 4 − 2ǫ are If we perform the change of basis with the systems of differential equations take the canonical form. Moreover, we can rationalise the system of differential equations and write it in terms of dlog forms through the following change of variables: Then, we find where we omitted any dependence on the kinematic variables other than x and At this point, we only need a boundary value of the I 4 integral. We choose to compute its value around the regular point x ∼ 1 (y ∼ 1):

B Infrared divergences and heavy-mass expansion
In this appendix, we review the classic result of [162] for the infrared divergences in gravity and consider its limit in the large mass expansion.

B.1 Weinberg's formula for infrared divergences of gravitational amplitudes
In [162], Weinberg presented a compact formula for the resummation of infrared divergences in gravitational amplitudes arising from the exchange of virtual soft gravitons, in addition to re-deriving similar formulae for photons in electrodynamics, reproducing (and in part upgrading) the work of [200,201].
There are two types of contribution: first, each pair of particles in the initial or final state contributes an infrared-divergent phase e −iW ij log( λ Λ ) , where Here is the relative velocity of any one particle in the rest frame of the other, and λ is an infrared cutoff (see below for a translation to dimensional regularisation). This phase is usually discarded in the computation of observables such as cross sections but will be important for us. In addition there is a divergent contribution of the type e i,j B ij log( λ Λ ) , where where η i = ±1 depending on whether the particle is outgoing (+) or incoming (−). In this case the sum is over all pairs of particles, including the case i = j.
If one of two particles in a pair is massless, say particleî, then βˆi j → 1. In this case Weinberg's formulae simplify to with the result being independent of the choice of µ after summing over j; indeed, a shift in µ 2 changes j Bˆi j by an amount proportional to pˆi· j η j = −ηˆip 2 i = 0. Note that we can combine Wˆi j and Bˆi j into one quantity valid for all kinematic regimes, by replacing the exponent by where as usual log − (s + iε) = log(s) − iπ for s > 0. 23 We also comment that to express Weinberg's formula in dimensional regularisation. we need to make the replacement Note that both sides of (B.6) are negative.

B.2 The large-m expansion of Weinberg's formula
We now apply Weinberg's formulae to the process we are describing. We will first discuss the phase and then the real contribution.
We will use (B.3), and compute the quantities β ij for the various cases. Because we want to perform a HEFT expansion, in order to be able to compare to the result of our calculation we need to expand p 1 ·p 2 and p ′ 1 ·p ′ 2 aroundp 1 ·p 2 , that is an expansion in m i . We begin with the contributions from the pairs (1, 2), (1 ′ , 2 ′ ). These take the form hence the contribution from pairs (1, 2), (1 ′ , 2 ′ ) is simply where W12 is obtained from (B.1) by replacing p 1 , p 2 withp 1 ,p 2 . When we expand the exponential in powers of G, its contribution will be of orderm 3 1m 3 2 in the large mass expansion. Then, there is no need to compute it -as we have explained in Appendix C, such hyper-classical contributions are obtained simply from the exponentiation of lowerorder amplitudes.
Next we consider the pairs (1 ′ , k), (2 ′ , k). Because one of the particles is massless we can use the first of (B.4). The result is then simply 2G(p ′ 1 + p ′ 2 )·k = 2G(m 1w1 +m 2w2 ), with the corresponding contribution to the phase being e −2iG(m 1w1 +m 2w2 ) log( λ Λ ) , (B.9) which, once expanded, will appear in the large-mass expansion of the amplitude at the order we are considering. Translating to dimensional regularisation, the expected infrared divergence at one loop is iG(m 1w1 +m 2w2 ) ǫ M

C Factorisation in impact parameter space
The purpose of this section is to show that also in the presence of radiation, HEFT diagrams that are two massive particle reducible factorise in impact parameter space. 24 To this end, consider the one-loop five-point diagram with two cut massive lines shown 24 See also Section 4 of [152] for a related discussion.
In conclusion, when transformed to impact parameter space, the particular two massive particle reducible diagram considered in (C.1) is a product of two tree-level amplitudes in impact parameter space. This shows that factorisation is made manifest by the HEFT expansion in impact parameter space also in the radiative case.
A short comment is in order here. In the four-point case (C.4), the delta functions impose that q lives in the (D−2)-dimensional subspace orthogonal top 1 andp 2 , hence M 4 ( b) will depend only on the projection b ⊥ of b living in the same subspace. For the case of M 5 ( b) in (8.24), the particular orthogonal subspace is slightly different because of the presence of k within the delta functions. This is a general feature of five-point kinematics and beyond.

D Details of the C 4 calculation
Here we present the details of the computation of C 4 . Starting from the cut diagram (5.37), if we simply plug in the HEFT amplitudes as before we obtain (D.1) An unpleasant feature of the above integrand is that it contains divergences of the form which come from the linearised massive propagators present in the six-point HEFT amplitude.
To deal with these divergences we perform the following steps: 1. First we consider C 4 without the delta function δ(v ! ·ℓ 1 ) and then expand the integrand in powers ofv 1 ·ℓ 1 .

2.
Then we only keep the (v 1 ·ℓ 1 ) 0 term in the expansion which has no divergences.
Despite this procedure being rather ad hoc, the resulting expression for C 4 merges exactly with the other cut diagrams. In addition, the resulting amplitude built using C 4 passes many nontrivial checks (see Section 7) including the cancellation of spurious poles involving the new contribution from this cut diagram.
To make the calculation of C 4 more rigorous here we outline how to calculate C 4 using the forward limit, which will turn out to be equivalent to the simplified procedure above. First, we find the cut tree-level amplitude with three pairs of massive particles E(q 1 , q 2 , q 3 ,p 1 ,p 2 ,p 3 ) = For each massive line, we define a momentum transfer q i := p i −p i ′ and hence momentum conservation for this amplitude can be written as q 1 +q 2 +q 3 = k. Both massive particles p 1 and p 3 have the same mass p 2 1 = p 2 1 ′ = p 2 3 = p 2 3 ′ = m 2 1 , and we can define barred masses in the usual way:m i := m 2 i − q 2 i /4. Thus the cut we are considering above is a cut in q 3 , and is homogeneous in the masses with scalingm 2 1m 2 2m 2 3 . Once we have this cut amplitude we take the forward limit on the two massive lines p 1 ′ and p 3 to form the loop diagram in (5.37). This method was also used in the worldline formalism [28] to compute classical radiation in dilaton gravity and Yang-Mills theory. The forward limit process is most clearly described by writing the amplitude in impact parameter space, which we obtain by performing a Fourier transform with respect to each q i 3 i=1 d D q i e iq i ·b i δ(p 1 ·q 1 )δ(p 2 ·q 2 )δ(p 3 ·q 3 ) δ (D) (q 1 + q 2 + q 3 − k) E(q 1 , q 2 , q 3 ,p 1 ,p 2 ,p 3 ) .

(D.7)
Note that this removes the appearance of q 3 in the momentum-conserving delta function. Next we take the forward limit by taking p 3 → p 1 ′ which amounts to the replacementp Since we are identifying particles 1 and 3 we will also identify their impact parameters b 3 → b 1 . Applying the forward limit we find which crucially is not singular. However, the above expression is no longer homogeneous in the massm 1 = m 2 1 − q 2 1 /4 and must be re-expanded in the large-m 1 limit. Before performing this expansion we note that now the integrals over q 1 and q 2 simply transform the HEFT amplitude with four scalars to impact parameter space, and hence can be stripped off, leaving us with 25 d D ℓ 1 δ p 1 ·ℓ 1 + ℓ 1 2 ·(ℓ 1 + q 1 ) δ (D) (q 1 + q 2 − k) E − ℓ 1 , q 2 , q 3 ,p 1 + q 1 2 + ℓ 1 2 ,p 2 ,p 1 + ℓ 1 2 , (D.10) where we have reparameterised the momentum as q 3 = −ℓ 3 = ℓ 1 + q 1 to make contact with our usual loop integration variables in (D.1). To recover C 4 the final step is to perform the large-m 1 expansion on (D.10). However, we will perform this expansion in a specific way in order to make contact with the simplified procedure for computing C 4 described above: 1. First we only expand E asm 1 → ∞ and never use the delta function in (D.10) to simplify the expression. To leading order, this yields the term E(−ℓ 1 , q 2 , q 3 ,p 1 ,p 2 ,p 1 ) which is O(m 4 1m 2 2 ). This is exactly the naive integrand in (D.1), without the delta function.
2. Next, we use the delta function in (D.10) to replace all instances ofp 1 ·ℓ 1 in E with ℓ 1 ·(ℓ 1 + q 1 )/2, which results in many terms with differing powers ofm 1 . When this replacement happens to a denominator, for example, 1 p 1 ·ℓ 1 = 2 ℓ 1 ·(ℓ 1 + q 1 ) , (D.11) the power ofm 1 is increased and a spurious massless pole is introduced. These poles must cancel, which is why expanding E to only the leading order is justified.

4.
Finally, we re-expand the expression δ(p 1 ·ℓ 1 )E in powers ofm 1 to find that the leading term is of O(m 3 1m 2 2 ), as expected.
The result of this process is exactly the same as the naive method described at the start of this appendix. As was the case for C 2 , we now need to perform additional loop momentum reparameterisations in order to write C 4 in the topology (5.8).
As a final remark, we note that it would be possible to compute the entirety of M in one fell swoop following the above method if we had instead started from the full tree-level six-scalar one-graviton amplitude at orderm 2 1m 2 2m 2 3 . This would bypass the use of generalised unitarity completely. One can compute this tree-level amplitude using the BCFW method presented in Section 4 and generalising it to six scalars. However, unsurprisingly, we found this amplitude contains a very large number of terms, so we find it more practical to split up the calculation using multiple cut diagrams and generalised unitarity.