Leading multi-soft limits from scattering equations

A Cachazo-He-Yuan (CHY) type formula is derived for the leading gluon, bi-adjoint scalar $\phi^3$, Yang-Mills-scalar and non-linear sigma model $m$-soft factors $S_m$ in arbitrary dimension. The general formula is used to evaluate explicit examples for up to three soft legs analytically and up to four soft legs numerically via comparison with amplitude ratios under soft kinematics. A structural pattern for gluon $m$-soft factor is inferred and a simpler formula for its calculation is conjectured. In four dimensions, a Cachazo-Svr\v{c}ek-Witten (CSW) recursive procedure producing the leading $m$-soft gluon factor in spinor helicity formalism is developed as an alternative, and Britto-Cachazo-Feng-Witten (BCFW) recursion is used to obtain the leading four-soft gluon factor for all analytically distinct helicity configurations.


Introduction
Investigation of soft factors has a rich history, reaching back to the contributions of Low, Weinberg and others [1][2][3][4][5][6][7][8][9][10]. Soft factorization is a universal property of scattering amplitudes. An n-point scattering amplitude A n depends on external momenta k µ i of the i = 1, 2, ..., n ingoing and outgoing scattering particles. If a subset of adjacent external momenta k µ j for ∀j = 1, 2, ..., m with m < n − 3 is taken to zero, for example parametrized as k µ j → τ k µ j and τ → 0, the amplitude is expected to factorize at leading order in τ into a soft factor S m times a lower point amplitude A n−m : A n → S m A n−m + sub-leading in τ. (1.1) Universality in this context means that S m is independent of the remaining lower point amplitude A n−m , such that S m is always the same whenever the same types of m external particles are taken soft within any original amplitude A n . More recently, interest in investigation of soft theorems was refueled [11][12][13][14] as Strominger et al. showed that soft-graviton theorems can be understood from the point of view of BMS symmetry [15][16][17][18][19]. Further study of leading and sub-leading soft theorems in Yang-Mills, gravity, string and supersymmetric theories ensued , partly based on the amplitude formulation due to Cachazo, He and Yuan (CHY) [22]. Double soft theorems have been considered in [53][54][55], and more recently [56][57][58][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75][76][77]. Construction rules for soft factors with multiple soft particles in N = 4 SYM theory appeared in [78]. Work on related topics was also done, like sub-leading collinear limits [80] and investigation of the current algebra at null infinity induced by soft gluon limits [81].
In this note we use the CHY formulation of scattering amplitudes [22,59] to derive the leading m-soft factor S m for gluons, bi-adjoint scalar φ 3 , Yang-Mills-scalar and non-linear sigma model.
We find the m-soft gluon factor in the case when external legs 1, 2, ..., m are soft to be given by the CHY type formula (3.19, 3.20, 3.21, 3.22). We then consider explicit examples, obtain analytic results in cases m = 1, 2, 3 , and check the cases m = 2, 3, 4 numerically via amplitude ratios in four dimensions obtained from the GGT package [79]. Based on these explicit examples, we infer and conjecture a general pattern for the m-soft gluon factor: is the only new contribution that has to be computed to construct S m at a given m.
As an alternative in four dimensions, we also develop a CSW type [20] automated recursive procedure that gives the leading m-soft gluon factor (compare with construction rules in [78]). Finally, we use BCFW recursion [21] to obtain all leading four-soft gluon factors with analytically distinct helicity combinations in four dimensions.
This work is organized as follows. In section 2 we recall the CHY formalism and introduce the soft limit. In section 3 we demonstrate the soft factorization of gluons at any m and obtain our general result. Explicit examples are worked out in section 4 and a simpler evaluation formula is conjectured. Multi-soft factors in scalar φ 3 , Yang-Mills-scalar and non-linear sigma model are discussed in section 5. Appendix A contains a CSW type recursive procedure for m-soft factors in four dimensions. Appendix B contains BCFW results for four-soft gluon factors in four dimensions. 1 The cases P (m+1−r) r,r+1,...,m,m+1 and P (i) i,i−1,...,2,1,n are obtained by simple index exchange after integration.

The CHY formulation of Yang-Mills and the soft limit
We start with the usual n-point formula for the tree-level gluon amplitude [22]: where the CHY integration measure dµ n and the Yang-Mills CHY integrand I Y M n are dµ n = d n σ σ ij σ jk σ ki vol (SL(2, C)) n a=1 a≠i,j,k Moduli differences are abbreviated as σ ab ≡ σ a − σ b and the matrix Ψ is given by with a, b ∈ {1, 2, ..., n}. The k µ are momenta of scattering particles and µ contain the corresponding polarization data. The indices 1 ≤ i < j < k ≤ n as well as 1 ≤ p < q ≤ n in (2.2) are chosen arbitrarily but fixed. Upper and lower indices on matrix Ψ denote removed columns and rows respectively. We would like to consider the case where m external legs with m < n − 3 are going soft simultaneously: As we take τ → 0, it is clear from the structure of matrix Ψ that at leading order in τ the Pfaffian factorizes as: 2 Pf(Ψ pq pq ) → Pf(ψ)Pf(Ψ p,q,1,2,...,m,n+1,n+2,...,n+m p,q,1,2,...,m,n+1,n+2,...,n+m τ =0 ) + subleading in τ, (2.5) possibly up to an overall sign. The 2m × 2m matrix ψ in the first Pfaffian on the right hand side of (2.5) is defined the same way as Ψ, except the indices a, b in the sub-matrices A, B, C are restricted to the subset a, b ∈ {1, 2, ..., m}. Here, to do the expansion along rows we employed the usual recursive formula for the Pfaffian of an anti-symmetric 2n × 2n matrix M : where m ij are elements of matrix M , θ x ≡ θ(x) is the Heaviside step function, and index i can be freely chosen. Alternatively, we could have noticed that τ → 0 reduces matrix Ψ pq pq at leading order to a block matrix structure, with several blocks equal to zero. Factorization (2.5) then directly follows from trivial Pfaffian factorization identities for block matrices.
Note that Pf(ψ) contains terms leading and/or sub-leading in τ , depending on whether it is evaluated on degenerate (σ ab = O(τ ) for some a, b) or non-degenerate (σ ab = O(1) for all a, b) solutions to the scattering equations. However, for our purposes it is only important that for all types of solutions Pf(ψ) contains all leading contributions.
The second Pfaffian on the right hand side of (2.5) is the one we expect in an (n − m)point amplitude as we take τ → 0. Furthermore, we can trivially rewrite and observe the following behavior in scattering equation delta functions The last equation holds since we necessarily have σ cb = O(1) for m + 1 ≤ c ≤ n due to the kinematics in all k µ c being generic and therefore producing non-degenerate configurations of σ c , while all k b = O(τ ) for the soft particles 1 ≤ b ≤ m tend to zero. The behavior of the first 1 ≤ a ≤ m delta functions in (2.8) is more subtle, since we can have σ ab = O(1) or σ ab = O(τ ) in this case. It will be investigated in detail in the next section.
Considering the above, we can structurally rewrite (2.1) at leading order in τ → 0 as Pf(ψ), (2.10) where dµ n−m and I Y M n−m are based on objects with indices in the range {m + 1, m + 2, ..., n}. Of course this alone does not provide a factorization yet, since S m still depends on σ n and σ m+1 , and the delta functions within still depend on all n momenta and σ-moduli. In the following we show that for any m the σ m+1 , ..., σ n dependence in S m drops out at leading order in τ and the amplitude indeed factorizes as A n → S m A n−m + sub-leading in τ . Furthermore, we find that S m only depends on polarizations µ 1 , µ 2 , ..., µ m as well as momenta k µ n , k µ 1 , k µ 2 , ..., k µ m+1 , and establish a CHY type formula for evaluating S m independently of the remaining factored amplitude A n−m .

Factorization of S m for Yang-Mills and the general result
Starting with S m in (2.10) we apply several transformations in order to more conveniently work with this expression. First we rewrite the delta functions making use of the general identity where • is a placeholder for some test function and we employ the specific m × m matrix which for our particular variables and functions of interest yields the effective relation Furthermore, we transform the moduli σ a into a new set of variables ρ and ξ i : which leads to a change of the integration measure as The transformation (3.5) is convenient, since σ a,a+1 = ξ a allows for more direct access to degenerate solutions σ a,a+1 = O(τ ) in the new ξ a variables. To keep expressions short, we will maintain the σ a notation while implying the substitution (3.5). With the above changes, S m becomes Pf(ψ). (3.7) Now consider keeping ρ fixed and integrating out the q = 1, 2, ..., m−1 constraints h q = 0 (which we will denote as {h} = 0) in the ξ variables. This introduces a Jacobian det(H) −1 with derivative matrix H ij = ∂ ξ i h j and a summation over all solutions to the set of m − 1 equations {h} = 0 in the ξ variables: Pf(ψ). (3.8) Clearly, here all expressions in the integrand can be effectively thought of as functions of the single variable ρ, since σ a = σ a (ρ, {ξ(ρ)}) for a ∈ {1, 2, ..., m} for each solution of {h} = 0 in ξ variables. Therefore, we can now map the single remaining delta function to a simple pole Pf(ψ), (3.9) and consider contour deformations away from the initial locus ∑ m a=1 ∑ n b=m+1 ka⋅k b σ ab = 0 in ρ. By simple power counting of poles we see that there is no pole and therefore no residue at infinity in ρ. As we deform the contour in ρ, the expressions {h} change dynamically since they depend on ρ directly and through ξ(ρ) variables. When we localize ρ at a pole contained in the integrand, the {h} = 0 constraints can get rescaled and simplified. However, since we are summing over the solutions, the set of constraints {h} = 0 has to stay analytic to leading order at the poles in ρ at all times. This implies i.e. that the Jacobian det(H) −1 can get rescaled and simplified due to the contour deformation, but may never diverge. This is a powerful constraint that allows us to find all integrand poles in ρ as follows.
Structurally, the only type of poles that exists in the integrand is of the shape 1 σ ab . As one such pole becomes localized, corresponding terms in the set of expressions {h} start to diverge. Maintaining analyticity at leading order of the divergence in one of the {h} = 0 constraints then demands that at least one different independent 1 σ cd pole must become localized as well simultaneously and at the same rate. 3 This second pole then threatens the analyticity in another {h} = 0 constraint which is affected only by this new divergence, etc. In this fashion a chain of relations occurs demanding that more and more poles must be localized at the same rate simultaneously until it is ensured that analyticity in all {h} = 0 constraints at leading order in the poles is preserved. Overall we realize that whenever a 1 σ ab pole is localized due to the d.o.f. in ρ, the ρ dependence in contributing {ξ(ρ)} solutions must be such that other (m − 1) independent poles become localized as well simultaneously to maintain analyticity in all the {h} = 0 constraints at leading order of divergence.
Equipped with the above observations, we must consider simultaneously localizing subsets of m independent 1 σ ab poles in the integrand, with a ≠ b pairs a, b ∈ {n, 1, 2, ..., m, m + 1}. The only part of the integrand which can diverge more or less dependent on the particular choice of m localized independent 1 σ ab poles is the Parke-Taylor like factor, while all other terms have a fixed scaling (for a given integer m). In the following we consider the case of highest divergence, where combinations of m poles that are present in the Parke-Taylor like factor are localized. There are m+1 m = m + 1 such pole combinations. We will see that this leads to a simple pole overall, such that any other combination of m localized poles does not develop an overall divergence or residue and thus does not contribute. Therefore, localizing m-pole combinations that are present in the Parke-Taylor like factor gives the only non-vanishing contributions.
In the ρ and ξ i variables the Parke-Taylor like factor reads: . where all appearing poles are localized except for: .
We choose to parametrize the m localized poles in the above three cases by a parameter ρ → 0 as follows: The new variablesξ i account for the original degrees of freedom of ξ i variables at leading order after localizingρ → 0. Note that in all three cases we have dρ = dρ, and the one pole that is not localized always directly reduces to 1 σ n,m+1 underρ → 0, which cancels the numerator in (3.10). In general, if we define 4 then, for all possible pole combinations, the behavior of (3.10) forρ → 0 is parametrized as where index r ∈ {1, 2, ..., m + 1} labels which one of the m + 1 poles in the denominator of (3.10) is not being localized. Similarly, for all m + 1 possible pole combinations we obtain with the same index r. Depending on the particular value of r we also get 5 where now H r and ψ r only contain terms supported on the localized poles appearing in the Parke-Taylor like factor (3.10) for each r. It is only at this point that the scattering equations {h r } = 0, their Jacobian 1 det(H r ) and all other terms become completely factorized from the remaining (n − m)-point amplitude A n−m . This means H r and ψ r only depend on momenta k µ n , k µ 1 , k µ 2 , ..., k µ m+1 and polarizations µ 1 , µ 2 , ..., µ m , as expected. Plugging the above findings into (3.9) and collecting the overall power ofρ we observe so that it is now trivial to compute the residues inρ, since for all r we just have a single simple pole atρ = 0. The result is Under closer inspection we note that the Pfaffian factorizes as Pf(ψ r ) = Pf ψ with definitions (3.22), again due to trivial factorization properties of Pfaffians of block matrices with some zero blocks. In principle, (3.18) is already the final completely factorized result. For convenience, we can rewrite it by reassembling the Jacobian and the sum over solutions back into a shape of delta function integrations. This leads to our final general formula: 6 where, identifying k µ 0 ≡ k µ n and keepingσ 0 ≡σ n =σ m+1 = 0 andσ m = 2 −σ 1 in mind, we have with θ x ≡ θ(x) being the Heaviside step function. We call the constraints h q,r = 0 the soft scattering equations.
can be written explicitly as and with indices in the range a, b ∈ {i, i + 1, ..., j}. This is the final result for the m-soft gluon theorem in CHY formulation. We emphasize that the result is correct to leading order in τ → 0. However, since (3.21) admits different solutions of typesσ a,b = O(1) andσ a,b = O(τ ), the integrations in (3.19) have to be evaluated before the result can be systematically expanded to leading order in τ .

Explicit examples and general pattern
In this section we work out examples for the first few soft factors S m . The factors S 1 , S 2 and S 3 are obtained analytically. The factor S 4 (and higher) involves solutions to soft scattering equations that cannot be solved in terms of radicals, therefore we verify the validity of S 4 numerically. Based on the considered examples, we infer a non-trivial structural pattern for the m-soft factors which we conjecture to hold for any m.

One-soft gluon factor S 1
For m = 1 there are no soft scattering equations (3.21) and no delta functions to integrate. The result is just directly given by the sum over r in (3.19): 7 which clearly is the correct Weinberg soft factor. 8 We see that the soft factor is composed out of two pieces such as: Anticipating the structure of higher m-soft factors, we also define Using (4.3) and (4.2) we can structurally write the Weinberg soft factor (4.1) as Based on this and further explicit results of this section, we propose in (4.21) that this structure generalizes and persists for all higher m-soft factors.
Restricting to four dimensions, we can convert the soft factor S 1 to spinor helicity formalism. We use the following standard dictionary to convert expressions of given helicity: where r i and r j label reference spinors assigned to spinor i and j respectively. With an appropriate choice of reference spinor, we see in four dimensions: which is the expected familiar single soft factor in spinor helicity formalism. For real momenta, S − 1 is given by complex conjugation of S + 1 . Here we have suppressed an overall factor of √ 2 in S + 1 per usual spinor helicity convention.

Two-soft gluons factor S 2
For m = 2, there is one soft scattering equation (3.21) for each r, and the number of solutions organizes as follows for the different solution types and different values of r: Adding up the contributions of all 5 solutions and expanding to leading order in τ , we obtain the following expression for S 2 : This agrees with the generalization (4.21) for m = 2. The quantities P Counting the powers of k 1 and k 2 we see that this expression diverges as τ −2 , as we expect for the two-soft gluon factor. The result (4.9) is gauge independent and reduces to the gauge fixed result found in [61] when we select the gauge 2 ⋅ k 3 = 0, 1 ⋅ k n = 0.
Restricting to four dimensions, converting to spinor helicity formalism by use of (4.5) and (4.6), and choosing appropriate reference spinors we get the following expression for the non-trivial helicity combination (+−) after some simplification via Schouten identities: which naturally agrees with the result found in [61]. The trivial helicity combination (++) reduces to the product of single soft factors S ++ 2 = ⟨n3⟩ ⟨n1⟩⟨12⟩⟨23⟩ as expected. Again, an overall factor of ( √ 2) 2 is suppressed in the above expressions per spinor helicity convention and the other helicity combinations can be obtained by complex conjugation.
We can additionally numerically test the above result in four dimensions. Making use of the GGT package provided in [79] to generate explicit lower point amplitudes, we can 9 Here, for brevity we use that 2(k1 + k2) ⋅ k3 ≈ s123 at leading order in τ .
form amplitude ratios that correspond to the soft factor in appropriate soft kinematics. 10 Keeping in mind the overall powers of √ 2 that are suppressed in spinor helicity, we expect to find the following relation at leading order in τ : .
Indeed, if we generate a numeric kinematic point where k µ 1 , k µ 2 have soft entries of order 10 −10 while the rest of the momenta have hard entries of order 10 0 , we can check that i.e.
hold at least to first 10 digits, reflecting that the leading soft factor receives a first correction at the next polynomially sub-leading power in τ . 11 Naturally, ratios of more complicated amplitudes yield the same agreement.
where we imply i ≠ j and i, j ∈ {1, 2}. Adding up the contributions of all 16 solutions and expanding to leading order in τ , we obtain the following expression for S 3 : 1,2,3,4 P (0) n − P This agrees with the generalization (4.21) for m = 3. As before, expressions of type P i,j and P (2) i,j,l are given by (4.3), (4.2) and (4.10), while the new contribution of type P (3) i,j,l,t can still be analytically computed to be: 12  10 Note that there is a Chop command in one of the routines of the GGT package, which does not work well with soft limit numerics and therefore needs to be removed. 11 To make sure that the comparison works properly, we use the same spinor conventions as the GGT package: 12 Again, we use that 2(k1 + k2 + k3) ⋅ k4 ≈ s1234 and similar at leading order in τ to keep notation short.
where we used the abbreviations Counting the powers of k 1 , k 2 and k 3 we see that this expression diverges as τ −3 , as we expect for the three-soft gluon factor.
Again, we can use (4.5) and (4.6) to pass to spinor helicity formalism if we restrict to four dimensions. In particular, the two non-trivial independent polarization combinations are (− + −) and (+ − −). For the case (− + −) we obtain, with appropriate choice of reference spinors and after some simplification via Schouten identities:  [14] . The trivial helicity configuration (+ + +) as expected reduces to S +++ 3 = ⟨n4⟩ ⟨n1⟩⟨12⟩⟨23⟩⟨34⟩ , and all other helicity configurations are obtained from the above by symmetry and complex conjugation. An overall factor of 2 3 2 is suppressed in the above expressions per spinor helicity convention.
As before, (4.12) is expected to hold. Making use of the GGT package [79] to generate explicit lower point amplitudes we can form ratios that correspond to the soft factor in appropriate soft kinematics. Generating a numeric kinematic point such that k µ 1 , k µ 2 and k µ 3 have soft entries of order 10 −10 while the rest of the momenta have hard entries of order 10 0 , we observe that i.e.
hold to at least the first 10 digits, after which the first sub-leading correction in τ becomes important. Again, ratios of more complicated amplitudes yield the same agreement.

Four-soft gluons factor S 4 and beyond
For m = 4, there are three soft scattering equations (3.21) for each r, and the number of solutions organizes as follows for the different solution types and different values of r: where we imply i ≠ j, i ≠ l, j ≠ l and i, j, l ∈ {1, 2, 3}. With the generalization (4.21) in mind, we expect that the contributions for cases r = 2, 3, 4 can be constructed from previously determined quantities (4.2), (4.10) and (4.16). That is easily verified numerically by obtaining and summing over explicit approximate solutions to the soft scattering equations (3.21) in some example kinematics. This confirms that the structure 1,2,3,4,5 P (0) n − P continues to hold. Trying to obtain P 1,2,3,4,5 for r = 1 (and r = 5) we discover that finding the 12 solutions of the typeξ 1 ∼ξ 2 ∼ξ 3 ∼ O(τ ) is equivalent to solving for the roots of two 6th degree polynomials. Therefore, an analytic solution cannot be obtained in this direct fashion.
Based on the knowledge of previous analytic results found so far, we could try to infer the pole structure of all the different terms appearing in P (4) 1,2,3,4,5 , effectively constructing the result without solving the soft scattering equations. This works reasonably well for some of the appearing terms such as 1 ⋅ k 2 2 ⋅ k 3 3 ⋅ k 4 4 ⋅ k 5 , for which the correct contribution can be guessed (and numerically checked) to be given by: However, P 1,2,3,4,5 also contains terms such as 3 ⋅ 4 1 ⋅ k 2 2 ⋅ k 5 or 1 ⋅ 2 3 ⋅ 4 for which the pole structure is unclear since these patterns did not appear before. Even though an analytic solution is thus not available, we can still check numerically that (3.19) is correct.
Using (4.5) and (4.6) to pass to spinor helicity formalism in four dimensions, (4.12) is again expected to hold. Therefore, we generate a numeric kinematic point such that k µ 1 , k µ 2 , k µ 3 and k µ 4 have soft entries of order 10 −10 while the rest of the momenta have hard entries of order 10 0 . Now we can solve (3.21) numerically and obtain the numeric soft factor S 4 as a sum over all 64 solutions. Subsequently, making use of the GGT package [79], we can generate explicit amplitude ratios and observe that e.g.
hold to at least the first 10 digits, after which the first sub-leading correction in τ becomes important. As before, ratios of more complicated amplitudes yield the same agreement. For even higher m, the soft scattering equations (3.21) become more and more complicated, so that even numeric evaluation becomes increasingly harder to do. However, in principle the m-soft gluon factor is always given by the CHY type expression summarized by (3.19), (3.21) and (3.22), valid to leading order in τ .

Conclusion and general structural pattern
The above findings are of interest since they prove the existence of a universal soft factor for any number of soft adjacent gluons and in principle provide a way to calculate these soft factors in arbitrary dimension. As a byproduct we obtained an explicit analytic result for the three-soft gluon factor for arbitrary polarizations and in arbitrary dimension, which to our knowledge is a new result.
Considering the particular results for m = 1, 2, 3, 4 discussed above, we can infer a generalization for the structural pattern at arbitrary m to be given by: r,r+1,...,m,m+1 P In this sense, it suffices to evaluate only the r = 1 summand of (3.19) to obtain all new information at a given m. 14 The above conjecture is inferred empirically, and it seems to be highly non-trivial to demonstrate the factorization of each summand of (3.19) into (4.21) analytically. While the structure of the Pfaffian admits such a factorization, the Parke-Taylor like factor as well as the multiplicative term remaining from the contour deformation in ρ are not convenient. This implies the necessity of a transformation along the lines of (3.1) with a non-trivial Jacobian, which is not easily guessed. We leave a general proof of the conjecture (4.21), (4.22) to future work.

Multi-soft factors in other theories
It is possible to directly apply the procedure described above to several other theories in CHY formulation. An important feature that largely governs the computations is the presence of at least one Parke-Taylor factor C ≡ 1 σ 12 σ 23 ...σ n1 (5.1) in the CHY integrand of the amplitude, such that the amplitude in question is color ordered. The theories considered in this section have this same feature. As further building blocks we will require the sub-matrix A defined in (2.3), the matrix Ψ n+1,n+2,...,n+q n+1,n+2,...,n+q which is (2.3) with rows and columns n + 1, n + 2, ..., n + q dropped, and the matrix where I a , I b are some internal space indices for scalar fields involved in the scattering process [59]. Since these indices have no non-trivial effect on the momentum dependence of soft factors, we will consider the simplest case where I a = I b for all particle labels a, b , such that δ Ia,I b = 1.

Multi-soft factors in bi-adjoint scalar φ 3 theory
The CHY formula for tree level scattering in bi-adjoint scalar φ 3 theory can be written as (2.1) [22] with I Y M n replaced by Starting with this integrand, the considerations in sections 2 and 3 go through in the same manner, such that we are left with the following general expression for the m-soft scalar factor with particles 1, 2, ..., m going soft: with dν r given in (3.20), and the identification σ 0 ≡ σ n . As in the gluon case, the soft scattering equations contained in dν r can be explicitly solved for the cases m = 1, 2, 3, with exactly the same solutions. At leading order in the soft limit this leads to It is worth noticing that all contributions to the soft factors at leading order in the soft limit are due to the two summands r = 1 and r = m + 1 only, while the intermediate summands are sub-leading. As before, the general expression S φ 3 m can be used to evaluate S φ 3 4 and higher soft factors numerically. We tested the results numerically against amplitude ratios in CHY formulation and found agreement.

Multi-soft factors in Yang-Mills-scalar theory
The CHY formula for tree level scattering in Yang-Mills-scalar theory is (2.1) with I Y M n replaced by j,n+1,n+2,...,n+q i,j,n+1,n+2,...,n+q ), where matrix χ is q × q dimensional (5.2), and 1 ≤ i < j ≤ n can be selected arbitrarily [59]. This corresponds to the first q of the scattering particles being scalars and the remaining n − q being gluons. Starting with this integrand, the considerations in sections 2 and 3 go through in the same manner. Soft gluon factors in this theory are exactly the same as in pure Yang-Mills. The general expression for the m-soft scalar factor with particles 1, 2, ..., m going soft amounts to: 15 with dν r given in (3.20), and the identification σ 0 ≡ σ n . The matrix A [i,j] was defined in (3.22), and the matrix χ [i,j] relates to χ in (5.2) the same as A [i,j] relates to A in (2.3).
As in the gluon case, the soft scattering equations contained in dν r can be explicitly solved for the cases m = 1, 2, 3, with exactly the same solutions. However, since Pf(χ [i,j] ) vanishes when χ [i,j] is of odd dimension, only soft factors with an even number m of soft scalars are non-zero and only summands of odd r contribute. At leading order in the soft limit this leads to This agrees with the result in [60]. As before, the general expression S Y M S m can be used to evaluate S Y M S 4 and higher soft factors numerically. We tested the results numerically against amplitude ratios in CHY formulation and found agreement.

Multi-soft factors in non-linear sigma model
The CHY formula for tree level scattering in non-linear sigma model is (2.1) with I Y M n replaced by where A i,j i,j is the matrix A defined in (2.3) with rows and columns i, j removed, and 1 ≤ i < j ≤ n can be selected arbitrarily [59]. Starting with this integrand, the considerations in sections 2 and 3 go through in the same manner. The general expression for the m-soft factor with particles 1, 2, ..., m going soft amounts to: with dν r given in (3.20), and the identification σ 0 ≡ σ n . The matrix A [i,j] was defined in (3.22). As in the gluon case, the soft scattering equations contained in dν r can be explicitly solved for the cases m = 1, 2, 3, with exactly the same solutions. However, since Pf(A [i,j] ) vanishes when A [i,j] is of odd dimension, only soft factors with an even number m of soft particles are non-zero and only summands of odd r contribute. At leading order in the soft limit this leads to This agrees with the result in [60]. As before, the general expression S N LSM m can be used to evaluate S N LSM 4 and higher soft factors numerically. We tested the results numerically against amplitude ratios in CHY formulation and found agreement. Additionally, our S N LSM 4 numerically agrees with the result found in [82]. 16

A CSW recursion for multi-gluon soft-factors in four dimensions
As an alternative to the construction rules presented in [78], we can set up a CSW type recursion [20] for the m-soft factors in four dimensions as follows. We start with the amplitude A (m) (k +1 n , k h 1 1 , ..., k hm m , k +1 m+1 ), where k h i i denotes the external momentum of the i-th particle with helicity h i ∈ {+1, −1}. Here we have cyclically rotated the n-th position to be the first, and suppressed all entries k h j j with m + 1 < j < n since they do not enter the soft factor that we want to extract from this amplitude. Since the helicities of particle n and m + 1 do not enter the soft factor, we can choose these helicities to be + without loss of generality. The superscript (m) keeps track of the number of adjacent external momenta that are taken soft.
In order to obtain the soft factor from CSW recursion, we have to generate all possible diagrams in MHV expansion. To do this recursively, we introduce the following two functions: as well as, making use of µ(x) ≡ mod(x − 1, l) + 1, the function: We supplement the above functions with the following resolution properties: which ensure that the explicit propagator momenta always are properly resolved in terms of external momenta. Naturally, the order of indices i, ..., j appearing in p(i, ..., j) and P 2 i,...,j is irrelevant and can be assumed to be sorted to make it easier to identify and group together identical expressions.
It is important to note that the sums in the functions (A.1) and (A.2) may contain summands that immediately vanish due to trivial helicity configurations of sub-amplitudes involved that enter the H function. 17 Setting such summands to zero directly without allowing for any recursion depth in such terms greatly speeds up the calculation.
Recursion by means of (A.1) and (A.2) with the above supplements will generate all possible diagrams in MHV expansion that contribute to leading order in the soft limit. However, the simple summation employed here comes at the expense of multiple counting for some of the resulting diagrams. The easiest way to remove the over-counting is to simply set the integer coefficient in front of each overall summand to 1 after the recursion has been completed and all terms have been properly grouped together: S ′ ≡ S with multiplicity of each overall summand set to 1, which implies that invariance of amplitudes under cyclic permutation of external legs is used to identify and group together equivalent terms in the expansion. This, as well as the entire recursive procedure, can be easily automated i.e. in Mathematica, such that the m-soft factor S m for any helicity configuration is automatically generated by the input: By trivial helicity configuration we mean amplitudes with none, or only one negative helicity gluon, as well as amplitudes with none, or only one positive helicity gluon (special care is required for 3-point amplitudes due to special kinematics).
Finally, to evaluate the soft factor we use the substitutions where entries like p(i, ..., j)⟩ are evaluated by the usual CSW replacement P i,...,j X] with reference spinor X]. Superficially, due to (A.9) it might seem that the soft factor depends on (n − 1)-st and (m + 2)-nd external momentum as well. However, just as in [78], this dependence always cancels out upon the CSW replacement of the shifted spinors at leading order in τ .
We have tested the above recursive procedure for soft factors S 1 , S 2 , ..., S 7 with various helicity configurations against appropriate amplitude ratios obtained from the GGT package [79], and found numerical agreement at leading order in τ . For example, our recursion takes about two minutes to generate the 2277 different analytic terms in the S −−−−−−+ 7 soft factor. If required, a trivial further expansion in τ can be used to isolate leading terms only.

B Four-soft gluons from BCFW
Naturally, it is also possible to apply BCFW recursion relations [21] to compute higher soft factors. Here we demonstrate the four-soft gluon calculation. We pick gluons 1, 2, 3, 4 to be soft and perform a [23⟩ BCFW shift, so that 2 →2 and 3 →3 with in each case. Here, A 3 , A 4 are mostly-soft-leg sub-amplitudes factored by BCFW, and S 2 , S 3 are two-and three-soft gluon factors that are extracted from the mostly-hard-leg sub-amplitudes factored by BCFW. The usual on-shell constraintsP 2 ⋯ = 0 provide the following z values to leading order in the soft limit: 18 , z B = − [12] [13] , z C = s 345 ⟨25⟩ [53] , In case when all four soft gluons have the same helicity, the four-soft factor trivially reduces to a product of consecutive soft factors. In the following, we specify explicit helicity configurations and obtain the results for all analytically distinct non-trivial helicity configurations.
Helicity configuration (− + ++): For the helicity configuration of soft gluons (1 − , 2 + , 3 + , 4 + ) we find: To see that the diagram C is zero, we use the fact that the soft factor is independent of the helicity of particle 5, thus we can choose it to be 5 + which leads to no non-vanishing helicity configurations for A 4 . In all other diagrams only one helicity configuration is non-vanishing. We tested the above result numerically against amplitude ratios and found agreement.

Helicity configuration (+ − ++):
For the helicity configuration of soft gluons (1 + , 2 − , 3 + , 4 + ) we find: Diagram C vanishes the same way as described above. In all other diagrams again only one helicity configuration is non-vanishing. We tested the above result numerically against amplitude ratios and found agreement. 18 We use the convention sij = ⟨ij⟩[ji], which with our spinor contraction conventions (⟨ij⟩ = λ 1 i λ 2 j − λ  In all diagrams again only one helicity configuration is non-vanishing. We tested the above result numerically against amplitude ratios and found agreement.  Diagram C vanishes the same way as described above. In all other diagrams again only one helicity configuration is non-vanishing. We tested the above result numerically against amplitude ratios and found agreement. For the helicity configuration of soft gluons (1 + , 2 − , 3 + , 4 − ) we find: In all diagrams again only one helicity configuration is non-vanishing. We tested the above result numerically against amplitude ratios and found agreement.