Testing subleading multiple soft graviton theorem for CHY prescription

In arXiv:1707.06803 we derived the subleading multiple soft graviton theorem in a generic quantum theory of gravity for arbitrary number of soft external gravitons and arbitrary number of finite energy external states carrying arbitrary mass and spin. In this paper we verify this explicitly using the CHY formula for tree level scattering amplitudes of arbitrary number of gravitons in Einstein gravity. We pay special care to fix the signs of the amplitudes and resolve an apparent discrepancy between our general results in arXiv:1707.06803 and previous results on soft graviton theorem from CHY formula.

Our main interest in this paper will be multiple soft graviton theorem -the case where multiple gravitons carry soft momentum. Even though, at the leading order in soft momenta, this was studied in the original paper of Weinberg [1,2], much of the analysis at the subleading order has been restricted to the case of two soft gravitons [45][46][47][48][49][50][51][52]. In a JHEP01(2018)090 previous paper [53] we gave a general formula for amplitudes with multiple soft gravitons in terms of amplitudes without soft gravitons in any quantum theory of gravity, for arbitrary mass and spin of the external states, following the general method developed in [29][30][31]. This result may be stated as follows. Let us suppose that we have n external finite energy particles carrying polarizations and momenta (ε a , p a ) for a = 1, · · · , n, and m soft gravitons carrying polarizations and momenta (ε n+r , τ k n+r ) for r = 1, · · · , m. ε a for 1 ≤ a ≤ n can be any tensor or spinor representation of the Lorentz group. Then for small τ the amplitude M n+m may be expressed in terms of the amplitude M n without the soft gravitons as follows: 1 = (p a · k) −1 (p a · k ′ ) −1 − (k · k ′ )(p a · ε · p a ) (p a · ε ′ · p a ) + 2(p a · k ′ )(p a · ε · p a )(p a · ε ′ · k) + 2(p a · k)(p a · ε ′ · p a )(p a · ε · k ′ ) − 2(p a · k)(p a · k ′ )(p a · ε · ε ′ · p a ) + (k · k ′ ) −1 − (k ′ · ε · ε ′ · p a )(k ′ · p a ) − (k · ε ′ · ε · p a )(k · p a ) + (k ′ · ε · ε ′ · p a )(k · p a ) + (k · ε ′ · ε · p a )(k ′ · p a ) − ε γδ ε 2γδ (k · p a )(k ′ · p a ) − 2(p a · ε · k ′ )(p a · ε ′ · k) + (p a · ε ′ · p a )(k ′ · ε · k ′ ) + (p a · ε · p a )(k · ε ′ · k) , (1.4) and J µν is the spin angular momentum which acts on the polarization vectors ε 1 , · · · ε n in M n . J µν is normalized such that acting on a covariant vector ε ρ , it gives Our goal in this paper will be to test this formula using the Cachazo-He-Yuan (CHY) formula [66][67][68][69][70][71][72] for tree level graviton scattering amplitude in Einstein's theory in any 1 In [53] the amplitudes Mn included the momentum conserving delta-functions in their definition. In this paper we shall use the CHY formula for Mn and will not include the momentum conserving delta functions in the definition of Mn. As has been discussed below (3.5), the soft graviton theorem takes the same form as in [53] for amplitudes without delta functions.

JHEP01(2018)090
dimensions. Our final formula, given in eq. (5.50), turns out to be in perfect agreement with (1.4). This provides a non-trivial test of the multiple soft graviton theorem derived in [53] and also of the CHY formula for graviton scattering amplitudes in Einstein theory.
Refs. [50,51] derived the double soft graviton theorem from CHY scattering equation. 2 In [53] we noticed that this differs from our general formula by a sign. In this paper we show by careful analysis that the sign that we get from the double soft limit of the CHY scattering formula agrees with that of [53] provided we normalize the amplitude so that the single soft theorem comes with the correct sign. We also find a few extra terms in the intermediate stages of analysis that were missed in the analysis of [50,51], but they cancel in the final result.
The rest of the paper is organized as follows. In section 2 we review the CHY formula for scattering amplitude of n gravitons. In section 3 we study the limit of this formula in which one graviton becomes soft, and derive the single soft graviton theorem. In this section we also fix the normalization of the amplitudes to ensure that the single soft graviton theorem comes with the conventional sign. In section 4 we study the limit of the CHY formula when two gravitons become soft and show that the result is consistent with the general formula found in [53] including the sign. One distinguishing feature of our analysis is that we work with general choice of gauge for various polarization tensors since the gauge choice that simplified previous analyses is not available for multiple soft gravitons. 3 Finally in section 5 we analyze CHY formula in the limit in which multiple gravitons become soft, and show that the result agrees with the general formula derived in [53]. Appendix A describes the counting of different types of solutions of CHY scattering equations for multiple soft gravitons.

Cachazo-He-Yuan prescription for amplitude calculation
For computation of tree level scattering amplitude of massless particles in arbitrary dimension, Cachazo-He-Yuan gave a prescription known as CHY prescription. Our goal in this section will be to describe this formula. However since this makes use of delta functions with holomorphic arguments, we shall first set up appropriate conventions to deal with these delta functions.

Conventions
In writing down the CHY formula for scattering amplitudes, one makes use of delta functions for holomorphic variables defined via the equation dσ δ(f (σ)) F (σ) = (α) (f ′ (σ (α) )) −1 F (σ (α) ) (2.1) 2 Double and multiple soft theorem for other soft particles from CHY scattering equations has been studied in [45,72,73]. 3 The only exception is that we do impose ε µ µ = 0 gauge condition on the polarization tensor ε of the soft gravitons. Using on-shell condition this also implies k µ εµν = 0 where k is the momentum of the soft graviton.
The absence of absolute value in the right hand side of (2.1) can some time cause confusion when we have multiple integration. Consider for example the integral dx dy δ(x − y) δ(x + y) F (x, y) . On the other hand if we use the first delta function to carry out the y integral first then the result is Therefore the two results do not agree. This ambiguity can be resolved if we regard the delta functions as grassmann odd objects so that exchanging their positions costs a sign. We also regard the integration measure as wedge products so that changing the order of doing the integration also costs a sign. For a given order of integration and delta functions we shall follow the convention that the last integration will be done first using the last delta function, and successive integrations will follow the same order. If we want to use a different delta function for the last integration, we need to first bring that particular delta function to the end picking up appropriate minus sign and then carry out the integration using (2.1) or (2.2).
It is easy to see that with this convention, under a change of variables the product of delta functions pick up the inverse of the Jacobian without absolute value where on both sides the delta functions in the product are arranged from left to right in the order of increasing i.

CHY formula
We are now ready to describe the CHY proposal [67][68][69][70][71][72]. According to this proposal, an n-point tree amplitude of massless particles can be derived from a sum over discrete set of points in the moduli space of an n-punctured Riemann sphere M 0,n . The position of JHEP01(2018)090 the punctures corresponding to these points in M 0,n are determined from the solutions of scattering equations: where {σ a } are the holomorphic coordinates of the punctures and {p a } are momenta of massless particles. Using the SL(2, C) invariance of M 0,n we can fix the positions of three punctures and use any (n − 3) of the n scattering equations as constraint relations. The CHY formula is given as, Here we used SL(2, C) invariance to fix three punctures p, q, r, and made use of the linear dependence of the scattering equations to remove three scattering equations for i, j, k particles. The expression (2.8) is independent of the choice of p, q, r and i, j, k if Using this requirement, together with multi-linearity in all the polarization vectors of ngravitons and gauge invariance of M n , CHY gave the following form of the integrand I n for Einstein gravity [67]: where Ψ is a 2n × 2n anti-symmetric matrix defined below and Ψ st st is obtained by removing s-th and t-th row from first n rows and removing s-th and t-th columns from first n columns of Ψ. The (−1) n factor was not present in the original formula, but we have included it in order to get soft graviton theorems with conventional signs for the leading soft graviton theorem. This factor has also appeared recently in [74]. The overall n independent normalization and sign of the amplitude will be irrelevant for our analysis. Ψ has the form: where A, B, C are n × n matrices defined as,

JHEP01(2018)090
Here for each a, ǫ a is a space-time vector. In the expression for I n , each ǫ a appears quadratically. In the final formula, we replace ǫ a,µ ǫ a,ν by ε a,µν and identify ε a as the polarization tensor for the a-th graviton.
Since the total number of independent scattering equations is equal to the total number of independent variables, the integration over σ a 's reduce to a sum over discrete set of points, it can be seen using simple counting that the total number of solutions to the scattering equation is (n − 3)! [66].
In the rest of the paper we shall use the CHY formula (2.8) to study soft limitslimits in which one or more external momenta become small. Our general strategy will be as follows. We first represent the integration over σ a 's as contour integrals using (2.2), taking the contours to lie close to the solutions of the scattering equations. In this region we can make appropriate approximation on the integrand using the soft limit. Once we make this approximation, we now deform the contours away from the solutions of the scattering equation, possibly through regions where the original approximation fails, but Cauchy's theorem guarantees that the value of the integral is unchanged. As we shall see, this gives an effective method for approximating the CHY formula in the soft limit.

Single soft-graviton theorem
In this section we shall review the derivation of subleading single soft graviton theorem from CHY prescription given in [7,12,16,17], taking special care about the signs. We shall see that the (−1) n factor in (2.10) is necessary for getting the conventional signs. We also work with general polarizations of external gravitons without making any gauge choice (other than the one given in footnote 3) as in [17] -this will be necessary for generalizing the analysis to multiple soft graviton theorem where special gauge choices of the type used e.g. in [12] is not possible. Single soft graviton theorem for general gauge choice has been derived in [17].

Statement of the theorem
We shall begin by stating the subleading single soft graviton theorem in the standard form so that we know what we want to prove. We consider the scattering amplitude of (n + 1) massless gravitons with (n+1)'th particle momentum soft and expand the amplitude about the soft momentum. Let p 1 , · · · p n denote the momenta of the finite energy particles, and p n+1 = τ k n+1 denote the momentum of the soft particle. We shall take the soft limit by taking τ → 0 limit at fixed k n+1 . We also denote by ε i the polarization tensor of the i-th graviton. Then the subleading soft theorem stated in (1.1) takes the form where, Here, J µν a is the total angular momentum acting on the a-th state, defined as, where J µν a is the spin angular momentum of the a th hard graviton, defined as There is one subtle point that must be mentioned here. The full amplitude contains a momentum conserving delta function, and the derivation of the multiple soft graviton theorem in [53] included this momentum conserving delta function in the definition of the amplitude. In that case the external momenta can be taken as independent and the derivatives with respect to momenta, present in the definition of J µν , are well defined. On the other hand, the CHY formula given in (2.8) does not have this delta function. It is however easy to see that as long as (2.8) satisfies the soft graviton theorem, the full amplitude including the delta function also satisfies the soft graviton theorem [6]. For this we multiply both sides of (3.1) by δ (D) (p 1 + · · · + p n + τ k n+1 ), and manipulate the right hand side as follows: where we have used the fact that n a=1 ε (n+1),µν p µ a k (n+1),ρ k n+1 .p a p ρ a ∂ ∂p a,ν δ (D) (p 1 + · · · + p n ) = 0 , as a consequence of the condition ε µ µ = 0 that we impose on the polarization tensors of the soft gravitons. Since the right hand side of (3.6) now takes the form of the right hand side of the soft graviton theorem with the momentum conserving delta function included in the definition of the amplitude, we see that once we prove eq. (3.1) for the CHY amplitudes, it also holds for the full amplitude including the momentum conserving delta functions. A similar result holds also for the multiple soft graviton theorem given in (1.1).
While proving (3.1), or more generally (1.1) for the CHY amplitude (2.8), we shall treat all the external momenta as independent while computing the derivatives of the amplitude with respect to the external momenta. However while verifying the equality of the two sides of the equation we shall make use of the conservation of total momenta, since the final identity we are interested in proving is not (3.1) or (1.1), but the ones obtained from these by multiplying both sides of the equation by momentum conserving delta functions.

Single soft limit of CHY formula
We now begin the analysis of the single soft limit of CHY formula. Let us introduce some compact notations, 4 and rewrite the CHY formula for the (n + 1)-point scattering amplitude as, While carrying out various manipulations below, we shall not be careful about orders of the delta functions. However it should be understood that while doing the final integration over the σ a 's, the delta functions must be brought to the same order in which they were arranged initially. With this understanding we expand the product of delta functions inside the integral (3.10) in a Taylor series in τ as: The integrand can also be expanded in powers of τ as, Now substituting (3.11) and (3.12) in (3.10) and keeping terms up to subleading order in τ , we obtain Now we shall analyze the three terms appearing in equation (3.13) one by one. 4 For single soft gravitons we only need f n n+1 , but later we shall make use of the definition of f n n+r for general r.

First term
Using the expressions (2.10), (2.11) and (2.12), we get, on putting τ = 0, (3.14) Here, the negative sign comes from the (−1) n factor in the definition of I n . Using this, the first term on the right hand side of (3.13) can be written as (3. 15) In going to the last line, we have used the contour integral representation of the delta function as described in (2.2), with the σ n+1 integration contour wrapping the solutions of the scattering equation of the (n + 1)-th particle. We now deform the contour and do the contour integration about the poles of rest of the integrand. The deformed contour will wrap the other poles in clockwise direction, and therefore the residue theorem will have an extra minus sign. We shall also need to take care of poles at infinity if they exist. If we denote collectively by {R i } all these poles, the first term becomes (3.16) Using momentum conservation n a=1 p a = −τ k n+1 and the on-shell conditions k 2 n+1 = 0, ǫ n+1 .k n+1 = 0, we can see that for large σ n+1 , the last term in (3.16) falls off as (σ n+1 ) −4 and the last but one term grows as (σ n+1 ) 2 . Therefore there is no pole at ∞. The rest of the poles inside the deformed contour are the simple poles at σ n+1 = σ a coming from the combination of the last two terms in (3.16). Computing the residue at each such pole, we obtain the result Replacing ǫ n+1,µ ǫ n+1,ν by ε n+1,µν in the above expression, we see that it agrees with the leading soft theorem (3.1) including the sign.

Second term
The second term on the right hand side of (3.13) is given by (3.18) Let us focus on the σ n+1 integration. The relevant integral is (3. 19) We can analyze this by deforming the σ n+1 integration contour in the same way as before.
The only extra complication is that along with the simple poles at σ n+1 = σ a (a = ℓ), we now also have a double pole at σ ℓ and we have to be a little more careful in evaluating the residue. This can be done, yielding the result Substituting this into (3.18) we get The derivative of the delta function can be obtained by using the representation (2.2) of the delta function. Instead of contour integral of 1/f we have contour integral of −1/f 2 . We shall see, however, that we do not need to make use of this explicit representation.

Third term
We now turn to the evaluation of the third term on the right hand side of (3.13) We begin by introducing some notations. We denote by Ψ the matrix Ψ st st defined in (2.11) and the paragraph above it for n finite energy external states, with some fixed choice of s and t in the range 1 ≤ s, t ≤ n. We shall denote by Ψ the same matrix for n + 1 external

JHEP01(2018)090
states with the (n + 1)-th state describing a soft graviton carrying momentum τ k n+1 . Ψ is a (2n − 2) × (2n − 2) matrix and Ψ is a 2n × 2n matrix. We shall denote by P the Pfaffian of Ψ and by P the Pfaffian of Ψ: (3.23) Furthermore we shall denote by P αβ = − P βα for α < β the Pfaffian of the matrix obtained by eliminating from P the α-th row and column and β-th row and column. Similarly P αβγδ is defined to be totally anti-symmetric in α, β, γ, δ, which, for α < β < γ < δ, is given by the Pfaffian of the matrix obtained by eliminating the α, β, γ, δ-th rows and columns of the matrix Ψ. Similar definitions and properties hold for P αβ and P αβγδ . Then we have and therefore Now from the definition (3.23) of the Pfaffian it follows that Here we have adopted the notation that in Ψ αβ and P αβ we shall let the indices run over 1 ≤ α, β ≤ (n + 1), α, β = s, t, and (n + 2) ≤ α, β ≤ 2n + 2. Therefore the last rows and columns will be called (2n + 2)-th rows and columns, even though they are actually 2n-th rows and columns. We have also implicitly assumed that the integers s and t are consecutive -otherwise there will be extra minus signs when one of α or β falls between s and t. Our final result will be valid even when s and t are not consecutive. Now using the explicit form of the (n + 1) × (n + 1) matrices A, B, C given in (2.12) for n + 1 gravitons, with p n+1 = τ k n+1 , we get,

JHEP01(2018)090
On the other hand using (3.26) and (2.11), we get Using (3.27), (3.28) we can rewrite (3.25) as, We shall now evaluate the various components P αβ that appear in the expression (3.29). For this it will be useful to express Ψ in the matrix form: where the indices a, b, c, d for the matrices A, B, C run from 1 to n (with a, b = s, t), and the (n + 1)-th components of these matrices have been written down explicitly. Also useful will be the explicit form of C (n+1)(n+1) : We begin with P a(n+1) . This is the Pfaffian of the matrix (3.30) without the a-th and (n + 1)-th rows and columns. Expanding this about the last column we get (3.32) We now notice that for τ = 0, the matrix obtained by eliminating the (n + 1)-th and (2n + 2)-th rows and columns of Ψ, is in fact the matrix Ψ. Therefore we can rewrite (3.32) as

JHEP01(2018)090
Similarly we have (3.34) To analyze P a(n+1+a) we expand it in the (n + 1)-th row. For τ = 0 only the (2n + 2)-th column contributes and gives In order to evaluate the right hand side of (3.29) we also need P | τ =0 . This is evaluated by expanding the Pfaffian of the matrix Ψ given in (3.30) in the (n + 1)-th row. The only contribution comes from the last element in the row, giving the result Substituting these results in (3.29) and using (3.27) we get the result for I (1) n+1 which can be substituted into (3.22). Using the definition of f n n+1 given in (3.8) and the definition (2.2) of the delta function, we now get The σ n+1 contour winds anti-clockwise around the zeroes of We can now deform the contour towards infinity, picking up residues from σ n+1 = σ a , -it is easy to see that there are no poles at ∞. For the coefficients of P ab , P (n+a)(n+b) and P a(n+b) for a = b, the integrand has single poles at σ n+1 = σ a and σ n+1 = σ b , and the result can be computed easily by picking up the residues. Finally for the coefficient of P a(n+a) the integrand has a double pole at σ n+1 = σ a , therefore to evaluate the residue we have to expand it in a Laurent series around the pole and pick the coefficient of the (σ n+1 − σ a ) −1

Comparison with soft graviton theorem
The first term A 1 given in (3.17) exactly matches the leading soft graviton result given by the first term on the right hand side of (3.1) once we make the identification ε a,µν = ǫ a,µ ǫ a,ν .
To check that the second and the third terms A 2 and A 3 given respectively in (3.21) and (3.38) agree with the known subleading soft graviton theorem, we first express (3.3) as Therefore we have, from (2.8)

JHEP01(2018)090
First let us compute the operation of S (1) orbital on the delta function, With this, the first term on the right hand side of equation ( This exactly agrees with A 2 given in equation (3.21).
It remains to show that the second term on the right hand side of equation (3.41) agrees with A 3 . Using the definition of I n given in (2.10) and that det(Ψ st st ) = det Ψ = P 2 , we get (3.44) In order to calculate the action of S (1) on this, it will be useful to record how S spin defined in (3.40) acts on the 'square root' ǫ of the graviton polarization tensor ε. This is done by identifying the action of J µν on ǫ as This allows us to express S and therefore

JHEP01(2018)090
Using (3.46) and (2.12) we get We now use a formula similar to (3.28) with (n + 1) replaced by n and ∂ τ replaced by S and use (3.47), (3.48) to express the second term on the right hand side of eq. (3.41) as This is identical to (3.38). This finishes the proof that the CHY prescription for graviton scattering amplitude is consistent with the subleading single soft graviton theorem.

Double soft graviton theorem
In this section we shall consider the limit of the graviton scattering amplitude when two of the gravitons carry soft momenta. We shall reproduce the result of [50,51] with opposite JHEP01(2018)090 sign and in general gauge. During this analysis we also find some additional terms that were left out in the analysis of [50,51], which nevertheless cancel at the end.

Double soft limit
We shall first follow [45,72] to analyze the solutions to the scattering equations of (n + 2) particles with momenta {p 1 , p 2 , . . . , p n , τ k n+1 , τ k n+2 } in the double soft limit τ → 0 keeping k n+1 , k n+2 fixed. The scattering equations of first n particles, (n + 1)-th particle and (n + 2)-th particle are given respectively by: where we have removed overall factors of τ from the last two equations. We also have momentum conservation relation, The solution to the scattering equations can be divided into two classes [45,72]: degenerate solutions where σ n+1 − σ n+2 ∼ τ and non-degenerate solutions for which σ n+1 − σ n+2 ∼ 1 in the τ → 0 limit. For both classes of solutions σ a 's for 1 ≤ a ≤ n remain finite distance away from each other and from σ n+1 , σ n+2 as τ → 0. In order to show that there are no other types of solutions we shall now show that the total number of degenerate and nondegenerate solutions add up to (n − 1)! -the actual number of solutions for scattering equation of (n + 2) particles [67]. At τ = 0, eq. (4.1) describes the scattering equations for n particles which have (n−3)! solutions for {σ 1 , σ 2 , . . . , σ n }. Generically these solutions are non-degenerate, with noncoincident σ a 's. Therefore to count the number of non-degenerate and degenerate solutions in the τ → 0 limit, we have to multiply (n − 3)! by the number of solutions for σ n+1 and σ n+2 for fixed σ 1 , · · · , σ n .
Let us first count the number of non-degenerate solutions. For fixed σ 1 , · · · , σ n , we can ignore the last term in (4.2) in τ → 0 limit and express this as a polynomial equation for σ n+1 : Naively this is of degree (n − 1), but momentum conservation (4.4) makes the coefficient of the (σ n+1 ) n−1 term −τ k n+1 .k n+2 which vanish in the τ → 0 limit. Therefore in this limit this is a polynomial equation of degree (n−2) and gives (n−2) solutions for σ n+1 . Similarly

Contribution from non-degenerate solutions
The contribution from the non-degenerate solutions to subleading order in τ may be written as,

JHEP01(2018)090
where it is understood that we pick contributions from those zeroes of the δ-functions for which σ n+1 − σ n+2 ∼ 1 as τ → 0. The product of first (n − 3) delta functions can be expanded in powers of τ as: (4.14) On the other hand the product of last two delta functions have the form: where f n n+r has been defined in (3.8). Finally the integrand I n+2 may be expanded as, Substituting all these expansions in the (n + 2)-point amplitude (4.13), we get For evaluating this, it will be convenient to introduce two soft parameters τ 1 and τ 2 instead of a single parameter τ and take the external momenta to be p 1 , · · · p n , τ 1 k n+1 , τ 2 k n+2 . I n+2 is then given by

JHEP01(2018)090
where P is the Pfaffian of the matrix In (4.20) we have set terms with two or more powers of τ to be zero. The values of C (n+1)(n+1) and C (n+2)(n+2) to linear order in τ are: At the end of the computation we shall take the limit τ 1 , τ 2 → τ → 0. We shall now analyze the different terms on the right hand side of (4.18). However, before we proceed to evaluate the integrals, we would like to remind the reader of the general strategy for analyzing these integrals that was described in the last paragraph of section 2. The expansion (4.18) of the amplitude in powers of τ is valid when the integration contours wrap around the non-degenerate solutions for which σ n+1 − σ n+2 ∼ 1. Once we have made this approximation we can deform the contours of σ n+1 and σ n+2 even to regions where |σ n+1 − σ n+2 | ≪ 1. Treating (4.18) as the integrand still gives the correct result by Cauchy's theorem even though the original integrand is no longer approximated by (4.18). A similar remark will hold for the contribution from degenerate solutions -we shall make our approximation in the region where the integration contour is close to the degenerate solution of the scattering equation, and once the approximation is made, we shall be free to deform the contour away from the region where the approximation is valid.
First two terms. In (4.18), the first two terms on the right hand side may be analyzed as in the case of single soft graviton by carrying out the integration over σ n+2 and σ n+1 independently. For this we note from (4.20) that 22) where P is the Pfaffian of the matrix Ψ for n-particle scattering amplitude without soft graviton. This gives (4.23)

JHEP01(2018)090
We can now proceed to evaluate the first two terms on the right hand side of (4.18) as in sections 3.2.1 and 3.2.2. The first term is given by where {A i } and {B i } are defined as the set of points satisfying We can now carry out integration over σ n+1 and σ n+2 independently by deforming the contours to ∞. The contribution comes from residues at the poles at σ a -which can be evaluated following the procedure described in section 3.2.1 -and the contribution from ∞ can be evaluated after using momentum conservation. The result is Using M n = Dσ δ (0) I n , this may be rewritten as where we define For evaluating the second term on the right hand side of (4.18) we note that δ (1) defined in (4.14) contains sum of two sets of terms -one set involving 1/σ l(n+1) and the other set involving 1/σ l(n+2) . First consider the set involving 1/σ l(n+1) . In this case we can carry out the integration over σ n+2 first following the procedure of section 3.2.1 and arrive at the product of S (0) n+2 times an integral of the form given in (3.18). The contribution from the pole at σ n+2 = ∞ can be ignored since that will produce an extra factor of τ and give a subsubleading contribution. The integration over σ n+1 is proportional to (3.18) and can be analyzed as in section 3.2.2, leading to the result (3.21) multiplied by S (0) n+2 . Using the JHEP01(2018)090 equality between (3.21), (3.43) and the first term on the right hand side of (3.41), this may be expressed as (4.29) where in the left hand side of this equation it is understood that only the orbital part of S (1) n+1 acts on δ (0) . The contribution from the second term in δ (1) introduced in (4.14) is given by an expression similar to that in (4.29) but with (n + 1) and (n + 2) interchanged. Therefore the total contribution from the second term on the right hand side of (4.18) may be expressed as (4.30) We shall now turn to the analysis of the last three terms in (4.18).
Third term. We now consider the third term on the right hand side of (4.18): To evaluate this we have to first evaluate I n+2 . We have the analog of (3.25): 32) where P is the Pfaffian of the matrix Ψ given in (4.20).
Let us now examine the τ 1 derivative of P . Since we have to set τ 2 = 0 at the end, we can make this replacement even before taking the τ 1 derivative. In this case the Pfaffian P of the matrix (4.20) reduces to One can now recognize this matrix as the same matrix that appears in (3.30). Therefore when the τ 1 derivative acts on the Pfaffian of this matrix, then at τ 1 = 0 it will have the same form as that in (3.38) after σ n+1 integration. The extra factor of C (n+2)(n+2) at τ 1 = 0, multiplied by a similar factor coming from the P factor in (4.32), will generate a factor of S (0) n+2 after integration over σ n+2 . There is an additional contribution from σ n+2 = ∞, but this is subsubleading and may be ignored. Using the equality of (3.38), (3.50) and second term on the right hand side of (3.41), the net contribution of this term to the right hand side of (4.18) can be shown to have the form: (4.34)

JHEP01(2018)090
A similar term with (n + 1) and (n + 2) exchanged comes from the ∂ P /∂τ 2 term in (4.32). The sum of these two give This cancels the last two terms in (4.30). This however is not the complete contribution from the third term in (4.18), since we still have to include the contribution where τ 1 derivative acts on the C (n+2)(n+2) factor in (4.33), and a similar contribution with (n + 1) and (n + 2) exchanged. Now from (4.21) we have Therefore at τ 1 = τ 2 = 0, the contribution to ∂ P /∂τ 1 , with the derivative acting on the first term in (4.33), reduces to where P is the Pfaffian of the matrix A ab −(C T ) ad C cb B cd . Substituting this into (4.32) and setting τ 1 = τ 2 = 0, using (4.22), and adding also the extra contribution from the ∂ P /∂τ 2 term in (4.32), we get the net extra contribution to I (1) n+2 to be Noting that 4 (−1) n (σ s − σ t ) −2 P 2 = I n , the contribution of (4.38) to (4.31) is given by We now regard the integrals over σ n+1 and σ n+2 as contour integrals as in (2.2). Let us examine the first term in the square bracket. It takes the form

JHEP01(2018)090
where, {A i } and {B i } have been defined in (4.25) and accordingly the σ n+1 and σ n+2 contours run anti-clockwise around the poles of the first and second factors of the integrand respectively. We shall now deform the σ n+2 contour away from the pole towards infinity. The only pole is at σ n+2 = σ n+1 and a possible contribution from σ n+2 = ∞. Let us first examine the contribution from the pole at σ n+1 . This is given by Next we analyze the contribution from σ n+2 = ∞. Using momentum conservation equation n a=1 p a = −τ (k n+1 + k n+2 ) this is given by (4.42) Now we deform the σ n+1 contour to ∞, picking up residues at σ a . The result is The residue at σ n+1 = ∞ has an additional factor of τ . Therefore its contribution is subsubleading. The contribution from the second term in the square bracket in (4.39) can be obtained from (4.41), (4.43) by exchange of (n + 1) and (n + 2) and is given by B 5 + B 6 , where and Fourth and fifth terms. The fourth term on the right hand side of (4.18) is given by: (4.46)

JHEP01(2018)090
We now deform the σ n+1 contour away from the poles {A i } towards infinity and pick up residues at σ n+1 = σ n+2 as well as the contribution from infinity. There are no poles at σ n+1 = σ a . The contribution from the pole at σ n+2 is given by (4.47) On the other hand the contribution from infinity is given by, after using momentum conservation, where in the second step we have performed integration over σ n+2 , picking up residues at σ b . In this case, as σ n+2 → ∞, the contribution to the integrand goes as τ , and therefore this does not contribute at the subleading order. The contribution from the fifth term on the right hand side of (4.18) can be evaluated in a similar manner, giving contributions similar to (4.47) and (4.48), with n + 1 and n + 2 exchanged: (4.49) and Adding the contributions B 0 to B 10 we get the total contribution from non-degenerate solutions: (4.51) -25 -

JHEP01(2018)090
The term in the first line was derived in [51]. The terms in the second and the third line vanish in the gauge chosen in [51]. The terms in the fourth and fifth lines were missed in the analysis of [51]. We shall see that these terms cancel similar terms arising in the analysis of the degenerate solutions. Using the gauge transformation laws of S n+1 M n analyzed in [53], one can verify that (4.51) is invariant under the gauge transformations of the soft graviton polarizations.

Contribution from degenerate solutions
When the integration contours of σ n+1 and σ n+2 wrap around the degenerate solutions for which |σ n+1 − σ n+2 | ∼ τ , we will carry out the integration in ρ and ξ variables introduced in (4.7). In terms of these variables the integration measure takes the form, with the understanding that ρ integration is done using the first delta function and ξ integration is done using the second delta function. Therefore the contribution to the (n + 2)-point amplitude from the degenerate solutions becomes, where it is understood that we sum over contributions from those zeroes of the second δ-function for which ξ ∼ τ . By carrying out the integration over ξ using the second delta function, this integral can be approximated as where (4.55)

JHEP01(2018)090
Now in the degeneration limit we can evaluate I n+2 by regarding this as the square of the Pfaffian of the matrix Ψ which now takes the form where we now have (4.57) Given the appearance of 1/τ in the (2n + 3)-(2n + 4)-th element of the matrix, one may wonder whether the Pfaffian may give a leading contribution of order 1/τ , thereby upsetting the counting of powers of τ . It is easy to see however that the coefficient of this element in the expansion of the Pfaffian, given by the Pfaffian obtained by eliminating the (2n + 3) and (2n + 4)-th rows and columns of this matrix, actually goes as τ and therefore does not upset the counting of powers of τ . With this understanding we can now compute the Pfaffian of Ψ, arriving at the result [72]: Substituting this into (4.53) we get (4.59) -27 -

JHEP01(2018)090
By making use of the δ-function, we can rewrite this as Note that there are many other ways of writing this expression -by replacing one or more factors of n b=1 We have chosen the one that will be most convenient for our analysis, but the final result does not depend on this choice.
We now represent the delta function as a contour integration as in (2.2) (4.61) and deform the ρ integration contour away from the zeroes of the argument of the delta function. This will generate three kinds of terms, the residues at ρ = σ a for 1 ≤ a ≤ n, residues at the zeroes of n a=1 k n+1 .p a ρ − σ a and n a=1 k n+2 .p a ρ − σ a -called {A i } and {B i } in (4.25) -and residue at ∞. We shall now analyze these three kinds of terms one by one.
Residue at σ a . This gives

JHEP01(2018)090
where the overall minus sign comes from the reversal of the integration contour. Expanding the square and using the convention we get where M(p a ; ε n+1 , k n+1 , ε n+2 , k n+2 ) (4.65) Residue at {A i } and {B i }. The contribution to (4.60) from the residues at A i and B i defined in (4.25) are given by (4.66) These terms were missed in the analysis of [50]. They cancel similar contributions (4.51) coming from non-degenerate solutions.

JHEP01(2018)090
Residue at ∞. Finally we can examine the residue of the integral (4.60) at ρ = ∞. This is given by (4.67) Using momentum conservation law n a=1 p a = −τ (k n+1 + k n+2 ) we can see that this is of order τ . Therefore this contribution is subsubleading and can be ignored.

Total contribution
The full amplitude corresponding to two soft gravitons is obtained by adding the contributions (4.51) from the non-degenerate solutions and the contributions (4.64) and (4.66) from the degenerate solutions. This is given by This agrees with the result of [53] reviewed in (1.1) for two soft gravitons. In [53] it was also shown that (4.68) is invariant under the gauge transformations of the soft graviton polarizations. Note that (4.66) cancels part of the contribution from (4.51). This suggests that there may be a better way of organising the calculation instead of representing it as a sum of contributions from degenerate and non-degenerate solutions.

Multiple soft graviton theorem
In this section we shall generalize the analysis of the previous section to the case where arbitrary number of gravitons become soft.

Degenerate and non-degenerate solutions
We assume that there are n + m number of gravitons and m of them become soft. We parametrize the m soft momenta as p µ a = τ k µ a , a = n + 1, · · · , n + m , (5.1) and take the soft limit by taking τ → 0 at fixed k a . The momentum conservation now takes the form p µ 1 + · · · + p µ n + τ (k µ n+1 + · · · + k µ n+m ) = 0 .
In this case the full solution space of the scattering equations gets divided into different sectors. These different sectors correspond to the case when a group of r 1 punctures associated with soft gravitons come within a distance of order τ of each other, another group of r 2 punctures associated with soft gravitons come within a distance of order τ of each other and so on. A detailed analysis of the number of solutions of each type can be found in appendix A. In this section our goal will be to prove that for the subleading multiple soft graviton amplitude, only two sectors contribute -non-degenerate solutions where all the punctures are finite distance away from each other and degenerate solutions where two of the punctures come within a distance τ of each other and all other punctures are finite distance away from each other.
To prove this, we note that the CHY formula for the amplitude is given by We now analyze the product of delta functions and focus on the last m of them corresponding to the soft gravitons. These m delta functions are given by Each of these delta functions gives a factor of 1 τ to the amplitude irrespective of how many of the m σ a 's come within a distance of order τ of each other. These are the only source of singular τ dependence coming from the delta functions. Therefore we get a net factor of τ −m from the delta functions.
Next, we consider the τ factors coming from the measure when r of the m σ a 's associated with the soft gravitons come within a distance of order τ of each other. Without loss of generality we can label these as σ n+m−r+q for 1 ≤ q ≤ r. We now make following redefinitions of the coordinates σ a = σ ′ a for 1 ≤ a ≤ n + m − r + 1 , σ n+m−r+q = σ ′ n+m−r+1 + τ ξ q , (q = 2, · · · , r) .

JHEP01(2018)090
We shall prove shortly that I n+m does not give rise to any singular behavior in the τ → 0 limit for finite σ ′ c , ξ q . Assuming this to be the case, we see that this contribution is leading only for r = 1 and can receive subleading contribution only for r = 1, 2. In other words, for the subleading soft graviton amplitude, the contributions come only from those solutions for which none of the punctures go close to each other (non degenerate solutions) or two of the punctures go close to each other (degenerate solutions).
Let us now prove that I n+m has a finite limit as τ → 0 with σ ′ c , ξ q fixed (i.e., when r punctures go close to each other). We define and express I n+m as where P is the Pfaffian of the matrix  · · · C p+r,p+r where the matrices A, B and C are defined as in (2.12) for m+n = p+r particles and A p , B p and C p denote the first p × p blocks of these matrices. Using the parametrization (5.5), the matrix (5.9) may be expressed as 

JHEP01(2018)090
We now note the following features of the matrix shown in (5.10): 1. (A p ) ab , (B p ) ab and (C p ) ab have finite τ → 0 limit.
2. There are four blocks of the matrix formed by the vertical double line and the horizontal double line. We note that (a) The upper left block has r rows and r columns that are proportional to τ .
(b) The upper right and lower left blocks have finite τ → 0 limit.
(c) In the lower right block, the diagonal matrix elements vanish, whereas the offdiagonal matrix elements are proportional to 1/τ .
We now note that the inverse powers of τ appear only in the elements of the lower right block of size r × r. Let us suppose first that r is even. Then we get a maximum contribution of τ −r/2 from the Pfaffian of the lower right block, and the coefficient of this term is proportional to the Pfaffian of the upper left block of (5.10). This matrix has r rows and columns with every element proportional to τ , and therefore its Pfaffian will be proportional to τ r/2 , cancelling the τ −r/2 factor. Therefore this term is finite as τ → 0.
If r is odd, then the singular terms in the lower right block can give a maximum contribution of τ −(r−1)/2 . This is given by a sum of terms, one of which is proportional to the Pfaffian of the matrix given by the lower right block of size (r − 1) × (r − 1), multiplied by the Pfaffian of the matrix obtained by eliminating the last (r − 1) rows and columns of Ψ. The other terms are related to this one by rearrangement of the last r rows and columns and can be analyzed similarly. It is easy to see that the matrix obtained by eliminating the last (r − 1) rows and columns of (5.10) has r rows and columns proportional to each other in the τ → 0 limit, and therefore its Pfaffian gives a factor of τ (r−1)/2 . This cancels the τ −(r−1)/2 factor, giving a finite τ → 0 limit.
Next we need to consider the possibility that we may not choose the maximally singular terms from the lower right block. One such term corresponds to choosing the Pfaffian of the matrix given by the last k × k block for some even integer k < r, multiplied by the Pfaffian of the matrix obtained by eliminating the last k rows and columns, but from the latter we do not pick any term that has 1/τ factor. The Pfaffian of the k × k matrix goes as τ −k/2 , whereas the matrix obtained by eliminating the last k rows and columns has r rows (and columns) given by linear combinations of (r − k) independent vectors in the τ → 0 limit. Therefore its Pfaffian goes as τ k/2 . This again shows that the product has a finite τ → 0 limit. The other terms of this kind are related to the one discussed above by rearrangement of the last r rows and columns and are therefore also finite as τ → 0.
This finishes our proof that I n+m remains finite in the τ → 0 limit.

Contribution from non-degenerate solutions
In this section we shall compute the contribution to the amplitude from non-degenerate solutions to the scattering equations. The amplitude is given by M n+m (p 1 , . . . , p n , τ k n+1 , · · · , τ k n+m ) = Dσ where f n n+r has been defined in (3.8). Finally, expansion of the integrand I n+m gives, n+m + O(τ 2 ) . (5.14)

JHEP01(2018)090
Then, up to subleading order, the non degenerate contribution becomes n+m .

JHEP01(2018)090
With this, (5.17) gives where P is the Pfaffian of (5.17) for m = 0. Therefore we have (5.20) To evaluate I (1) n+m we use the analogue of (4.32): ∂ P ∂τ r τ 1 =τ 2 =...=τm=0 . (5.21) The strategy for evaluating P ∂ P ∂τ r {τq=0} will be to first set τ q = 0 for q = r in the matrix (5.17) and then expand the Pfaffian successively about the rows (n + 1) to (n + m) except the (n + r)-th row. This gives where P (n) r is the Pfaffian of the matrix Ψ for (n + 1) graviton scattering with the first n gravitons carrying momenta p 1 , · · · , p n and polarizations ε 1 , · · · ε n and the last one carrying soft momentum τ r k n+r and polarization ε n+r . The overall sign in (5.22) will not be needed for our analysis. Using We shall now proceed to evaluate the right hand side of (5.15).
First term. The first term on the right hand side of (5.15) may be expressed as where for fixed r, {A r,i } are defined as the set of points satisfying We can now carry out integration over all {σ n+r } independently by deforming the contours to ∞. The contribution comes from residues at the poles at σ a -which can be evaluated following the procedure described in section 3.2.1 -and the contribution from ∞ can be evaluated after using momentum conservation. The result up to subleading order is Second term. The second term on the right hand side of (5.15) is given by (5.28)

JHEP01(2018)090
The integrals over {σ n+q } for q = v are of the form (3.15) and can be analyzed as in section 3.2.1 to produce factors of −S (0) n+q . The contribution from the pole at σ n+q = ∞ is subsubleading and can be ignored. The remaining integration over σ n+v has the form given in (3.18) and can be analyzed as in section 3.2.2. The final result, given in (3.21), can in turn be related, using the equality of (3.21), (3.43) and the first term on the right hand side of (3.41), to where S n+v has been defined in (4.28). Then the result up to subleading order is (5.30) Third term. We shall now consider the third term in the right hand side of (5.15): Substituting the first term of I (1) n+m given in (5.24) into (5.31) and taking the sum over r out of the integration we get, We can easily see that the integration over each σ n+q for q = r generates a factor of −S n+q . On the other hand integration over σ n+r has exactly the structure of third term (3.22) of single soft graviton case with I (1) n+1 given in (3.25). Using the equality of (3.22), (3.38), (3.50) and the second term on the right hand side of (3.41), the expression (5.32) can be written as: (5.33)

JHEP01(2018)090
On the other hand, substituting the second term on the right hand side of (5.24) into (5.31) we get: We shall evaluate the integration over σ n+u by deforming it to ∞. During this deformation we shall encounter poles at σ n+u = σ n+r for 1 ≤ r ≤ m, r = u, and at σ n+u = ∞ but there are no poles at {σ a } for 1 ≤ a ≤ n. The contribution to (5.34) from the residue at the pole at {σ n+r : 1 ≤ r ≤ m} is given by: On the other hand the contribution to (5.34) from the residue at the pole at σ n+u = ∞ is given by: where in the last line we have used the definition of derivative of delta function inside the contour integration as δ ′ (f ) = − 1 f 2 . For evaluating the integral over σ n+v , we shall deform the contour away from {A v,i } towards the infinity and pick the residues at σ n+v = σ n+u as well as the contribution from infinity. There is no pole at σ n+v = σ a . The contribution from the pole at σ n+u is given by On the other hand the contribution from infinity is given by, after using momentum conservation,

JHEP01(2018)090
In the τ → 0 limit, all the r 1 punctures (σ n+1 , · · · , σ n+r 1 ) come within a distance of order τ of each other, and therefore ξ a ∼ 1. Our goal now is to find out how many solutions exist for the variables (σ n+1 , ξ 2 , · · · , ξ r 1 ). We shall now prove that the number of solutions for this set is (n − 2)(r 1 − 1)!. To do this, we write the scattering equations for the group S 1 in terms of the variables in (A.3) to leading order in τ : Adding all the equations and using momentum conservation, we find that all the ξ a 's drop out of the equation and we get a polynomial equation of degree n − 2 for σ n+1 . Hence it has n − 2 solutions. We can now recursively prove that the number of solutions for the set (ξ 2 , · · · , ξ r 1 ) is (r 1 − 1)!. We consider a situation in which the last momentum k n+r 1 is softer than the other r 1 − 1 soft momenta. We then replace k n+r 1 by ω ℓ n+r 1 for a new soft parameter ω and fixed ℓ n+r 1 , and take the limit ω → 0. 5 In this limit, the terms involving k n+r 1 in (A.4)-(A.6) can be ignored and therefore they reduce to the analog of eqs. (A.4)-(A.7) for a cluster of (r 1 − 1) punctures. These by assumption have (r 1 − 2)! solutions for the variables (ξ 2 , · · · , ξ r 1 −1 ). On the other hand, after we factor out an overall factor of ω from the last equation (A.7), it reduces to a polynomial equation for ξ r 1 of degree r 1 − 1 for a given solution of (σ 1 , · · · , σ n+1 , ξ 2 , · · · , ξ r−1 ). Therefore ξ r 1 has (r 1 − 1) solutions. Thus, there are a total of (r 1 − 1)! solutions for the set (ξ 2 , · · · , ξ r ). For finite ω the solutions will change, but we expect their number to remain the same. Since for fixed σ 1 , · · · σ n , σ n+1 has (n − 2) solutions, this proves our claim that the number of solutions for the set (σ n+1 , · · · , σ n+r 1 ) is (n − 2)(r 1 − 1)!.
We can repeat the analysis for the remaining group of punctures with the result that the number of solutions for the punctures in the group S i is (n − 2)(r i − 1)!. The total number of solutions is, thus, given by where the product runs over all clusters containing punctures carrying soft momentum.

JHEP01(2018)090
To see that it gives the correct number of total solutions (n + m − 3)! for the n + m number of gravitons, we need to sum over all the possibilities. This is given by,