---
author:
- Michael Zlotnikov
title: 'Leading multi-soft limits from scattering equations'
---

Introduction
============

Investigation of soft factors has a rich history, reaching back to the contributions of Low, Weinberg and others [@LowFourPt1; @LowFourPt2; @LowFourPt3; @LowTheorem; @Weinberg1; @weinberg; @OtherSoftPhotons1; @Gross:1968in; @OtherSoftPhotons2; @OtherSoftPhotons3]. Soft factorization is a universal property of scattering amplitudes. An $n$-point scattering amplitude $A_n$ depends on external momenta $k_i^\mu$ of the $i=1,2,...,n$ ingoing and outgoing scattering particles. If a subset of adjacent external momenta $k_j^\mu$ for $\forall j=1,2,...,m$ with $m<n-3$ is taken to zero, for example parametrized as $k_j^\mu\to\tau k_j^\mu$ and $\tau\to 0$, the amplitude is expected to factorize at leading order in $\tau$ into a soft factor $S_m$ times a lower point amplitude $A_{n-m}$: $$\begin{aligned}
A_n\to S_m A_{n-m}+\text{sub-leading in }\tau.\end{aligned}$$ 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 [@strominger1; @strominger2; @strominger3; @Cachazo:2014fwa] as Strominger et al. showed that soft-graviton theorems can be understood from the point of view of BMS symmetry [@BMS1; @BMS2; @ExtendedBMS1; @ExtendedBMS2; @ExtendedBMS3]. Further study of leading and sub-leading soft theorems in Yang-Mills, gravity, string and supersymmetric theories ensued [@Cachazo:2013hca; @schwab; @nima; @1405.1015; @1405.1410; @Broedel:2014fsa; @Bern:2014vva; @Casali:2014xpa; @Larkoski:2014hta; @1406.7184; @Zlotnikov:2014sva; @Kalousios:2014uva; @1408.4179; @Lipstein:2015rxa; @Avery:2015gxa; @1509.07840; @Bianchi:2015lnw; @Avery:2015iix; @Bianchi:2016tju; @DiVecchia:2016amo; @1509.01406; @1605.09094; @1608.00685; @1610.03481; @1702.03934; @1703.00024; @1706.00759; @Chakrabarti:2017ltl; @Mao:2017wvx; @1611.02172; @1611.07534], partly based on the amplitude formulation due to Cachazo, He and Yuan (CHY) [@Cachazo:2013hca]. Double soft theorems have been considered in [@weinberg2; @JS; @adler], and more recently [@ArkaniHamed:2008gz; @Chen:2014cuc; @Chen:2014xoa; @Cachazo:2015ksa; @Volovich:2015yoa; @Klose:2015xoa; @Du:2015esa; @DiVecchia:2015bfa; @1507.08829; @1505.05854; @1502.05258; @1406.4172; @1406.5155; @1411.6661; @DiVecchia:2015srk; @Low:2015ogb; @He:2016vfi; @Saha:2016kjr; @Saha:2017yqi; @1604.02834; @1705.06175]. Construction rules for soft factors with multiple soft particles in $\mathcal{N}=4$ SYM theory appeared in [@Georgiou:2015jfa]. Work on related topics was also done, like sub-leading collinear limits [@Nandan:2016ohb] and investigation of the current algebra at null infinity induced by soft gluon limits [@McLoughlin:2016uwa].

In this note we use the CHY formulation of scattering amplitudes [@Cachazo:2013hca; @Cachazo:2014xea] to derive the leading $m$-soft factor $S_m$ for gluons, bi-adjoint scalar $\phi^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 (\[softfacres\], \[softmeasure\], \[softscat\], \[softps\]). 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 [@Dixon:2010ik]. Based on these explicit examples, we infer and conjecture a general pattern for the $m$-soft gluon factor: $$\begin{aligned}
\label{smres0}
S^{gluon}_m&=\sum_{r=1}^{m+1}(-1)^{r+1}P^{(m+1-r)}_{r,r+1,...,m,m+1}P^{(r-1)}_{r-1,r-2,...,1,n}\,,\end{aligned}$$ where $P^{(0)}_{m+1}=P^{(0)}_{n}\equiv 1$, and $P^{(i)}_{1,2,...,i,i+1}$, with $d\nu_1$ and $\psi_{[1,i]}^{(i+1)}$ defined in (\[softmeasure\]) and (\[softps\]), is[^1] $$\begin{aligned}
\label{Pdef0}
P^{(i)}_{1,2,...,i,i+1}=&\int d\nu_1\frac{1}{\prod_{c=2}^{i+1}\bar\sigma_{c-1,c}}\text{Pf}\left(\psi_{[1,i]}^{(i+1)}\right).\end{aligned}$$ If all $P^{(i)}_{1,2,...,i,i+1}$ with $i<m$ are known from calculations of lower soft factors, then $P^{(m)}_{1,2,...,m,m+1}$ is the only new contribution that has to be computed to construct $S_m$ at a given $m$.

The leading $m$-soft factor in bi-adjoint scalar $\phi^3$, Yang-Mills-scalar and non-linear sigma model theories involves the same integration measure $d\nu_r$ as in (\[softfacres\]), while the integrands are different: (\[SmScalar\]), (\[SmYMS\]) and (\[SmNLSM\]).

As an alternative in four dimensions, we also develop a CSW type [@Cachazo:2004kj] automated recursive procedure that gives the leading $m$-soft gluon factor (compare with construction rules in [@Georgiou:2015jfa]). Finally, we use BCFW recursion [@Britto:2004ap] to obtain all leading four-soft gluon factors with analytically distinct helicity combinations in four dimensions.

This work is organized as follows. In section \[CHYsoft\] we recall the CHY formalism and introduce the soft limit. In section \[factorization\] we demonstrate the soft factorization of gluons at any $m$ and obtain our general result. Explicit examples are worked out in section \[examples\] and a simpler evaluation formula is conjectured. Multi-soft factors in scalar $\phi^3$, Yang-Mills-scalar and non-linear sigma model are discussed in section \[othertheories\]. Appendix \[CSWrecur\] contains a CSW type recursive procedure for $m$-soft factors in four dimensions. Appendix \[foursoftBCFW\] contains BCFW results for four-soft gluon factors in four dimensions.

The CHY formulation of Yang-Mills and the soft limit {#CHYsoft}
====================================================

We start with the usual $n$-point formula for the tree-level gluon amplitude [@Cachazo:2013hca]: $$\begin{aligned}
\label{ampl}
A_n=&\int d\mu_n ~{\cal I}_n^{YM},\end{aligned}$$ where the CHY integration measure $d\mu_n$ and the Yang-Mills CHY integrand ${\cal I}_n^{YM}$ are $$\begin{aligned}
\label{CHYmeasure}
d\mu_n=&\int d^n\sigma\frac{\sigma_{ij}\sigma_{jk}\sigma_{ki}}{\text{vol}\left(SL(2,C)\right)}\prod_{\substack{a=1\\a\neq i,j,k}}^n\delta\left(\sum_{\substack{b=1\\b\neq a}}^n\frac{k_a\cdot k_b}{\sigma_{ab}}\right)~~~,~~~{\cal I}_n^{YM}=\frac{2\frac{(-1)^{p+q}}{\sigma_{pq}}\text{Pf}(\Psi^{pq}_{pq})}{\sigma_{12}\sigma_{23}...\sigma_{n1}}.\end{aligned}$$ Moduli differences are abbreviated as $\sigma_{ab}\equiv \sigma_a-\sigma_b$ and the matrix $\Psi$ is given by $$\begin{aligned}
\label{psimat}
\Psi=\left({{A}\atop{C}}~{{-C^T}\atop{B}}\right),~~A=\left\{{{\frac{k_a\cdot k_b}{\sigma_{ab}}}\atop{0}}{{;}\atop{;}}{{a\neq b}\atop{a=b}}\right.,~~B=\left\{{{\frac{\epsilon_a\cdot \epsilon_b}{\sigma_{ab}}}\atop{0}}{{;}\atop{;}}{{a\neq b}\atop{a=b}}\right.,~~C=\left\{\begin{matrix}
\frac{\epsilon_a\cdot k_b}{\sigma_{ab}}&;\, a\neq b\\
-\sum_{\substack{c=1\\c\neq a}}^n\frac{\epsilon_a\cdot k_c}{\sigma_{ac}}&;\, a=b
\end{matrix}\right.,\end{aligned}$$ with $a,b\in\{1,2,...,n\}$. The $k^\mu$ are momenta of scattering particles and $\epsilon^\mu$ contain the corresponding polarization data. The indices $1\leq i<j<k\leq n$ as well as $1\leq p<q\leq n$ in (\[CHYmeasure\]) are chosen arbitrarily but fixed. Upper and lower indices on matrix $\Psi$ 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: $$\begin{aligned}
\label{taulim}
k^\mu_q\to \tau k^\mu_q,~~~\tau\to 0,~~~\text{for}~~~q\in \{1,2,...,m\}.\end{aligned}$$ As we take $\tau\to 0$, it is clear from the structure of matrix $\Psi$ that at leading order in $\tau$ the Pfaffian factorizes as:[^2] $$\begin{aligned}
\label{pffactor}
\text{Pf}(\Psi^{pq}_{pq})\to\text{Pf}(\psi)\text{Pf}(\Psi^{p,q,1,2,...,m,n+1,n+2,...,n+m}_{p,q,1,2,...,m,n+1,n+2,...,n+m}|_{\tau=0})+\text{subleading in }\tau,\end{aligned}$$ possibly up to an overall sign. The $2m\times 2m$ matrix $\psi$ in the first Pfaffian on the right hand side of (\[pffactor\]) is defined the same way as $\Psi$, except the indices $a,b$ in the sub-matrices $A,B,C$ are restricted to the subset $a,b\in\{1,2,...,m\}$. Here, to do the expansion along rows we employed the usual recursive formula for the Pfaffian of an anti-symmetric $2n\times 2n$ matrix $M$: $$\begin{aligned}
{\rm Pf}\left(M\right)=\sum_{{j=1}\atop{j\neq i}}^{2n}(-1)^{i+j+1+\theta_{i-j}}m_{ij}{\rm Pf}\left(M^{ij}_{ij}\right),\end{aligned}$$ where $m_{ij}$ are elements of matrix $M$, $\theta_x\equiv\theta(x)$ is the Heaviside step function, and index $i$ can be freely chosen.\
Alternatively, we could have noticed that $\tau\to 0$ reduces matrix $\Psi^{pq}_{pq}$ at leading order to a block matrix structure, with several blocks equal to zero. Factorization (\[pffactor\]) then directly follows from trivial Pfaffian factorization identities for block matrices.\
Note that $\text{Pf}(\psi)$ contains terms leading and/or sub-leading in $\tau$, depending on whether it is evaluated on degenerate ($\sigma_{ab}=O(\tau)$ for some $a,b$) or non-degenerate ($\sigma_{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 $\text{Pf}(\psi)$ contains all leading contributions.\
The second Pfaffian on the right hand side of (\[pffactor\]) is the one we expect in an $(n-m)$-point amplitude as we take $\tau\to0$. Furthermore, we can trivially rewrite $$\begin{aligned}
\frac{1}{\sigma_{12}\sigma_{23}...\sigma_{n1}}=\frac{\sigma_{n,m+1}}{\sigma_{n1}\sigma_{12}...\sigma_{m,m+1}}\cdot\frac{1}{\sigma_{n,m+1}\sigma_{m+1,m+2}...\sigma_{n-1,n}},\end{aligned}$$ and observe the following behavior in scattering equation delta functions $$\begin{aligned}
\label{deltafact1}
\prod_{a=1}^n\delta\left(\sum_{{b=1}\atop{b\neq a}}^n\frac{k_a\cdot k_b}{\sigma_{ab}}\right)=\prod_{a=1}^m\delta\left(\sum_{{b=1}\atop{b\neq a}}^n\frac{k_a\cdot k_b}{\sigma_{ab}}\right)\prod_{c=m+1}^n\delta\left(\sum_{{b=m+1}\atop{b\neq c}}^n\frac{k_c\cdot k_b}{\sigma_{cb}}+O(\tau)\right).\end{aligned}$$ The last equation holds since we necessarily have $\sigma_{cb}=O(1)$ for $ m+1\leq c\leq n$ due to the kinematics in all $k^\mu_c$ being generic and therefore producing non-degenerate configurations of $\sigma_c$, while all $k_b=O(\tau)$ for the soft particles $1\leq b\leq m$ tend to zero. The behavior of the first $1\leq a\leq m$ delta functions in (\[deltafact1\]) is more subtle, since we can have $\sigma_{ab}=O(1)$ or $\sigma_{ab}=O(\tau)$ in this case. It will be investigated in detail in the next section.\
Considering the above, we can structurally rewrite (\[ampl\]) at leading order in $\tau\to 0$ as $$\begin{aligned}
A_n\to&\int d\mu_{n-m} ~S_m~{\cal I}_{n-m}^{YM}+\text{sub-leading in }\tau,\\
\label{sdef}
S_m=&\int d^m\sigma\prod_{a=1}^m\delta\left(\sum_{{b=1}\atop{b\neq a}}^n\frac{k_a\cdot k_b}{\sigma_{ab}}\right)\frac{\sigma_{n,m+1}}{\sigma_{n1}\sigma_{12}...\sigma_{m,m+1}}\text{Pf}(\psi),\end{aligned}$$ where $d\mu_{n-m}$ and ${\cal I}_{n-m}^{YM}$ 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 $\sigma_n$ and $\sigma_{m+1}$, and the delta functions within still depend on all $n$ momenta and $\sigma$-moduli. In the following we show that for any $m$ the $\sigma_{m+1},...,\sigma_n$ dependence in $S_m$ drops out at leading order in $\tau$ and the amplitude indeed factorizes as $A_n\to S_m A_{n-m}+$ sub-leading in $\tau$. Furthermore, we find that $S_m$ only depends on polarizations $\epsilon_1^\mu,\epsilon_2^\mu,...,\epsilon_m^\mu$ as well as momenta $k_n^\mu,k_1^\mu,k_2^\mu,...,k_{m+1}^\mu$, 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 {#factorization}
============================================================

Starting with $S_m$ in (\[sdef\]) 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 $$\begin{aligned}
\label{deltatrafo}
\int d^m x\prod_{i=1}^m\delta\left(f_i(\vec x)\right)\,\bullet~ =\int d^m x\det(M)\prod_{i=1}^m\delta\left(\sum_{j=1}^m M_{ij}f_j(\vec x)\right)\,\bullet~,\end{aligned}$$ where $\bullet$ is a placeholder for some test function and we employ the specific $m\times m$ matrix $$\begin{aligned}
M=\left(\begin{matrix}
1 & 1 & 1 & ... & 1 & 1 & 1 \\
1 & -1 & 0 & ... & 0 & 0 & 0 \\
0 & 1 & -1 & ... & 0 & 0 & 0 \\
\ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots \\
0 &0 & 0 & ... & 1 & -1 & 0 \\
0 & 0 & 0 & ... & 0 & 1 & -1 
\end{matrix}\right),~~~~~~\det(M)=(-1)^{m+1}m,\end{aligned}$$ which for our particular variables and functions of interest yields the effective relation $$\begin{aligned}
\label{deltraf}
\prod_{a=1}^m\delta\left(\sum_{{b=1}\atop{b\neq a}}^n\frac{k_a\cdot k_b}{\sigma_{ab}}\right)=&(-1)^{m+1}m\,\delta\left(\sum_{a=1}^m\sum_{b=m+1}^n\frac{k_a\cdot k_b}{\sigma_{ab}}\right)\prod_{q=1}^{m-1}\delta\left(h_q\right),\\
\label{hdef}
h_q=& \sum_{{a=1}\atop{a\neq q}}^n\frac{k_q\cdot k_a}{\sigma_{qa}}-\sum_{{b=1}\atop{b\neq q+1}}^n\frac{k_{q+1}\cdot k_b}{\sigma_{q+1,b}}.\end{aligned}$$ Furthermore, we transform the moduli $\sigma_a$ into a new set of variables $\rho$ and $\xi_i$: $$\begin{aligned}
\label{sigtraf}
\sigma_q=\rho-\sum_{a=1}^{q-1}\frac{\xi_a}{2}+\sum_{b=q}^{m-1}\frac{\xi_b}{2},\end{aligned}$$ which leads to a change of the integration measure as $$\begin{aligned}
d\sigma_1\wedge d\sigma_2\wedge ...\wedge d\sigma_m =(-1)^{m+1}d\rho\wedge d\xi_1\wedge d\xi_2\wedge...\wedge d\xi_{m-1}.\end{aligned}$$ The transformation (\[sigtraf\]) is convenient, since $\sigma_{a,a+1}=\xi_a$ allows for more direct access to degenerate solutions $\sigma_{a,a+1}=O(\tau)$ in the new $\xi_a$ variables. To keep expressions short, we will maintain the $\sigma_a$ notation while implying the substitution (\[sigtraf\]). With the above changes, $S_m$ becomes $$\begin{aligned}
S_m=\int d\rho\, d^{m-1}\xi ~m\,\delta\left(\sum_{a=1}^m\sum_{b=m+1}^n\frac{k_a\cdot k_b}{\sigma_{ab}}\right)\prod_{q=1}^{m-1}\delta\left(h_q\right)\frac{\sigma_{n,m+1}}{\sigma_{n1}\sigma_{12}...\sigma_{m,m+1}}\text{Pf}(\psi).\end{aligned}$$

Now consider keeping $\rho$ fixed and integrating out the $q=1,2,...,m-1$ constraints $h_q=0$ (which we will denote as $\{h\}=0$) in the $\xi$ variables. This introduces a Jacobian $\det(H)^{-1}$ with derivative matrix $H_{ij}=\partial_{\xi_i}h_j$ and a summation over all solutions to the set of $m-1$ equations $\{h\}=0$ in the $\xi$ variables: $$\begin{aligned}
S_m=\sum_{{\{h\}=0}\atop{\text{solutions}}}\int d\rho \frac{m}{\det(H)}\,\delta\left(\sum_{a=1}^m\sum_{b=m+1}^n\frac{k_a\cdot k_b}{\sigma_{ab}}\right)\frac{\sigma_{n,m+1}}{\sigma_{n1}\sigma_{12}...\sigma_{m,m+1}}\text{Pf}(\psi).\end{aligned}$$ Clearly, here all expressions in the integrand can be effectively thought of as functions of the single variable $\rho$, since $\sigma_a=\sigma_a\left(\rho,\{\xi(\rho)\}\right)$ for $a\in\{1,2,...,m\}$ for each solution of $\{h\}=0$ in $\xi$ variables. Therefore, we can now map the single remaining delta function to a simple pole $$\begin{aligned}
\label{scontour}
S_m=\sum_{{\{h\}=0}\atop{\text{solutions}}}\oint \frac{d\rho}{2\pi i} \frac{m}{\det(H)}\frac{1}{\sum_{a=1}^m\sum_{b=m+1}^n\frac{k_a\cdot k_b}{\sigma_{ab}}}\frac{\sigma_{n,m+1}}{\sigma_{n1}\sigma_{12}...\sigma_{m,m+1}}\text{Pf}(\psi),\end{aligned}$$ and consider contour deformations away from the initial locus $\sum_{a=1}^m\sum_{b=m+1}^n\frac{k_a\cdot k_b}{\sigma_{ab}}=0$ in $\rho$.

By simple power counting of poles we see that there is no pole and therefore no residue at infinity in $\rho$. As we deform the contour in $\rho$, the expressions $\{h\}$ change dynamically since they depend on $\rho$ directly and through $\xi(\rho)$ variables. When we localize $\rho$ 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 $\rho$ 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 $\rho$ as follows.

Structurally, the only type of poles that exists in the integrand is of the shape $1/\sigma_{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/\sigma_{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/\sigma_{ab}$ pole is localized due to the d.o.f. in $\rho$, the $\rho$ dependence in contributing $\{\xi(\rho)\}$ 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/\sigma_{ab}$ poles in the integrand, with $a\neq b$ pairs $a,b\in\{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/\sigma_{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 $\binom{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 $\rho$ and $\xi_i$ variables the Parke-Taylor like factor reads: $$\begin{aligned}
\label{ptlike}
\frac{\sigma_{n,m+1}}{\sigma_{n1}\sigma_{12}...\sigma_{m,m+1}}=\frac{\sigma_{n,m+1}}{\left(\sigma_{n}-\rho-\sum_{i=1}^{m-1}\frac{\xi_i}{2}\right)\xi_1\xi_2...\xi_{m-1}\left(\rho-\sum_{i=1}^{m-1}\frac{\xi_i}{2}-\sigma_{m+1}\right)}.\end{aligned}$$ Structurally, there are three different classes of $m$-poles combinations that can occur, namely where all appearing poles are localized except for: $$\begin{aligned}
\label{poletypes}
  1.)&~ \text{the pole }\frac{1}{\sigma_{n}-\rho-\sum_{i=1}^{m-1}\frac{\xi_i}{2}},\notag\\
  2.)&~ \text{xor a single pole }\frac{1}{\xi_i}\text{  out of }i\in\{1,2,...,m-1\},\\
  3.)&~ \text{xor the pole }\frac{1}{\rho-\sum_{i=1}^{m-1}\frac{\xi_i}{2}-\sigma_{m+1}}.\notag\end{aligned}$$ We choose to parametrize the $m$ localized poles in the above three cases by a parameter $\bar\rho\to0$ as follows: $$\begin{aligned}
  1.)&~ \rho=\bar\rho+\sigma_{m+1},~~~\xi_j=\bar\rho\,\bar\xi_j\text{ for all }j\in\{1,2,...,m-1\},\notag\\
  2.)&~ \rho=\bar\rho+\frac{1}{2}(\sigma_{m+1}+\sigma_{n}),~~~\xi_i=\sigma_{n}-\sigma_{m+1}+\bar\rho\,\bar\xi_i,~~~\text{ and }~~~\xi_j=\bar\rho\,\bar\xi_j\text{ for all }j\neq i,\\
  3.)&~ \rho=\bar\rho+\sigma_{n},~~~\xi_j=\bar\rho\,\bar\xi_j\text{ for all }j\in\{1,2,...,m-1\}.\notag\end{aligned}$$ The new variables $\bar\xi_i$ account for the original degrees of freedom of $\xi_i$ variables at leading order after localizing $\bar\rho\to0$. Note that in all three cases we have $d\rho=d\bar\rho$, and the one pole that is not localized always directly reduces to $1/\sigma_{n,m+1}$ under $\bar\rho\to0$, which cancels the numerator in (\[ptlike\]). In general, if we define [^4] $$\begin{aligned}
\label{barsigtraf}
\bar\sigma_n=\bar\sigma_0=\bar\sigma_{m+1}\equiv 0~~~\text{and}~~~\bar\sigma_q&=1-\sum_{a=1}^{q-1}\frac{\bar\xi_a}{2}+\sum_{b=q}^{m-1}\frac{\bar\xi_b}{2}~~~\text{for}~~~~q\in\{1,2,...,m\},\\
d\bar\xi_1\wedge d\bar\xi_2\wedge...\wedge d\bar\xi_{m-1}&=2 \,d\bar\sigma_1\wedge d\bar\sigma_2\wedge...\wedge d\bar\sigma_{m-1}\end{aligned}$$ then, for all possible pole combinations, the behavior of (\[ptlike\]) for $\bar\rho\to0$ is parametrized as $$\begin{aligned}
\frac{\sigma_{n,m+1}}{\sigma_{n1}\sigma_{12}...\sigma_{m,m+1}}= \frac{1}{\bar\rho^m\prod_{{a=1}\atop{a\neq r}}^{m+1}\bar\sigma_{a-1,a}}+O\left(\frac{1}{\bar\rho^{m-1}}\right),\end{aligned}$$ where index $r\in\{1,2,...,m+1\}$ labels which one of the $m+1$ poles in the denominator of (\[ptlike\]) is not being localized. Similarly, for all $m+1$ possible pole combinations we obtain $$\begin{aligned}
\frac{1}{\sum_{a=1}^m\sum_{b=m+1}^n\frac{k_a\cdot k_b}{\sigma_{ab}}}=\frac{1}{\bar\rho\left(\sum_{a=1}^{r-1}\frac{k_a\cdot k_n}{\bar\sigma_a}+\sum_{b=r}^{m}\frac{k_b\cdot k_{m+1}}{\bar\sigma_b}\right)}+O(\bar\rho^0),\end{aligned}$$ with the same index $r$. Depending on the particular value of $r$ we also get[^5] $$\begin{aligned}
\frac{1}{\det(H)}=\frac{\bar\rho^{2m-2}}{\det(H_r)}+O\left(\bar\rho^{2m-1}\right)~~~\text{and}~~~\text{Pf}(\psi)=\frac{1}{\bar\rho^m}\text{Pf}(\psi_r)+O\left(\frac{1}{\bar\rho^{m-1}}\right),\end{aligned}$$ where now $H_r$ and $\psi_r$ only contain terms supported on the localized poles appearing in the Parke-Taylor like factor (\[ptlike\]) 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 $\psi_r$ only depend on momenta $k_n^\mu,k_1^\mu,k_2^\mu,...,k_{m+1}^\mu$ and polarizations $\epsilon_1^\mu,\epsilon_2^\mu,...,\epsilon_{m}^\mu$, as expected.

Plugging the above findings into (\[scontour\]) and collecting the overall power of $\bar\rho$ we observe $$\begin{aligned}
S_m=\sum_{r=1}^{m+1}\sum_{{\{h_r\}=0}\atop{\text{solutions}}}\oint \frac{d\bar\rho}{2\pi i} \left(\frac{1}{\bar\rho}\frac{m}{\det(H_r)}\frac{1}{\sum_{a=1}^{r-1}\frac{k_a\cdot k_n}{\bar\sigma_a}+\sum_{b=r}^{m}\frac{k_b\cdot k_{m+1}}{\bar\sigma_b}}\frac{1}{\prod_{{c=1}\atop{c\neq r}}^{m+1}\bar\sigma_{c-1,c}}\text{Pf}(\psi_r)+O(\bar\rho^0)\right),\notag\end{aligned}$$ so that it is now trivial to compute the residues in $\bar\rho$, since for all $r$ we just have a single simple pole at $\bar\rho =0$. The result is $$\begin{aligned}
\label{preres}
S_m=\sum_{r=1}^{m+1}\sum_{{\{h_r\}=0}\atop{\text{solutions}}}\frac{m}{\det(H_r)}\frac{1}{\sum_{a=1}^{r-1}\frac{k_a\cdot k_n}{\bar\sigma_a}+\sum_{b=r}^{m}\frac{k_b\cdot k_{m+1}}{\bar\sigma_b}}\frac{1}{\prod_{{c=1}\atop{c\neq r}}^{m+1}\bar\sigma_{c-1,c}}\text{Pf}(\psi_r).\end{aligned}$$ Under closer inspection we note that the Pfaffian factorizes as $\text{Pf}(\psi_r)=\text{Pf}\left(\psi_{[1,r-1]}^{(n)}\right)\text{Pf}\left(\psi_{[r,m]}^{(m+1)}\right)$ with definitions (\[softps\]), again due to trivial factorization properties of Pfaffians of block matrices with some zero blocks.

In principle, (\[preres\]) 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] $$\begin{aligned}
\label{softfacres}
S_m=&\sum_{r=1}^{m+1}\int d\nu_r\frac{1}{\prod_{\substack{c=1\\c\neq r}}^{m+1}\bar\sigma_{c-1,c}}\text{Pf}\left(\psi_{[1,r-1]}^{(n)}\right)\text{Pf}\left(\psi_{[r,m]}^{(m+1)}\right),\\
\label{softmeasure}
d\nu_r\equiv&~\prod_{i=1}^{m-1}d\bar\sigma_i\prod_{q=1}^{m-1}\delta(h_{q,r})\frac{2m}{\sum_{a=1}^{r-1}\frac{k_a\cdot k_n}{\bar\sigma_a}+\sum_{b=r}^{m}\frac{k_b\cdot k_{m+1}}{\bar\sigma_b}},\end{aligned}$$ where, identifying $k_0^\mu\equiv k_n^\mu$ and keeping $\bar\sigma_0\equiv\bar\sigma_n=\bar\sigma_{m+1}= 0$ and $\bar\sigma_m=2-\bar\sigma_1$ in mind, we have $$\begin{aligned}
\label{softscat}
h_{q,r}=\sum_{a=q}^{q+1}\sum_{\substack{b=0\\b\neq a}}^{m+1}(-1)^{a-q}\frac{k_a\cdot k_b}{\bar\sigma_{ab}}\theta_{(r-a-\frac{1}{2})(r-b-\frac{1}{2})}, \end{aligned}$$ with $\theta_x\equiv\theta(x)$ being the Heaviside step function. We call the constraints $h_{q,r}=0$ the soft scattering equations. The $2(j-i+1)\times 2(j-i+1)$ matrix $\psi_{[i,j]}^{(w)}$ can be written explicitly as $$\begin{aligned}
\label{softps}
\psi_{[i,j]}^{(w)}=&\left({{A_{[i,j]}}\atop{C_{[i,j]}^{(w)}}}~{{-(C_{[i,j]}^{(w)})^T}\atop{B_{[i,j]}}}\right),~~~\text{with $(j-i+1)\times (j-i+1)$ sub-matrices}\\
A_{[i,j]}=&\left\{{{\frac{k_a\cdot k_b}{\bar\sigma_{ab}}}\atop{0}}{{;}\atop{;}}{{ a\neq b }\atop{a=b}}\right.~~,~~B_{[i,j]}=\left\{{{\frac{\epsilon_a\cdot \epsilon_b}{\bar\sigma_{ab}}}\atop{0}}{{;}\atop{;}}{{a\neq b}\atop{a=b}}\right.~~,~~C_{[i,j]}^{(w)}=\left\{\begin{matrix}
\frac{\epsilon_a\cdot k_b}{\bar\sigma_{ab}}&;\, a\neq b\\
-\frac{\epsilon_a\cdot k_w}{\bar\sigma_{a}}-\sum_{\substack{q=1\\q\neq a}}^n\frac{\epsilon_a\cdot k_q}{\bar\sigma_{aq}}&;\, a=b
\end{matrix}\right.,\notag\end{aligned}$$ and with indices in the range $a,b\in\{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 $\tau\to0$. However, since (\[softscat\]) admits different solutions of types $\bar\sigma_{a,b}=O(1)$ and $\bar\sigma_{a,b}=O(\tau)$, the integrations in (\[softfacres\]) have to be evaluated before the result can be systematically expanded to leading order in $\tau$.

Explicit examples and general pattern {#examples}
=====================================

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 (\[softscat\]) and no delta functions to integrate. The result is just directly given by the sum over $r$ in (\[softfacres\]):[^7] $$\begin{aligned}
\label{s1res}
S_1=2\frac{\epsilon_1\cdot k_2}{s_{12}}-2\frac{\epsilon_1\cdot k_n}{s_{1n}},\end{aligned}$$ which clearly is the correct Weinberg soft factor.[^8] We see that the soft factor is composed out of two pieces such as: $$\begin{aligned}
\label{s1piece}
P^{(1)}_{1,2}\equiv 2\frac{\epsilon_1\cdot k_2}{s_{12}}.\end{aligned}$$ Anticipating the structure of higher $m$-soft factors, we also define $$\begin{aligned}
\label{s0piece}
P^{(0)}_{m+1}=P^{(0)}_{n}\equiv 1.\end{aligned}$$ Using (\[s0piece\]) and (\[s1piece\]) we can structurally write the Weinberg soft factor (\[s1res\]) as $$\begin{aligned}
\label{s1res2}
S_1&=P^{(1)}_{1,2}P^{(0)}_{n}-P^{(0)}_{2}P^{(1)}_{1,n}.\end{aligned}$$ Based on this and further explicit results of this section, we propose in (\[smres\]) 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: $$\begin{aligned}
\label{momspin}
k_i\cdot k_j = \frac{1}{2}\langle ij\rangle [ji]&,&\epsilon^{+}_i\cdot k_j =\frac{[ij]\langle jr_i\rangle}{\sqrt{2}\langle r_ii\rangle}&,&\epsilon^{-}_i\cdot k_j =\frac{\langle ij\rangle[jr_i]}{\sqrt{2}[ir_i]},\\
\label{polspin}
\epsilon^{+}_i\cdot\epsilon^{-}_j =\frac{\langle jr_i\rangle[ir_j]}{[jr_j]\langle r_ii\rangle}&,&\epsilon^{+}_i\cdot\epsilon^{+}_j =\frac{\langle r_ir_j\rangle[ji]}{\langle r_ii \rangle\langle r_jj\rangle}&,&\epsilon^{-}_i\cdot\epsilon^{-}_j =\frac{\langle ij\rangle[r_jr_i]}{[ir_i][ jr_j]},\end{aligned}$$ 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: $$\begin{aligned}
S^+_1=\frac{\langle n2\rangle}{\langle n1\rangle\langle 12\rangle},\end{aligned}$$ 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 $\sqrt{2}$ in $S^+_1$ per usual spinor helicity convention.

Two-soft gluons factor $S_2$
----------------------------

For $m=2$, there is one soft scattering equation (\[softscat\]) for each $r$, and the number of solutions organizes as follows for the different solution types and different values of $r$: $$\begin{aligned}
\begin{array}{c|c|c|c}
\text{ solution type }&~r=1~&~r=2~&~r=3~\\ \hline
\bar\xi_1\sim O(1)    &1&1&1\\ \hline
\bar\xi_1\sim O(\tau) &1&0&1
\end{array}.\end{aligned}$$ Adding up the contributions of all $5$ solutions and expanding to leading order in $\tau$, we obtain the following expression for $S_2$: $$\begin{aligned}
\label{s2res}
S_2&=P^{(2)}_{1,2,3}P^{(0)}_{n}-P^{(1)}_{2,3}P^{(1)}_{1,n}+P^{(0)}_{3}P^{(2)}_{2,1,n}.\end{aligned}$$ This agrees with the generalization (\[smres\]) for $m=2$. The quantities $P^{(0)}_{i}$ and $P^{(1)}_{i,j}$ are given by (\[s0piece\]), (\[s1piece\]), and the new contribution of type $P^{(2)}_{i,j,l}$ reads:[^9] $$\begin{aligned}
\label{s2piece}
P^{(2)}_{1,2,3}=\frac{s_{13} \epsilon _1\cdot \epsilon _2}{s_{123} s_{12}}-\frac{s_{23} \epsilon _1\cdot \epsilon _2}{s_{123} s_{12}}-\frac{4 \epsilon _1\cdot
   k_3 \epsilon _2\cdot k_1}{s_{123} s_{12}}+\frac{4 \epsilon _1\cdot k_2 \epsilon _2\cdot k_3}{s_{123} s_{12}}+\frac{4 \epsilon _1\cdot k_3 \epsilon _2\cdot
   k_3}{s_{123} s_{23}}.\end{aligned}$$ Counting the powers of $k_1$ and $k_2$ we see that this expression diverges as $\tau^{-2}$, as we expect for the two-soft gluon factor. The result (\[s2res\]) is gauge independent and reduces to the gauge fixed result found in [@Volovich:2015yoa] when we select the gauge $\epsilon_2\cdot k_3=0,~\epsilon_1\cdot k_n=0$.

Restricting to four dimensions, converting to spinor helicity formalism by use of (\[momspin\]) and (\[polspin\]), and choosing appropriate reference spinors we get the following expression for the non-trivial helicity combination $(+-)$ after some simplification via Schouten identities: $$\begin{aligned}
S^{+-}_2=\frac{\langle n2\rangle}{\langle n1\rangle\langle 12\rangle}\frac{[13]}{[12][23]}\left(1+\frac{\langle n1\rangle[13]\langle 32\rangle}{s_{123}\langle n2\rangle}+\frac{[1n]\langle n2\rangle[23]}{s_{n12}[13]}\right),\end{aligned}$$ which naturally agrees with the result found in [@Volovich:2015yoa]. The trivial helicity combination $(++)$ reduces to the product of single soft factors $S^{++}_2={\scriptstyle\frac{\langle n3\rangle}{\langle n1\rangle\langle 12\rangle\langle 23\rangle}}$ as expected. Again, an overall factor of $(\sqrt{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 [@Dixon:2010ik] to generate explicit lower point amplitudes, we can form amplitude ratios that correspond to the soft factor in appropriate soft kinematics.[^10] Keeping in mind the overall powers of $\sqrt{2}$ that are suppressed in spinor helicity, we expect to find the following relation at leading order in $\tau$: $$\begin{aligned}
\label{numrel}
|S_m|=|\frac{(\sqrt{2})^{m}A_n(1,2,...,n)}{A_{n-m}(m+1,m+2,...,n)}|.\end{aligned}$$ Indeed, if we generate a numeric kinematic point where $k_1^\mu,k_2^\mu$ 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. $$\begin{aligned}
 |S^{++}_2|=|\frac{2 A_6(1^+,2^+,3^+,4^+,5^-,6^-)}{A_4(3^+,4^+,5^-,6^-)}|,~~~\text{or}~~~|S^{+-}_2|=|\frac{2 A_6(1^+,2^-,3^+,4^+,5^-,6^-)}{A_4(3^+,4^+,5^-,6^-)}|,\end{aligned}$$ 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 $\tau$.[^11] Naturally, ratios of more complicated amplitudes yield the same agreement.

Three-soft gluons factor $S_3$
------------------------------

For $m=3$, there are two soft scattering equations (\[softscat\]) for each $r$, and the number of solutions organizes as follows for the different solution types and different values of $r$: $$\begin{aligned}
\begin{array}{c|c|c|c|c}
\text{ solution type }&~r=1~&~r=2~&~r=3~&~r=4~\\ \hline
\bar\xi_1\sim\bar\xi_2\sim O(1)  &2&1&1&2\\ \hline
\bar\xi_i\sim O(1),~\bar\xi_j\sim O(\tau) &2&1&1&2\\ \hline
\bar\xi_1\sim\bar\xi_2\sim O(\tau) &2&0&0&2
\end{array},\end{aligned}$$ where we imply $i\neq j$ and $i,j\in\{1,2\}$. Adding up the contributions of all $16$ solutions and expanding to leading order in $\tau$, we obtain the following expression for $S_3$: $$\begin{aligned}
\label{s3res}
S_3&=P^{(3)}_{1,2,3,4}P^{(0)}_{n}-P^{(2)}_{2,3,4}P^{(1)}_{1,n}+P^{(1)}_{3,4}P^{(2)}_{2,1,n}-P^{(0)}_{4}P^{(3)}_{3,2,1,n}.\end{aligned}$$ This agrees with the generalization (\[smres\]) for $m=3$. As before, expressions of type $P^{(0)}_{i}$, $P^{(1)}_{i,j}$ and $P^{(2)}_{i,j,l}$ are given by (\[s0piece\]), (\[s1piece\]) and (\[s2piece\]), while the new contribution of type $P^{(3)}_{i,j,l,t}$ can still be analytically computed to be:[^12] $$\begin{aligned}
\label{s3piece}
P^{(3)}_{1,2,3,4}&=\frac{1}{s_{12}}(w_{312}-u_{312}-u_{213}-v_{312}-v_{213})+\frac{1}{s_{23}}(w_{231}-u_{231}-u_{132}-v_{231}-v_{132})+\notag\\
&+\left(\frac{1}{s_{12}}+\frac{1}{s_{23}}\right)(u_{123}+u_{321}+v_{123}+v_{321}-w_{123})+\frac{8 \epsilon _1\cdot k_4 \epsilon _2\cdot k_4 \epsilon _3\cdot
   k_4}{s_{34} s_{234} s_{1234}}+\\
    &+\frac{8 \epsilon _1\cdot k_4 \left(\epsilon _2\cdot k_3 \epsilon _3\cdot
   k_4- \epsilon _3\cdot k_2\epsilon _2\cdot k_4\right)}{s_{23}
   s_{234} s_{1234}}+\frac{8 \left(\epsilon _1\cdot k_2 \epsilon _2\cdot k_4-
   \epsilon _2\cdot k_1\epsilon _1\cdot k_4\right) \epsilon _3\cdot k_4}{s_{12} s_{34}
   s_{1234}}+\notag\\
&+\frac{2 \epsilon _1\cdot \epsilon _2 \epsilon_3\cdot k_4}{s_{12} s_{1234}}\left(\frac{2
   s_{13}}{s_{123}}+\frac{2
   s_{14}}{s_{34}}-\frac{s_{1234}}{s_{34}}\right)+\frac{4
    \epsilon _2\cdot \epsilon _3 \epsilon _1\cdot k_4}{s_{23}
   s_{1234}}\left(\frac{s_{13}}{s_{123}}-\frac{s_{34}}{s_{234
   }}\right)+\frac{4 \epsilon _3\cdot \epsilon _1 \epsilon _2\cdot k_4}{s_{123}
   s_{1234}},\notag\end{aligned}$$ where we used the abbreviations $$\begin{aligned}
u_{ijl}\equiv \frac{4 \epsilon _i\cdot k_j  \epsilon_j\cdot \epsilon _l}{s_{ijl}}\left(\frac{1}{3}-\frac{ s_{l4}}{s_{ijl4}}\right),~~v_{ijl}\equiv \frac{8 \epsilon _i\cdot k_j \epsilon _j\cdot k_l \epsilon _l\cdot
   k_4}{s_{ijl} s_{ijl4}},~~w_{ijl}\equiv \frac{8 \epsilon _i\cdot k_j \epsilon _j\cdot k_4 \epsilon _l\cdot
   k_j}{s_{ijl} s_{ijl4}}.\notag\end{aligned}$$ Counting the powers of $k_1,k_2$ and $k_3$ we see that this expression diverges as $\tau^{-3}$, as we expect for the three-soft gluon factor.

Again, we can use (\[momspin\]) and (\[polspin\]) 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: $$\begin{aligned}
S^{-+-}_{3}=&\frac{[n2]}{[n1][12]}\frac{\langle13\rangle}{\langle12\rangle\langle23\rangle}\frac{[24]}{[23][34]}\left(1-\left[\frac{\langle1n\rangle[n2]\langle23\rangle}{s_{n123}\langle13\rangle}+\frac{[2n]\langle n|k_1+k_3|2]\langle23\rangle[34]}{s_{123}s_{n123}[24]}\right.\right.\\
&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~+\left.\left.\frac{[n1]\langle13\rangle[32]}{s_{123}[n2]}+\frac{\langle1n\rangle[n2]\langle23\rangle~[23]\langle3n\rangle[n4]}{\langle13\rangle s_{n12}s_{n123}[24]}+{{  {\tiny{\arraycolsep=0.3\arraycolsep\ensuremath{\left\{\begin{matrix}n&\leftrightarrow&4\\1&\leftrightarrow&3\end{matrix}\right\}}}}}}\right]\right).\notag\end{aligned}$$ Similarly, the case $(+--)$ with an appropriate choice of reference spinors and after some simplification via Schouten identities yields $$\begin{aligned}
S^{+--}_{3}&=\frac{\langle n2\rangle}{\langle n1\rangle\langle12\rangle}\frac{[14]}{[12][23][34]}\left(1-\frac{\langle n1\rangle[14]\langle42\rangle}{s_{1234}\langle n2\rangle}-\frac{[1n]\langle n|k_2+k_3|4]}{s_{n123}[14]}-\frac{[1n]\langle n2\rangle[23]\langle3n\rangle[n4]}{s_{n12}s_{n123}[14]}\right.\notag\\
&\left.-\frac{\langle n1\rangle[1|k_2+k_3|4\rangle[43]\langle32\rangle}{s_{123}s_{1234}\langle n2\rangle}-\frac{s_{n1}[12]\langle23\rangle[34]}{s_{123}s_{n123}[14]}+\frac{\langle n1\rangle[13]\langle32\rangle~[1n]\langle n3\rangle[34]}{\langle n2\rangle s_{123}s_{n123}[14]}\right).\end{aligned}$$ The trivial helicity configuration $(+++)$ as expected reduces to $S^{+++}_3={\scriptstyle\frac{\langle n4\rangle}{\langle n1\rangle\langle 12\rangle\langle 23\rangle\langle 34\rangle}}$, 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, (\[numrel\]) is expected to hold. Making use of the GGT package [@Dixon:2010ik] 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^\mu,k_2^\mu$ and $k_3^\mu$ 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. $$\begin{aligned}
|S^{-++}_3|=|\frac{2^{3/2} A_7(1^-,2^+,3^+,4^+,5^+,6^-,7^-)}{A_4(4^+,5^+,6^-,7^-)}|,~~~|S^{+-+}_3|=|\frac{2^{3/2} A_7(1^+,2^-,3^+,4^+,5^+,6^-,7^-)}{A_4(4^+,5^+,6^-,7^-)}|,~~~\text{etc.}\notag\end{aligned}$$ hold to at least the first 10 digits, after which the first sub-leading correction in $\tau$ 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 (\[softscat\]) for each $r$, and the number of solutions organizes as follows for the different solution types and different values of $r$: $$\begin{aligned}
\begin{array}{c|c|c|c|c|c}
\text{ solution type }&~r=1~&~r=2~&~r=3~&~r=4~&~r=5~\\ \hline
\bar\xi_1\sim\bar\xi_2\sim\bar\xi_3\sim O(1)  &5&2&1&2&5\\ \hline
\bar\xi_i\sim\bar\xi_j\sim O(1),~\bar\xi_l\sim O(\tau) &8&2&2&2&8\\ \hline
\bar\xi_i\sim O(1),~\bar\xi_j\sim\bar\xi_l\sim O(\tau) &5&2&1&2&5\\ \hline
\bar\xi_1\sim\bar\xi_2\sim \bar\xi_3\sim O(\tau) &6&0&0&0&6
\end{array},\end{aligned}$$ where we imply $i\neq j,~i\neq l,~j\neq l$ and $i,j,l\in\{1,2,3\}$. With the generalization (\[smres\]) in mind, we expect that the contributions for cases $r=2,3,4$ can be constructed from previously determined quantities (\[s1piece\]), (\[s2piece\]) and (\[s3piece\]). That is easily verified numerically by obtaining and summing over explicit approximate solutions to the soft scattering equations (\[softscat\]) in some example kinematics. This confirms that the structure $$\begin{aligned}
\label{s4res}
S_4&=P^{(4)}_{1,2,3,4,5}P^{(0)}_{n}-P^{(3)}_{2,3,4,5}P^{(1)}_{1,n}+P^{(2)}_{3,4,5}P^{(2)}_{2,1,n}-P^{(1)}_{4,5}P^{(3)}_{3,2,1,n}+P^{(0)}_{5}P^{(4)}_{4,3,2,1,n}\end{aligned}$$ continues to hold. Trying to obtain $P^{(4)}_{1,2,3,4,5}$ for $r=1$ (and $r=5$) we discover that finding the $12$ solutions of the type $\bar\xi_1\sim\bar\xi_2\sim \bar\xi_3\sim O(\tau)$ is equivalent to solving for the roots of two $6$th 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 $\epsilon _1\cdot k_2 \epsilon _2\cdot k_3 \epsilon _3\cdot k_4 \epsilon _4\cdot k_5$, for which the correct contribution can be guessed (and numerically checked) to be given by: $$\begin{aligned}
16 \frac{ \epsilon _1\cdot k_2
   \epsilon _2\cdot k_3 \epsilon _3\cdot k_4 \epsilon _4\cdot k_5}{s_{1234} s_{12345}}\left(\left(\frac{1}{s_{12}}+\frac{1}{s_{23}}\right)\frac{1}{s_{123}}+\frac{1}{s_{12}
   s_{34}}+\left(\frac{1}{s_{23}}+\frac{1}{s_{34}}\right)\frac{1}{s_{234}}\right),\notag\end{aligned}$$ or terms like $\epsilon _1\cdot k_2 \epsilon _2\cdot k_3\epsilon _3\cdot \epsilon _4 $ with the correct guess for the contribution being: $$\begin{aligned}
8\frac{ 
     \epsilon
   _1\cdot k_2 \epsilon _2\cdot k_3\epsilon _3\cdot \epsilon _4}{s_{1234}}\left(\frac{1}{4}-\frac{s_{45}}{s_{12345}}\right)\left(\left(\frac{1}{s_{12}}+\frac{1}{s_{23}}\right)\frac{1}{s_{123}}+\frac{1}{s_{12}
   s_{34}}+\left(\frac{1}{s_{23}}+\frac{1}{s_{34}}\right)\frac{1}{s_{234}}\right).\notag\end{aligned}$$ However, $P^{(4)}_{1,2,3,4,5}$ also contains terms such as $\epsilon _3\cdot \epsilon _4 \epsilon _1\cdot k_2 \epsilon _2\cdot k_5$ or $\epsilon _1\cdot \epsilon _2 \epsilon _3\cdot \epsilon _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 (\[softfacres\]) is correct.

Using (\[momspin\]) and (\[polspin\]) to pass to spinor helicity formalism in four dimensions, (\[numrel\]) is again expected to hold. Therefore, we generate a numeric kinematic point such that $k_1^\mu,k_2^\mu,k_3^\mu$ and $k_4^\mu$ have soft entries of order $10^{-10}$ while the rest of the momenta have hard entries of order $10^0$. Now we can solve (\[softscat\]) numerically and obtain the numeric soft factor $S_4$ as a sum over all $64$ solutions. Subsequently, making use of the GGT package [@Dixon:2010ik], we can generate explicit amplitude ratios and observe that e.g. $$\begin{aligned}
|S^{-+++}_4|=|\frac{4 A_8(1^-,2^+,3^+,4^+,5^+,6^+,7^-,8^-)}{A_4(5^+,6^+,7^-,8^-)}|,~|S^{-+-+}_4|=|\frac{4 A_8(1^-,2^+,3^-,4^+,5^+,6^+,7^-,8^-)}{A_4(5^+,6^+,7^-,8^-)}|,~\text{etc.}\notag\end{aligned}$$ hold to at least the first 10 digits, after which the first sub-leading correction in $\tau$ becomes important. As before, ratios of more complicated amplitudes yield the same agreement.

For even higher $m$, the soft scattering equations (\[softscat\]) 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 (\[softfacres\]), (\[softscat\]) and (\[softps\]), valid to leading order in $\tau$.

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: $$\begin{aligned}
\label{smres}
S_m&=\sum_{r=1}^{m+1}(-1)^{r+1}P^{(m+1-r)}_{r,r+1,...,m,m+1}P^{(r-1)}_{r-1,r-2,...,1,n}.\end{aligned}$$ In essence, if all soft factors $S_a$ with $a<m$ for a fixed $m$ are known, then all contributions to $S_m$ with $1<r<m+1$ are constructed from the lower point results, while the summand[^13] $r=1$ equals the only previously unknown contribution $P^{(m)}_{1,2,...,m,m+1}$. In general we define $P^{(0)}_{m+1}=P^{(0)}_{n}\equiv 1$ and $$\begin{aligned}
\label{Pdef}
P^{(i)}_{1,2,...,i,i+1}=&\int d\nu_1\frac{1}{\prod_{c=2}^{i+1}\bar\sigma_{c-1,c}}\text{Pf}\left(\psi_{[1,i]}^{(i+1)}\right).\end{aligned}$$ In this sense, it suffices to evaluate only the $r=1$ summand of (\[softfacres\]) 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 (\[softfacres\]) into (\[smres\]) 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 $\rho$ are not convenient. This implies the necessity of a transformation along the lines of (\[deltatrafo\]) with a non-trivial Jacobian, which is not easily guessed. We leave a general proof of the conjecture (\[smres\]), (\[Pdef\]) to future work.

Multi-soft factors in other theories {#othertheories}
====================================

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 $$\begin{aligned}
\mathcal{C}\equiv\frac{1}{\sigma_{12}\sigma_{23}...\sigma_{n1}}\end{aligned}$$ 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 (\[psimat\]), the matrix $\Psi_{n+1,n+2,...,n+q}^{n+1,n+2,...,n+q}$ which is (\[psimat\]) with rows and columns $n+1,n+2,...,n+q$ dropped, and the matrix $$\begin{aligned}
\label{chimatrix}
\chi=\left\{{{\frac{\delta^{I_a,I_b}}{\sigma_{ab}}}\atop{0}}{{;}\atop{;}}{{a\neq b}\atop{a=b}}\right. ,\end{aligned}$$ where $I_a,I_b$ are some internal space indices for scalar fields involved in the scattering process [@Cachazo:2014xea]. 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 $\delta^{I_a,I_b}=1$.

Multi-soft factors in bi-adjoint scalar $\phi^3$ theory
-------------------------------------------------------

The CHY formula for tree level scattering in bi-adjoint scalar $\phi^3$ theory can be written as (\[ampl\]) [@Cachazo:2013hca] with $\mathcal{I}_n^{YM}$ replaced by $$\begin{aligned}
\mathcal{I}_n^{\phi^3}=\mathcal{C}^2.\end{aligned}$$ Starting with this integrand, the considerations in sections \[CHYsoft\] and \[factorization\] 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: $$\begin{aligned}
\label{SmScalar}
S^{\phi^3}_m=\sum_{r=1}^{m+1}\int d\nu_r \frac{1}{\prod_{\substack{c=1\\c\neq r}}^{m+1}(\bar\sigma_{c-1,c})^2},\end{aligned}$$ with $d\nu_r$ given in (\[softmeasure\]), and the identification $\sigma_0\equiv\sigma_n$. As in the gluon case, the soft scattering equations contained in $d\nu_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 $$\begin{aligned}
S^{\phi^3}_1=&\frac{1}{k_n\cdot k_1}+\frac{1}{k_1\cdot k_2},\\
S^{\phi^3}_2=&\frac{1}{k_1\cdot k_2}\left(\frac{1}{k_n\cdot (k_1+k_2)}+\frac{1}{(k_1+k_2)\cdot k_3}\right),\\
S^{\phi^3}_3=&\frac{2}{s_{123}}\left(\frac{1}{k_1\cdot k_2}+\frac{1}{k_2\cdot k_3}\right)\left(\frac{1}{k_n\cdot (k_1+k_2+k_3)}+\frac{1}{(k_1+k_2+k_3)\cdot k_4}\right).\end{aligned}$$ 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^{\phi^3}_m$ can be used to evaluate $S^{\phi^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 (\[ampl\]) with $\mathcal{I}_n^{YM}$ replaced by $$\begin{aligned}
\mathcal{I}_n^{YMS}=2\,\mathcal{C}\,\text{Pf}(\chi)\frac{(-1)^{i+j}}{\sigma_{ij}}\text{Pf}(\Psi_{i,j,n+1,n+2,...,n+q}^{i,j,n+1,n+2,...,n+q}),\end{aligned}$$ where matrix $\chi$ is $q\times q$ dimensional (\[chimatrix\]), and $1\leq i<j\leq n$ can be selected arbitrarily [@Cachazo:2014xea]. 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 \[CHYsoft\] and \[factorization\] 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] $$\begin{aligned}
\label{SmYMS}
S^{YMS}_m=\sum_{r=1}^{m+1}\int d\nu_r \frac{1}{\prod_{\substack{c=1\\c\neq r}}^{m+1}\bar\sigma_{c-1,c}}\text{Pf}(\chi_{[1,r-1]})\text{Pf}(\chi_{[r,m]})\text{Pf}(A_{[1,r-1]})\text{Pf}(A_{[r,m]}),\end{aligned}$$ with $d\nu_r$ given in (\[softmeasure\]), and the identification $\sigma_0\equiv\sigma_n$. The matrix $A_{[i,j]}$ was defined in (\[softps\]), and the matrix $\chi_{[i,j]}$ relates to $\chi$ in (\[chimatrix\]) the same as $A_{[i,j]}$ relates to $A$ in (\[psimat\]). As in the gluon case, the soft scattering equations contained in $d\nu_r$ can be explicitly solved for the cases $m=1,2,3$, with exactly the same solutions. However, since $\text{Pf}(\chi_{[i,j]})$ vanishes when $\chi_{[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 $$\begin{aligned}
S^{YMS}_2=&\frac{1}{2k_1\cdot k_2}\left(\frac{k_n\cdot (k_2-k_1)}{k_n\cdot (k_1+k_2)}+\frac{(k_1-k_2)\cdot k_3}{(k_1+k_2)\cdot k_3}\right).\end{aligned}$$ This agrees with the result in [@Cachazo:2015ksa]. As before, the general expression $S^{YMS}_m$ can be used to evaluate $S^{YMS}_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 (\[ampl\]) with $\mathcal{I}_n^{YM}$ replaced by $$\begin{aligned}
\mathcal{I}_n^{NLSM}=\mathcal{C}\frac{4}{(\sigma_{ij})^2}\text{Pf}(A_{i,j}^{i,j})^2,\end{aligned}$$ where $A_{i,j}^{i,j}$ is the matrix $A$ defined in (\[psimat\]) with rows and columns $i,j$ removed, and $1\leq i<j\leq n$ can be selected arbitrarily [@Cachazo:2014xea].\
Starting with this integrand, the considerations in sections \[CHYsoft\] and \[factorization\] go through in the same manner. The general expression for the $m$-soft factor with particles $1,2,...,m$ going soft amounts to: $$\begin{aligned}
\label{SmNLSM}
S^{NLSM}_m=\sum_{r=1}^{m+1}\int d\nu_r \frac{1}{\prod_{\substack{c=1\\c\neq r}}^{m+1}\bar\sigma_{c-1,c}}\text{Pf}(A_{[1,r-1]})^2\text{Pf}(A_{[r,m]})^2,\end{aligned}$$ with $d\nu_r$ given in (\[softmeasure\]), and the identification $\sigma_0\equiv\sigma_n$. The matrix $A_{[i,j]}$ was defined in (\[softps\]). As in the gluon case, the soft scattering equations contained in $d\nu_r$ can be explicitly solved for the cases $m=1,2,3$, with exactly the same solutions. However, since $\text{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 $$\begin{aligned}
S^{NLSM}_2=&\frac{1}{2}\left(\frac{k_n\cdot (k_2-k_1)}{k_n\cdot (k_1+k_2)}+\frac{(k_1-k_2)\cdot k_3}{(k_1+k_2)\cdot k_3}\right).\end{aligned}$$ This agrees with the result in [@Cachazo:2015ksa]. As before, the general expression $S^{NLSM}_m$ can be used to evaluate $S^{NLSM}_4$ and higher soft factors numerically. We tested the results numerically against amplitude ratios in CHY formulation and found agreement. Additionally, our $S^{NLSM}_4$ numerically agrees with the result found in [@Du:2016njc].[^16]

CSW recursion for multi-gluon soft-factors in four dimensions {#CSWrecur}
=============================================================

As an alternative to the construction rules presented in [@Georgiou:2015jfa], we can set up a CSW type recursion [@Cachazo:2004kj] for the $m$-soft factors in four dimensions as follows. We start with the amplitude $A^{(m)}(k_n^{+1},k_1^{h_1},...,k_m^{h_m},k_{m+1}^{+1})$, where $k_i^{h_i}$ denotes the external momentum of the $i$-th particle with helicity $h_i\in\{+1,-1\}$. Here we have cyclically rotated the $n$-th position to be the first, and suppressed all entries $k_j^{h_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: $$\begin{aligned}
\label{funS}
&\bold{S}\left(A^{(m)}(k_{q_1}^{h_{q_1}},k_{q_2}^{h_{q_2}},...,k_{q_l}^{h_{q_l}})\right)=\\
&=\left\{\begin{matrix}
 {\substack{\displaystyle\sum_{\nu=\pm 1}\sum_{i=1}^{l-1}\sum_{{\substack{j=i+1\\j-i<l-1}}}^{l}\bold{H}\left(A_{j-i+2}(k_{q_i}^{h_{q_i}},...,k_{q_j}^{h_{q_j}},k_{p(q_i,...,q_j)}^{-\nu})\right)\frac{1}{P^2_{q_i,...,q_j}}\times\\\displaystyle~~~~~~~~~~\times\bold{S}\left(A^{(m)}(k_{q_{1}}^{h_{q_{1}}},...,k_{q_{i-1}}^{h_{q_{i-1}}},k_{p(q_i,...,q_j)}^{+\nu},k_{q_{j+1}}^{h_{q_{j+1}}},...,k_{q_l}^{h_{q_l}})\right)}}&~~~;~~~\displaystyle\text{ if }\sum_{a=1}^lh_{q_a}<l,\\
~&~\\
\displaystyle A^{(m)}(k_{q_1}^{h_{q_1}},k_{q_2}^{h_{q_2}},...,k_{q_l}^{h_{q_l}})&~~~;~~~\displaystyle~\text{ otherwise, }~
\end{matrix}\right.\notag\end{aligned}$$ as well as, making use of $\mu(x)\equiv\text{mod}(x-1,l)+1$, the function: $$\begin{aligned}
\label{funH}
&\bold{H}\left(A_l(k_{q_1}^{h_{q_1}},k_{q_2}^{h_{q_2}},...,k_{q_l}^{h_{q_l}})\right)=\\
&=\left\{\begin{matrix}
 {\substack{\displaystyle\sum_{i=1}^{l}\sum_{{\substack{j=i+1}}}^{i+l-3}\bold{H}\left(A_{j-i+2}(k_{q_{\mu(i)}}^{h_{q_{\mu(i)}}},...,k_{q_{\mu(j)}}^{h_{q_{\mu(j)}}},k_{p(q_{\mu(i)},...,q_{\mu(j)})}^{-1})\right)\frac{1}{P^2_{q_{\mu(i)},...,q_{\mu(j)}}}\times\\\displaystyle~~~~~~~~~~\times\bold{H}\left(A_{l+i-j}(k_{q_{\mu(j+1)}}^{h_{q_{\mu(j+1)}}},...,k_{q_{\mu(l+i-1)}}^{h_{q_{\mu(l+i-1)}}},k_{p(q_{\mu(j+1)},...,q_{\mu(l+i-1)})}^{+1})\right)}}&~~;~~\displaystyle\text{ if }\sum_{a=1}^lh_{q_a}<l-4,\\
~&~\\
\displaystyle A_l(k_{q_1}^{h_{q_1}},k_{q_2}^{h_{q_2}},...,k_{q_l}^{h_{q_l}})&;~~~\,~\displaystyle~\text{ otherwise. }~
\end{matrix}\right.\notag\end{aligned}$$ We supplement the above functions with the following resolution properties: $$\begin{aligned}
p(i,...,j,p(a,...,b),u,...,v)=&p(i,...,j,a,...,b,u,...,v),\\
p(i,...,j,r,a,...,b,r,u,...,v)=&p(i,...,j,a,...,b,u,...,v),\\
P^2_{i,...,j,p(a,...,b),u,...,v}=&P^2_{i,...,j,a,...,b,u,...,v},\\
P^2_{i,...,j,r,a,...,b,r,u,...,v}=&P^2_{i,...,j,a,...,b,u,...,v},\end{aligned}$$ 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 (\[funS\]) and (\[funH\]) may contain summands that immediately vanish due to trivial helicity configurations of sub-amplitudes involved that enter the $\bold{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 (\[funS\]) and (\[funH\]) 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: $$\begin{aligned}
\bold{S}'\equiv\bold{S}\text{ with multiplicity of each overall summand set to }1,\end{aligned}$$ 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: $$\begin{aligned}
S_m=\bold{S}'\left(A^{(m)}(k_{n}^{+1},k_{1}^{h_{1}},...,k_{m}^{h_{m}},k_{q_{m+1}}^{+1})\right).\end{aligned}$$ Finally, to evaluate the soft factor we use the substitutions $$\begin{aligned}
\label{subhard}
&A^{(m)}(k_{q_1}^{+1},k_{q_2}^{+1},...,k_{q_l}^{+1})\to\frac{\langle n-1,n\rangle\langle n,m+1\rangle\langle m+1,m+2\rangle}{\langle n-1,q_1\rangle\left(\prod_{i=1}^{l-1}\langle q_{i},q_{i+1}\rangle\right)\langle q_l,m+2\rangle},\\
\label{subsoft}
&A_l(k^{+1}_{q_1},...,k_{q_{i-1}}^{+1},k_{q_i}^{-1},k_{q_{i+1}}^{+1},...,k_{q_{j-1}}^{+1},k_{q_j}^{-1},k_{q_{j+1}}^{+1},...,k^{+1}_{q_l})\to\frac{\langle q_i,q_j\rangle^4}{\langle q_l,q_1\rangle\prod_{i=1}^{l-1}\langle q_{i},q_{i+1}\rangle},\end{aligned}$$ where entries like $|p(i,...,j)\rangle$ are evaluated by the usual CSW replacement $P_{i,...,j}|X]$ with reference spinor $|X]$. Superficially, due to (\[subhard\]) it might seem that the soft factor depends on $(n-1)$-st and $(m+2)$-nd external momentum as well. However, just as in [@Georgiou:2015jfa], this dependence always cancels out upon the CSW replacement of the shifted spinors at leading order in $\tau$.

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 [@Dixon:2010ik], and found numerical agreement at leading order in $\tau$. 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 $\tau$ can be used to isolate leading terms only.

Four-soft gluons from BCFW {#foursoftBCFW}
==========================

Naturally, it is also possible to apply BCFW recursion relations [@Britto:2004ap] 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\rangle$ BCFW shift, so that $2\to\hat 2$ and $3\to\hat 3$ with $$\begin{aligned}
|\hat 2\rangle=|2\rangle~~~,~~~|\hat 2]=|2]+z|3]~~~,~~~|\hat 3\rangle=|3\rangle - z |2\rangle~~~,~~~|\hat 3]=|3].\end{aligned}$$ It is trivial to see that under this shift only the following four diagrams could possibly contribute to the leading soft factor with any helicity configuration: $$\begin{aligned}
S_{4,A}=&A_4(n,1,\hat 2,-\hat P_{n12})\frac{1}{s_{n12}}S_2(\hat P_{n12},\hat 3,4,5),\\
S_{4,B}=&A_3(1,\hat 2,-\hat P_{12})\frac{1}{s_{12}}S_3(n,\hat P_{12},\hat 3,4,5),\\
S_{4,C}=&A_4(-\hat P_{345},\hat 3,4,5)\frac{1}{s_{345}}S_2(n,1,\hat 2,\hat P_{345}),\\
S_{4,D}=&A_3(-\hat P_{34},\hat 3,4)\frac{1}{s_{34}}S_3(n,1,\hat 2,\hat P_{34},5),\end{aligned}$$ while the complete four-soft gluon factor is given by $$\begin{aligned}
S_4=S_{4,A}+S_{4,B}+S_{4,C}+S_{4,D}\end{aligned}$$ 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 constraints $\hat P^2_{\cdots}=0$ provide the following $z$ values to leading order in the soft limit:[^18] $$\begin{aligned}
z_A=\frac{-s_{n12}}{\langle 2n\rangle[n3]}~~~,~~~z_B=-\frac{[12]}{[13]}~~~,~~~z_C=\frac{s_{345}}{\langle25\rangle [53]}~~~,~~~z_D=\frac{\langle34\rangle}{\langle24\rangle}.\end{aligned}$$ 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: $$\begin{aligned}
S_{4,A}^{-+++}=&~\frac{[3\, n]^3 \langle 1\, n\rangle ^3 \langle 5\, n\rangle }{s_{n12} s_{n123} \langle 1\, 2\rangle  \langle 4\, 5\rangle  \langle n|k_{12}|3] \langle 4|k_{n123}k_{n1}|2\rangle },\\
S_{4,B}^{-+++}=&~\frac{[2\, 3]^3 \langle n\, 5\rangle }{s_{123} [1\, 2] \langle 4\, 5\rangle  \langle 4|k_{23}|1] \langle n|k_{12}|3]},\\
S_{4,C}^{-+++}=&~0,\\
S_{4,D}^{-+++}=&~\frac{\langle n\, 5\rangle  \langle 4|k_{23}|n]{}^3}{\langle 2\, 3\rangle  \langle 3\, 4\rangle  \langle 4\, 5\rangle  [n\, 1] \langle 4|k_{23}|1] \langle 4|k_{123}|n] \langle 2|k_{n1}k_{n123}|4\rangle }+\\
&+\frac{\langle 1\, 5\rangle ^3 [n\, 5]}{s_{12345} \langle 1\, 2\rangle  \langle 2\, 3\rangle  \langle 3\, 4\rangle  \langle 4\, 5\rangle  \langle 5|k_{1234}|n]}+\frac{\langle n\, 5\rangle  \langle 1|k_{234}|n]{}^3}{s_{1234} s_{\text{n1234}} \langle 1\, 2\rangle  \langle 2\, 3\rangle  \langle 3\, 4\rangle  \langle 4|k_{123}|n] \langle 5|k_{1234}|n]}.\notag\end{aligned}$$ 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: $$\begin{aligned}
&S_{4,A}^{+-++}=\frac{[3\, n]^3 \langle 2\, n\rangle ^4 \langle 5\, n\rangle }{s_{\text{n12}} s_{\text{n123}} \langle 1\, 2\rangle  \langle 4\, 5\rangle  \langle 1\, n\rangle  \langle n|k_{12}|3] \langle 4|k_{\text{n123}}k_{\text{n1}}|2\rangle },\\
&S_{4,B}^{+-++}=\frac{[1\, 3]^4 \langle n\, 5\rangle }{s_{123} [1\, 2] [2\, 3] \langle 4\, 5\rangle  \langle 4|k_{23}|1] \langle n|k_{12}|3]},\\
&S_{4,C}^{+-++}=~0,\\
&S_{4,D}^{+-++}=\frac{\langle 5\, n\rangle }{\langle 2\, 3\rangle  \langle 3\, 4\rangle  \langle 5|k_{234}|1]}\left(\frac{[1\, 5]^3 \langle 2\, 5\rangle ^4}{s_{12345} s_{2345} \langle 4\, 5\rangle  \langle 2|k_{345}k_{12345}|n\rangle }+\frac{\langle 2|k_{34}|1]{}^4}{s_{1234} s_{234} \langle 4|k_{23}|1] \langle 2|k_{34}k_{1234}|n\rangle }\right)\notag\\
&+\frac{\langle 2\, n\rangle ^3}{\langle 1\, 2\rangle  \langle 2\, 3\rangle  \langle 3\, 4\rangle  \langle 1\, n\rangle  \langle 2|k_{\text{n1}}k_{\text{n1234}}|5\rangle }\left(\frac{\langle 2\, 5\rangle ^3 [n\, 5]}{\langle 4\, 5\rangle  \langle 2|k_{345}k_{12345}|n\rangle }+\frac{\langle 2\, n\rangle  \langle 5\, n\rangle  \langle 2|k_{34}|n]{}^3}{s_{\text{n1234}} \langle 2|k_{\text{n1}}k_{\text{n123}}|4\rangle  \langle n|k_{1234}k_{34}|2\rangle }\right).\end{aligned}$$ 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.

#### Helicity configuration $(+--+)$:

$~$\
For the helicity configuration of soft gluons $(1^+,2^-,3^-,4^+)$ we find: $$\begin{aligned}
&S_{4,A}^{+--+}=\frac{\langle 2\, n\rangle ^3}{s_{\text{n12}} \langle 1\, 2\rangle  \langle 1\, n\rangle  \langle n|k_{12}|3] \langle 2|k_{\text{n1}}k_{\text{n1234}}|5\rangle }\left(\frac{[4\, n]^3 \langle 2\, n\rangle  \langle 5\, n\rangle }{s_{\text{n1234}} [3\, 4]}+\right.\\
&~~~~~~~~~~+\left.\frac{[5\, n] \langle 2|k_{\text{n1}}k_{\text{n123}}|5\rangle {}^3}{\langle 4\, 5\rangle  \langle 2|k_{\text{n1}}k_{\text{n123}}|4\rangle  \left(s_{345} [3\, n] \langle 2\, n\rangle +s_{\text{n12}} [3\, 5] \langle 2\, 5\rangle \right)}\right)\notag\\
&S_{4,B}^{+--+}=\frac{\langle n|k_{23}|1]{}^3 }{s_{123} [1\, 2] [2\, 3] \langle n|k_{12}|3] \left(\langle 5|k_{1234}|n] \langle
   n|k_{23}|1]-\langle 5\, n\rangle  [1|k_{23}k_{123}|n]\right)}\times\\
&~~~~~~~~~~\times\left(\frac{[4\, n]^3 \langle 5\, n\rangle  \langle n|k_{23}|1]}{s_{\text{n123}} s_{\text{n1234}} \langle n|k_{123}|4]}+\frac{[5\, n] \langle 5|k_{23}|1]{}^3}{\langle 4\, 5\rangle \langle 4|k_{23}|1] \left([4\, 5]  \langle 5\, n\rangle  \langle 4|k_{23}|1]+  \langle 5|k_{23}|1]  \langle n|k_{1234}|5]\right)}\right)\notag\\
&~~~~~~~~~~+\frac{1}{[1\, 2] [2\, 3] \langle 5|k_{234}|1]}\left(\frac{[1\, 4]^4 \langle 5\, n\rangle }{s_{1234} [3\, 4] \langle n|k_{123}|4]}+\right.\notag\\
&~~~~~~~~~~\left.+\frac{[1\, 5]^3 \langle n\, 5\rangle  \langle 5|k_{23}|1]{}^4}{s_{12345} \langle 4\, 5\rangle  \langle 4|k_{23}|1] [1|k_{2345}k_{45}|3] \left([4\, 5] \langle 5\, n\rangle  \langle 4|k_{23}|1]+\langle 5|k_{23}|1] \langle n|k_{1234}|5]\right)}\right),\notag\\
&S_{4,C}^{+--+}=\frac{[4\, 5]^3 \langle 2\, 5\rangle ^3 }{s_{345} [3\, 4] \langle 2|k_{34}|5] \langle 2|k_{345}k_{12345}|n\rangle }\left(\frac{[1\, 5]^3 \langle 2\, 5\rangle  \langle n\, 5\rangle }{s_{12345} s_{2345} [1|k_{2345}k_{45}|3]}+\right.\\
&~~~~~~~~~~\left.+\frac{[n\, 5] \langle 2\, n\rangle ^3}{\langle 1\, 2\rangle  \langle 1\,
   n\rangle  \left(s_{345} [3\, n] \langle 2\, n\rangle +s_{\text{n12}} [3\, 5] \langle 2\, 5\rangle \right)}\right),\notag\\
&S_{4,D}^{+--+}=\frac{\langle 2\, 3\rangle ^3}{s_{234} \langle 3\, 4\rangle  \langle 4|k_{23}|1] \langle n|k_{1234}|5]}\left(\frac{[n\, 5] \langle n|k_{234}|1]{}^3}{s_{1234} s_{\text{n1234}} \langle 2|k_{34}k_{1234}|n\rangle }+\frac{[1\, 5]^3 \langle n\, 5\rangle }{s_{12345} \langle 2|k_{34}|5]}\right)\\
&~~~~~~~~~~+\frac{\langle 2\, 3\rangle ^3 [n\, 5] \langle 2\, n\rangle ^3}{\langle 1\, 2\rangle  \langle 3\, 4\rangle  \langle 1\, n\rangle  \langle 2|k_{34}|5] \langle 4|k_{\text{n123}}k_{\text{n1}}|2\rangle  \langle 2|k_{34}k_{1234}|n\rangle }.\notag\end{aligned}$$ In all diagrams again 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: $$\begin{aligned}
&S_{4,A}^{--++}=\frac{\langle 1\, 2\rangle ^3 [3\, n]^3 \langle 5\, n\rangle }{s_{\text{n12}} s_{\text{n123}} \langle 4\, 5\rangle  \langle 1\, n\rangle  \langle n|k_{12}|3] \langle 4|k_{\text{n123}}k_{\text{n1}}|2\rangle },\\
&S_{4,B}^{--++}=\frac{1}{s_{123} [1\, 2] [2\, 3] \langle 4|k_{23}|1] \langle 5|k_{1234}|n]}\left(\frac{\langle n\, 5\rangle  [n|k_{1234}k_{12}|3]{}^3}{s_{1234} s_{\text{n1234}} \langle 4|k_{123}|n]}+\frac{[5\, n] \langle 5|k_{12}|3]{}^3}{s_{12345} \langle 4\, 5\rangle }\right)\\
&~~~~~~~~~~+\frac{[3\, n]^3 \langle 5\, n\rangle }{s_{\text{n123}} [1\, 2] [2\, 3] \langle 4\, 5\rangle  [1\, n] \langle 4|k_{123}|n]},\notag\\
&S_{4,C}^{--++}=~0,\\
&S_{4,D}^{--++}=\frac{[n\, 5]}{\langle 2\, 3\rangle  \langle 3\, 4\rangle  \langle 5|k_{234}|1]}\left(\frac{\langle 2\, 5\rangle ^3}{s_{2345} \langle 4\, 5\rangle  [1\, n]}+\frac{\langle 2|k_{34}k_{1234}|5\rangle {}^3}{s_{1234} s_{12345} s_{234} \langle 4|k_{23}|1] \langle 5|k_{1234}|n]}\right)\\
&~~~~~~~~~~+\frac{\langle n\, 5\rangle  \langle 2|k_{34}|n]{}^3}{s_{234} s_{\text{n1234}} \langle 2\, 3\rangle  \langle 3\, 4\rangle  [1\, n] \langle 4|k_{23}|1] \langle 5|k_{1234}|n]}.\notag\end{aligned}$$ 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.

#### Helicity configuration $(+-+-)$:

$~$\
For the helicity configuration of soft gluons $(1^+,2^-,3^+,4^-)$ we find: $$\begin{aligned}
&S_{4,A}^{+-+-}=\frac{[3\, n]^3 \langle 2\, n\rangle ^4}{s_{\text{n12}} \langle 1\, 2\rangle  \langle 1\, n\rangle  [5|k_{\text{n1234}}k_{\text{n12}}|3] \langle n|k_{12}|3]}\times\\
&~~~~~~~~~~\times\left(\frac{[3\, n] [5\, n] \langle 4\, n\rangle ^3}{s_{\text{n123}} s_{\text{n1234}} \langle 4|k_{\text{n123}}k_{\text{n1}}|2\rangle }+\frac{[3\, 5]^3 \langle n\, 5\rangle }{[3\, 4] [4\, 5] \left(s_{345} [3\, n] \langle 2\, n\rangle +s_{\text{n12}} [3\, 5] \langle 2\, 5\rangle \right)}\right),\notag\\
&S_{4,B}^{+-+-}=\frac{[1\, 3]^4 [3\, 5]^3 \langle n\, 5\rangle }{[1\, 2] [2\, 3] [3\, 4] [4\, 5] [3|k_{12}k_{1234}|5] [1|k_{2345}k_{45}|3] \langle n|k_{12}|3]}+\\
&~~~~~~~~~~+\frac{[1\, 3]^4}{s_{123} [1\, 2] [2\, 3] \langle 4|k_{23}|1] \langle n|k_{1234}|5]}\left(\frac{\langle 5\, n\rangle  \langle 4|k_{123}|5]{}^3}{s_{1234} s_{12345} [3|k_{12}k_{1234}|5]}+\frac{[5\, n] \langle 4\, n\rangle ^3}{s_{\text{n1234}} \langle n|k_{12}|3]}\right),\notag\\
&S_{4,C}^{+-+-}=\frac{[3\, 5]^4 \langle 2\, 5\rangle ^3}{s_{345} [3\, 4] [4\, 5] \langle 2|k_{34}|5] \langle n|k_{12345}k_{345}|2\rangle }\times\\
&~~~~~~~~~~\times\left(\frac{[1\, 5]^3 \langle 2\, 5\rangle  \langle 5\, n\rangle }{s_{12345} s_{2345} [1|k_{2345}k_{45}|3]}+\frac{[5\, n] \langle 2\, n\rangle ^3}{\langle 1\, 2\rangle  \langle 1\, n\rangle  \left(s_{345} [3\, n] \langle 2\, n\rangle +s_{\text{n12}} [3\, 5] \langle 2\, 5\rangle \right)}\right),\notag\\
&S_{4,D}^{+-+-}=\frac{\langle 2\, 4\rangle ^4 [n\, 5] \langle 2\, n\rangle ^3}{\langle 1\, 2\rangle  \langle 2\, 3\rangle  \langle 3\, 4\rangle  \langle 1\, n\rangle  \langle 2|k_{34}|5] \langle 2|k_{\text{n1}}k_{\text{n123}}|4\rangle \langle n|k_{1234}k_{34}|2\rangle }+\\
&~~~~~~~~~~+\frac{\langle 2\, 4\rangle ^4}{s_{234} \langle 2\, 3\rangle  \langle 3\, 4\rangle  \langle 4|k_{23}|1] \langle n|k_{1234}|5]}\left(\frac{[n\, 5] \langle n|k_{234}|1]{}^3}{s_{1234} s_{\text{n1234}} \langle 2|k_{34}k_{1234}|n\rangle }+\frac{[1\, 5]^3 \langle n\, 5\rangle }{s_{12345} \langle 2|k_{34}|5]}\right).\notag\end{aligned}$$ In all diagrams again only one helicity configuration is non-vanishing. We tested the above result numerically against amplitude ratios and found agreement.

The author would like to thank Anastasia Volovich for suggesting this topic, initial collaboration, and useful comments regarding the draft. The author also thanks Song He for interesting discussions and encouragement on pursuing the topic during IGST16 conference. This work is supported by Simons Investigator Award \#376208 of A. Volovich.

[99]{} F. E. Low, “Scattering of light of very low frequency by systems of spin 1/2,” Phys. Rev.  [**96**]{}, 1428 (1954). M. Gell-Mann and M. L. Goldberger, “Scattering of low-energy photons by particles of spin 1/2,” Phys. Rev.  [**96**]{}, 1433 (1954).

F. E. Low, “Bremsstrahlung of very low-energy quanta in elementary particle collisions,” Phys. Rev.  [**110**]{}, 974 (1958). S. Weinberg, “Photons and Gravitons in s Matrix Theory: Derivation of Charge Conservation and Equality of Gravitational and Inertial Mass,” Phys. Rev.  [**135**]{}, B1049 (1964). S. Weinberg, “Infrared photons and gravitons,” Phys. Rev.  [**140**]{}, B516 (1965).

T. H. Burnett and N. M. Kroll, “Extension of the low soft photon theorem,” Phys. Rev. Lett.  [**20**]{}, 86 (1968). D. J. Gross and R. Jackiw, “Low-Energy Theorem for Graviton Scattering,” Phys. Rev.  [**166**]{}, 1287 (1968). S. Saito, “Low-energy theorem for Compton scattering,” Phys. Rev.  [**184**]{}, 1894 (1969). J. S. Bell and R. Van Royen, “On the low-burnett-kroll theorem for soft-photon emission,” Nuovo Cim. A [**60**]{}, 62 (1969). V. Del Duca, “High-energy Bremsstrahlung Theorems for Soft Photons,” Nucl. Phys. B [**345**]{}, 369 (1990).

A. Strominger, “On BMS Invariance of Gravitational Scattering,” arXiv:1312.2229 \[hep-th\]. T. He, V. Lysov, P. Mitra and A. Strominger, “BMS supertranslations and Weinberg’s soft graviton theorem,” arXiv:1401.7026 \[hep-th\]. D. Kapec, V. Lysov, S. Pasterski and A. Strominger, “Semiclassical Virasoro Symmetry of the Quantum Gravity S-Matrix,” arXiv:1406.3312 \[hep-th\]. F. Cachazo and A. Strominger, “Evidence for a New Soft Graviton Theorem,” arXiv:1404.4091 \[hep-th\].

H. Bondi, M. G. J. van der Burg and A. W. K. Metzner, “Gravitational waves in general relativity. 7. Waves from axisymmetric isolated systems,” Proc. Roy. Soc. Lond. A [**269**]{}, 21 (1962). R. K. Sachs, “Gravitational waves in general relativity. 8. Waves in asymptotically flat space-times,” Proc. Roy. Soc. Lond. A [**270**]{}, 103 (1962). G. Barnich and C. Troessaert, “Symmetries of asymptotically flat 4 dimensional spacetimes at null infinity revisited,” Phys. Rev. Lett.  [**105**]{}, 111103 (2010) \[arXiv:0909.2617 \[gr-qc\]\]. G. Barnich and C. Troessaert, “BMS charge algebra,” JHEP [**1112**]{}, 105 (2011) \[arXiv:1106.0213 \[hep-th\]\]. G. Barnich and C. Troessaert, “Comments on holographic current algebras and asymptotically flat four dimensional spacetimes at null infinity,” JHEP [**1311**]{}, 003 (2013) \[arXiv:1309.0794 \[hep-th\]\]. F. Cachazo, P. Svrcek and E. Witten, “MHV vertices and tree amplitudes in gauge theory,” JHEP [**0409**]{}, 006 (2004) doi:10.1088/1126-6708/2004/09/006 \[hep-th/0403047\]. R. Britto, F. Cachazo and B. Feng, “New recursion relations for tree amplitudes of gluons,” Nucl. Phys. B [**715**]{}, 499 (2005) doi:10.1016/j.nuclphysb.2005.02.030 \[hep-th/0412308\]. F. Cachazo, S. He and E. Y. Yuan, “Scattering of Massless Particles in Arbitrary Dimensions,” Phys. Rev. Lett.  [**113**]{}, no. 17, 171601 (2014) doi:10.1103/PhysRevLett.113.171601 \[arXiv:1307.2199 \[hep-th\]\]. B. U. W. Schwab and A. Volovich, “Subleading soft theorem in arbitrary dimension from scattering equations,” arXiv:1404.7749 \[hep-th\]. N. Afkhami-Jeddi, “Soft Graviton Theorem in Arbitrary Dimensions,” arXiv:1405.3533 \[hep-th\]. Z. Bern, S. Davies and J. Nohle, “On Loop Corrections to Subleading Soft Behavior of Gluons and Gravitons,” Phys. Rev. D [**90**]{}, no. 8, 085015 (2014) doi:10.1103/PhysRevD.90.085015 \[arXiv:1405.1015 \[hep-th\]\]. S. He, Y. t. Huang and C. Wen, “Loop Corrections to Soft Theorems in Gauge Theories and Gravity,” JHEP [**1412**]{}, 115 (2014) doi:10.1007/JHEP12(2014)115 \[arXiv:1405.1410 \[hep-th\]\]. J. Broedel, M. de Leeuw, J. Plefka and M. Rosso, “Constraining subleading soft gluon and graviton theorems,” arXiv:1406.6574 \[hep-th\]. Z. Bern, S. Davies, P. Di Vecchia and J. Nohle, “Low-Energy Behavior of Gluons and Gravitons from Gauge Invariance,” arXiv:1406.6987 \[hep-th\]. E. Casali, “Soft sub-leading divergences in Yang-Mills amplitudes,” arXiv:1404.5551 \[hep-th\]. A. J. Larkoski, “Conformal Invariance of the Subleading Soft Theorem in Gauge Theory,” arXiv:1405.2346 \[hep-th\]. C. D. White, “Diagrammatic insights into next-to-soft corrections,” Phys. Lett. B [**737**]{}, 216 (2014) doi:10.1016/j.physletb.2014.08.041 \[arXiv:1406.7184 \[hep-th\]\]. M. Zlotnikov, JHEP [**1410**]{}, 148 (2014) doi:10.1007/JHEP10(2014)148 \[arXiv:1407.5936 \[hep-th\]\]. C. Kalousios and F. Rojas, JHEP [**1501**]{}, 107 (2015) doi:10.1007/JHEP01(2015)107 \[arXiv:1407.5982 \[hep-th\]\]. Y. J. Du, B. Feng, C. H. Fu and Y. Wang, “Note on Soft Graviton theorem by KLT Relation,” JHEP [**1411**]{}, 090 (2014) doi:10.1007/JHEP11(2014)090 \[arXiv:1408.4179 \[hep-th\]\]. A. E. Lipstein, “Soft Theorems from Conformal Field Theory,” JHEP [**1506**]{}, 166 (2015) doi:10.1007/JHEP06(2015)166 \[arXiv:1504.01364 \[hep-th\]\]. S. G. Avery and B. U. W. Schwab, “Burg-Metzner-Sachs symmetry, string theory, and soft theorems,” Phys. Rev. D [**93**]{}, 026003 (2016) doi:10.1103/PhysRevD.93.026003 \[arXiv:1506.05789 \[hep-th\]\]. Y. t. Huang and C. Wen, “Soft theorems from anomalous symmetries,” JHEP [**1512**]{}, 143 (2015) doi:10.1007/JHEP12(2015)143 \[arXiv:1509.07840 \[hep-th\]\]. M. Bianchi and A. L. Guerrieri, “On the soft limit of closed string amplitudes with massive states,” Nucl. Phys. B [**905**]{}, 188 (2016) doi:10.1016/j.nuclphysb.2016.02.005 \[arXiv:1512.00803 \[hep-th\]\]. S. G. Avery and B. U. W. Schwab, “Residual Local Supersymmetry and the Soft Gravitino,” Phys. Rev. Lett.  [**116**]{}, no. 17, 171601 (2016) doi:10.1103/PhysRevLett.116.171601 \[arXiv:1512.02657 \[hep-th\]\]. M. Bianchi and A. L. Guerrieri, “On the soft limit of tree-level string amplitudes,” arXiv:1601.03457 \[hep-th\]. P. Di Vecchia, R. Marotta and M. Mojaza, “Subsubleading soft theorems of gravitons and dilatons in the bosonic string,” JHEP [**1606**]{}, 054 (2016) doi:10.1007/JHEP06(2016)054 \[arXiv:1604.03355 \[hep-th\]\]. M. Campiglia and A. Laddha, “Asymptotic symmetries of gravity and soft theorems for massive particles,” JHEP [**1512**]{}, 094 (2015) doi:10.1007/JHEP12(2015)094 \[arXiv:1509.01406 \[hep-th\]\]. M. Campiglia and A. Laddha, “Sub-subleading soft gravitons: New symmetries of quantum gravity?,” Phys. Lett. B [**764**]{}, 218 (2017) doi:10.1016/j.physletb.2016.11.046 \[arXiv:1605.09094 \[gr-qc\]\]. M. Campiglia and A. Laddha, “Sub-subleading soft gravitons and large diffeomorphisms,” JHEP [**1701**]{}, 036 (2017) doi:10.1007/JHEP01(2017)036 \[arXiv:1608.00685 \[gr-qc\]\]. P. Di Vecchia, R. Marotta and M. Mojaza, “Soft behavior of a closed massless state in superstring and universality in the soft behavior of the dilaton,” JHEP [**1612**]{}, 020 (2016) doi:10.1007/JHEP12(2016)020 \[arXiv:1610.03481 \[hep-th\]\]. A. Luna, S. Melville, S. G. Naculich and C. D. White, “Next-to-soft corrections to high energy scattering in QCD and gravity,” JHEP [**1701**]{}, 052 (2017) doi:10.1007/JHEP01(2017)052 \[arXiv:1611.02172 \[hep-th\]\]. H. Elvang, C. R. T. Jones and S. G. Naculich, “Soft Photon and Graviton Theorems in Effective Field Theory,” arXiv:1611.07534 \[hep-th\]. A. Sen, “Soft Theorems in Superstring Theory,” arXiv:1702.03934 \[hep-th\]. A. Sen, “Subleading Soft Graviton Theorem for Loop Amplitudes,” arXiv:1703.00024 \[hep-th\].

A. Laddha and A. Sen, “Sub-subleading Soft Graviton Theorem in Generic Theories of Quantum Gravity,” arXiv:1706.00759 \[hep-th\]. S. Chakrabarti, S. P. Kashyap, B. Sahoo, A. Sen and M. Verma, “Subleading Soft Theorem for Multiple Soft Gravitons,” arXiv:1707.06803 \[hep-th\]. P. Mao and H. Ouyang, “Note on soft theorems and memories in even dimensions,” arXiv:1707.07118 \[hep-th\]. S. Weinberg, “Pion scattering lengths,” Phys. Rev. Lett.  [**17**]{}, 616 (1966). R. Jackiw and L. Soloviev, “Low-energy-theorem approach to single-particle singularities in the presence of massless bosons,” Phys. Rev.  [**173**]{}, 1485 (1968).

S. L. Adler, “Consistency conditions on the strong interactions implied by a partially conserved axial vector current,” Phys. Rev.  [**137**]{}, B1022 (1965).

N. Arkani-Hamed, F. Cachazo and J. Kaplan, “What is the Simplest Quantum Field Theory?,” JHEP [**1009**]{}, 016 (2010) \[arXiv:0808.1446 \[hep-th\]\]. W. M. Chen, Y. t. Huang and C. Wen, “From U(1) to E8: soft theorems in supergravity amplitudes,” JHEP [**1503**]{}, 150 (2015) \[arXiv:1412.1811 \[hep-th\]\].

W. M. Chen, Y. t. Huang and C. Wen, “New fermionic soft theorems,” arXiv:1412.1809 \[hep-th\]. F. Cachazo, S. He and E. Y. Yuan, “Scattering Equations and Matrices: From Einstein To Yang-Mills, DBI and NLSM,” JHEP [**1507**]{}, 149 (2015) doi:10.1007/JHEP07(2015)149 \[arXiv:1412.3479 \[hep-th\]\]. F. Cachazo, S. He and E. Y. Yuan, “New Double Soft Emission Theorems,” arXiv:1503.04816 \[hep-th\]. A. Volovich, C. Wen and M. Zlotnikov, “Double Soft Theorems in Gauge and String Theories,” JHEP [**1507**]{}, 095 (2015) doi:10.1007/JHEP07(2015)095 \[arXiv:1504.05559 \[hep-th\]\]. T. Klose, T. McLoughlin, D. Nandan, J. Plefka and G. Travaglini, “Double-Soft Limits of Gluons and Gravitons,” JHEP [**1507**]{}, 135 (2015) doi:10.1007/JHEP07(2015)135 \[arXiv:1504.05558 \[hep-th\]\]. Y. J. Du and H. Luo, “On single and double soft behaviors in NLSM,” JHEP [**1508**]{}, 058 (2015) doi:10.1007/JHEP08(2015)058 \[arXiv:1505.04411 \[hep-th\]\]. P. Di Vecchia, R. Marotta and M. Mojaza, “Double-soft behavior for scalars and gluons from string theory,” JHEP [**1512**]{}, 150 (2015) doi:10.1007/JHEP12(2015)150 \[arXiv:1507.00938 \[hep-th\]\]. B. U. W. Schwab, “Subleading Soft Factor for String Disk Amplitudes,” JHEP [**1408**]{}, 062 (2014) doi:10.1007/JHEP08(2014)062 \[arXiv:1406.4172 \[hep-th\]\]. M. Bianchi, S. He, Y. t. Huang and C. Wen, “More on Soft Theorems: Trees, Loops and Strings,” Phys. Rev. D [**92**]{}, no. 6, 065022 (2015) doi:10.1103/PhysRevD.92.065022 \[arXiv:1406.5155 \[hep-th\]\]. B. U. W. Schwab, “A Note on Soft Factors for Closed String Scattering,” JHEP [**1503**]{}, 140 (2015) doi:10.1007/JHEP03(2015)140 \[arXiv:1411.6661 \[hep-th\]\]. P. Di Vecchia, R. Marotta and M. Mojaza, “Soft theorem for the graviton, dilaton and the Kalb-Ramond field in the bosonic string,” JHEP [**1505**]{}, 137 (2015) doi:10.1007/JHEP05(2015)137 \[arXiv:1502.05258 \[hep-th\]\]. M. Bianchi and A. L. Guerrieri, “On the soft limit of open string disk amplitudes with massive states,” JHEP [**1509**]{}, 164 (2015) doi:10.1007/JHEP09(2015)164 \[arXiv:1505.05854 \[hep-th\]\]. A. L. Guerrieri, “Soft behavior of string amplitudes with external massive states,” Nuovo Cim. C [**39**]{}, no. 1, 221 (2016) doi:10.1393/ncc/i2016-16221-2 \[arXiv:1507.08829 \[hep-th\]\]. P. Di Vecchia, R. Marotta and M. Mojaza, “Soft Theorems from String Theory,” Fortsch. Phys.  [**64**]{}, 389 (2016) doi:10.1002/prop.201500068 \[arXiv:1511.04921 \[hep-th\]\]. I. Low, “Double Soft Theorems and Shift Symmetry in Nonlinear Sigma Models,” Phys. Rev. D [**93**]{}, no. 4, 045032 (2016) doi:10.1103/PhysRevD.93.045032 \[arXiv:1512.01232 \[hep-th\]\]. S. He, Z. Liu and J. B. Wu, “Scattering Equations, Twistor-string Formulas and Double-soft Limits in Four Dimensions,” JHEP [**1607**]{}, 060 (2016) doi:10.1007/JHEP07(2016)060 \[arXiv:1604.02834 \[hep-th\]\]. A. P. Saha, “Double Soft Theorem for Perturbative Gravity,” JHEP [**1609**]{}, 165 (2016) doi:10.1007/JHEP09(2016)165 \[arXiv:1607.02700 \[hep-th\]\]. A. P. Saha, “Double Soft Theorem for Perturbative Gravity II: Some Details on CHY Soft Limits,” arXiv:1702.02350 \[hep-th\].

S. He, Z. Liu and J. B. Wu, “Scattering Equations, Twistor-string Formulas and Double-soft Limits in Four Dimensions,” JHEP [**1607**]{}, 060 (2016) doi:10.1007/JHEP07(2016)060 \[arXiv:1604.02834 \[hep-th\]\]. P. Di Vecchia, R. Marotta and M. Mojaza, “Double-soft behavior of the dilaton of spontaneously broken conformal invariance,” arXiv:1705.06175 \[hep-th\].

G. Georgiou, “Multi-soft theorems in Gauge Theory from MHV Diagrams,” JHEP [**1508**]{}, 128 (2015) doi:10.1007/JHEP08(2015)128 \[arXiv:1505.08130 \[hep-th\]\]. L. J. Dixon, J. M. Henn, J. Plefka and T. Schuster, “All tree-level amplitudes in massless QCD,” JHEP [**1101**]{}, 035 (2011) doi:10.1007/JHEP01(2011)035 \[arXiv:1010.3991 \[hep-ph\]\].

D. Nandan, J. Plefka and W. Wormsbecher, “Collinear limits beyond the leading order from the scattering equations,” JHEP [**1702**]{}, 038 (2017) doi:10.1007/JHEP02(2017)038 \[arXiv:1608.04730 \[hep-th\]\]. T. McLoughlin and D. Nandan, “Multi-Soft gluon limits and extended current algebras at null-infinity,” arXiv:1610.03841 \[hep-th\]. Y. J. Du and H. Luo, “Leading order multi-soft behaviors of tree amplitudes in NLSM,” JHEP [**1703**]{}, 062 (2017) doi:10.1007/JHEP03(2017)062 \[arXiv:1611.07479 \[hep-th\]\].

[^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.

[^2]: To see this, make the substitution (\[taulim\]) and expand the Pfaffian along rows and/or columns $1,2,...,m,n+1,n+2,...,n+m$. Retain only leading summands under $\tau\to0$, keeping in mind that solutions with $\sigma_{ab}=O(\tau)$ or $\sigma_{ab}=O(1)$ for $a,b\in\{1,2,...,m\}$ are possible. Finally, reassemble the remaining coefficients into Pf$(\psi)$.

[^3]: This is the case since $\frac{k_a\cdot k_b}{\bar\sigma_{ab}}=0$ for generic momenta only has the solution $|\bar\sigma_{ab}|=\infty$, which is non-analytic, while in the case of at least two summands $\frac{k_a\cdot k_b}{\bar\sigma_{ab}}+\frac{k_c\cdot k_d}{\bar\sigma_{cd}}=0$ finite solutions for the $\bar\sigma_i$ exist such that analyticity is preserved.

[^4]: Note that only $m-1$ of the $\bar\sigma_q$ are now linearly independent since we have $\bar\sigma_{m}=2-\bar\sigma_1$.

[^5]: Recall that $H$ is the derivative matrix of scattering equations. This means it is composed of elements $k_i\cdot k_j/\sigma_{ij}^2$ and their sums. While $i,j\in\{1,2,...,n\}$ initially, localizing the poles from the Parke-Taylor like factor (\[ptlike\]) as described above removes all dependence on $\sigma_{m+1},...,\sigma_n$. This factorizes the scattering equations and their Jacobian from the remaining $(n-m)$-point amplitude.

[^6]: Note the convention $\text{Pf}\left(\psi_{[i,j]}^{(w)}\right)\equiv 1$ when $i>j$.

[^7]: Recall that we imply $\bar\sigma_m=2-\bar\sigma_1$, which for $m=1$ reduces to $\bar\sigma_1=1$.

[^8]: The $s_{ij}=(k_i+k_j)^2$ is the usual Mandelstam variable.

[^9]: Here, for brevity we use that $2(k_1+k_2)\cdot k_3\approx s_{123}$ at leading order in $\tau$.

[^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: $\lambda_i^1=\sqrt{k_i^0+k_i^3},~\lambda_i^2=(k_i^1+i k_i^2)/\sqrt{k_i^0+k_i^3}$ and $\tilde\lambda_i=(\lambda_i)^*$.

[^12]: Again, we use that $2(k_1+k_2+k_3)\cdot k_4\approx s_{1234}$ and similar at leading order in $\tau$ to keep notation short.

[^13]: Or alternatively the summand $r=m+1$, which is related by simple index exchange.

[^14]: There seems to be no obstruction to assuming that a similar pattern should appear for soft theorems e.g. in the other theories discussed below as well, where appropriate.

[^15]: Again, we introduce the convention $\text{Pf}(\chi_{[i,j]})=\text{Pf}(A_{[i,j]})\equiv 1$ when $i>j$.

[^16]: Note a typo in eq. (4.10) of [@Du:2016njc]: The numerator of last expression on the first line should involve $q_5\cdot k_1$ instead of $q_4\cdot k_1$.

[^17]: 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).

[^18]: We use the convention $s_{ij}=\langle ij\rangle[ji]$, which with our spinor contraction conventions ($\langle ij\rangle=\lambda_i^1\lambda_j^2-\lambda_i^2\lambda_j^1$ and $[ij]=\tilde\lambda_i^2\tilde\lambda_j^1-\tilde\lambda_i^1\tilde\lambda_j^2$) corresponds to $(+,-,-,-)$ Minkowski metric signature.
