Bootstrapping solutions of scattering equations

The scattering equations are a set of algebraic equations connecting the kinematic space of massless particles and the moduli space of Riemann spheres with marked points. We present an efficient method for solving the scattering equations based on the numerical algebraic geometry. The cornerstone of our method is the concept of the physical homotopy between different points in the kinematic space, which naturally induces a homotopy of the scattering equations. As a result, the solutions of the scattering equations with different points in the kinematic space can be tracked from each other. Finally, with the help of soft limits, all solutions can be bootstrapped from the known solution for the four-particle scattering.


Introduction
The scattering equations are a set of algebraic equations connecting the kinematic space K n of n massless particles spanned by linearly independent Lorentz invariants s ij := 2k i · k j and the moduli space of Riemann spheres with n marked points M 0,n [1][2][3][4][5][6], where the unknowns z a ∈ CP 1 denote punctures on the Riemann sphere, and the kinematical invariants s ab satisfy momentum conservation and on-shell conditions, i.e. b =a s ab = −s aa = 0. This system has a global SL(2, C) symmetry, and thus only n−3 out of the n equations are independent. It has been proven that the number of independent solutions to the scattering equations is (n−3)! [6,7].
In a new formalism developed by Cachazo, He and Yuan (CHY) [8][9][10][11], the tree-level S-matrix in massless field theories is expressed as a multiple integral over the M 0,n . The integral is fully localized to the zeroes of the scattering equations and can be written as a sum over residues where z (i) ≡ (z (i) 1 , . . . , z (i) n ) stands for the i th solution and det ′ Φ is the relevant Jacobian determinant (see e.g. [8,12] for its explicit expression). The scattering equations are theory-independent, while the function I n encodes dynamics of the specific theory. We do not show the precise form of I n for any theory, since this paper focuses mainly on the scattering equations.
The scattering equations are universal for all theories, and are fundamental objects in quantum field theory as well as string theory. They shed new light on the perturbative S-matrix and even the theory itself. Due to the universality of the scattering equations, the CHY formalism provides an elegant representation for exposing relations between different theories [8][9][10], such as the famous double-copy relation between amplitudes in Yang-Mills theory and Einstein gravity [37][38][39]. Very recently, it was observed that the scattering equations can be interpreted geometrically as a diffeomorphism from the worldsheet associahedron to the kinematic associahedron [40]. This gives new insight into the origin of the scattering equations and the CHY formalism [41].
Due to the importance, it is crucially important to solve the scattering equations. Notwithstanding efforts have been made to solve the scattering equations or evaluate the CHY formulas [6,7,[42][43][44][45][46][47][48][49][50][51][52], a good method is still missing. In this paper, we close this gap: we develop an efficient technique to solve the scattering equations based on the numerical algebraic geometry.
Homotopy continuation Let us give a brief introduction to homotopy continuation [53,54], which is the primary method in numerical algebraic geometry that we will use throughout this paper. In order to solve a system of equations p(z) = 0 with p ≡ (p 1 , . . . , p N ) and z ≡ (z 1 , . . . , z N ), the basic idea is to introduce a continuous deformation (homotopy) p(z) → p(z, t), t ∈ [0, 1], that connects the target system p(z, 1) = p(z) with a start system p(z, 0) = q(z) whose solutions z(0) are known. Then the solutions z(1) of the target system can be obtained from z(0) via smooth paths as the continuation parameter t varies from 0 to 1. To be explicit, constructing a differentiable homotopy p(z, t) and differentiating it with respect to t lead to a system of ordinary differential equations (ODEs) on z = z(t) as follows: Viewing this as a system of linear equations on dz i /dt, it can be transformed into the following standard form: where terms in parentheses should be understood as matrices. Providing the initial condition z(0), the desired solutions z(1) of the target system can be obtained by integrating the system of the ODEs (4). Usually, numerical algorithms for initial value problems [55] are applied to obtain an estimate for z(1). This approximated solution serves as the initial guess of the true solution, and are fed to the Newton method to further improve its precision [54]. The homotopy continuation method described above has been well-studied, in particular on polynomial systems, during the past decades. Therefore this technique can be straightforwardly applied to the system of the scattering equations, since it is equivalent to a system of polynomial equations as follows [7]: with where s A := i<j∈A s ij , and three punctures have been fixed as (z 1 , z 2 , z n ) → (0, 1, ∞) by SL(2, C) invariance. Following the homotopy continuation method, a frequently used homotopy is: The advantage of such construction is that the start system has (n−3)! known solutions and the number of solutions remains unchanged for any regular t. Although such a homotopy can be used to solve the scattering equations (5) in principle, with some experimentations, we found that it is highly inefficient. One reason is on the technical side, saying that the complexity of evaluating ODEs (3) corresponding to the polynomial system is too high. Another reason is that the initial system is significantly different from the target system, thus implies that a lot of steps are spent to reach the target system.
In this paper, we extend the homotopy continuation method to solve the fractional scattering equations (1) by establishing an appropriate homotopy.
Instead of constructing the homotopy for the system of the scattering equations directly, we propose the physical homotopy in the kinematic space, i.e. S → S t , where S is a point in K n . The momentum conservation and on-shell conditions hold for S t at any t. More explicitly, a simple construction is: wheres ij and s ij are two sets of Mandelstam variables belonging to the physical region of interest in K n . Clearly, as long as on-shell conditions and momentum conservation are satisfied fors ij and s ij , they are satisfied for s ij (t). We define the kinematic homotopy 1 as a oneparameter smooth path in the kinematic space K n , like (7). The physical kinematic homotopy connects different points in K n , and this may be used to establish the connection between the physics quantities evaluated at different points. The kinematic homotopy S t naturally induces a homotopy for the scattering equations Since the physical homotopy preserves on-shellness and momentum conservation, the system has exact (n−3)! solutions for any regular t. To proceed, let us use the SL(2, C) redundancy to fix three punctures, for example The last equation f n = 0 is then trivially satisfied [6]. Differentiating other equations with respect to t gives the following system of ODEs: A perfect property is that the matrix Φ(t) has exactly rank n−3 at any t [6]. This ensures that there is no singularity in our algorithm. To improve numerical stability, we retain all (n−1) equations except f n = 0 which is satisfied trivially, and employ matrix decomposition methods [55] to generate the standard form, like eq. (4). Therefore, once the solutions of the scattering equations fors ij is known, the solutions for s ij can be obtained by numerically integrating the ODEs. However, so far the start solutions (the solutions of the scattering equations for kinematical invariants s ij (0) = s ij ) are not readily available yet. We would like to emphasize that it is highly non-trivial to obtain the start solutions, in particular when the multiplicity n is large. In order to initiate our program, we develop an algorithm based on the properties of the scattering equations in some special kinematical regions as well as the homotopy continuation technique. This algorithm will be described in detail in the following.
We employ the homotopy (7) again, i.e. s ij (t) = (1−t)ŝ ij + ts ij . Here the kinematical invariantsŝ ij satisfŷ s 1i > 0,ŝ 2i > 0,ŝ ij > 0, i, j ∈ {3, . . . , n−1}, (11) which are referred to as the positive region denoted by K + n in [56]. A remarkable property is that all (n−3)! solutions of the scattering equations in K + n are real [56]. More interestingly, after using the gauge fixing condition given previously, all puncturs (z 3 , . . . , z n−1 ) live inside the interval (0, 1) and distinct from each other for each solution. It is clear that due to this feature, the scattering equations in K + n can be solved much more easily, compared to generic kinematic regions. As will be detailed below, all (n−3)! real solutions can be obtained using the homotopy continuation technique too. 2 Once these solutions are readily available, they will serve as start solutions, and we can use the homotopy (7) and integrate the system of the ODEs (9) to generate the solutions for general kinematicss ij . It is also worth stressing that we can encounter singularities if we still adopt the real contour for t from 0 to 1, since the starting and target points live in unphysical and physical regions of K n respectively. A solution to avoiding the singularities is to employ a complex contour for t. In our program, we choose a simple contour consisting of two line segments in the complex t plane: 0 → 0.5 + 0.5 i → 1.
Now the final task is to obtain all solutions to the scattering equations for one point in K + n . Inspired by the soft limit of the scattering equations, we propose the following homotopy 3 All the remaining kinematic invariants can be easily obtained via on-shell conditions and momentum conservation. Clearly, this homotopy preserves the "positivity" of the kinematic region K + n . Another remarkable property is that, in the limit t → 0 which defines the soft limit k n−1 → 0, the kinematic space of n particles is reduced to (n−1)-particle one which is still in positive region. In this limit, f n−1 (t) is invariant up to a factor t, i.e., while other equations become nothing but the system of scattering equations associated with (n−1) particles without the soft leg in K + n−1 . In order to solve the scattering equations in K + n−1 , we can use the inverse soft homotopy (12) recursively until the four-particle case, whose unique solution is known, i.e. z 3 = −s 12 /s 13 with gauge fixing (z 1 , z 2 , z 4 ) → (0, 1, ∞). The equation corresponding to the soft particlẽ f n−1 = 0 (referred to as the soft equation) is equivalent to a polynomial equation of degree n−3 in z n−1 . For each solution of the scattering equations for the (n−1)point system without the soft particle, the n−3 zeroes of the soft equationf n−1 (z n−1 ) = 0 are distributed in the n−3 sub-intervals of (0, 1), separated by z 3 , z 4 , · · · , z n−2 . Thus simple numerical techniques such as the bisection method can be applied to obtain all n−3 roots. For using the inverse soft homotopy (12) each time, the similar method can be used to solve the soft equation. Here it should be noted that the f n−1 (t) is always replaced bỹ f n−1 (t) when we employ the inverse soft homotopy (12). Finally, we can obtain all (n−3)! solutions to the scattering equations for one point in K + n . With the start solutions from solving the scattering equations in K + n , by integrating the corresponding differential equations given in (9), we can obtain the solutions to the scattering equations for one point in K n .
To summarise, we have proposed a homotopy continuation method to solve the scattering equations and given a workable framework in detail. As shown schematically below (superscript (s) stands for the soft limit), our method consist of two main steps.

− −→ K n
Step II (14) The first step is to obtain the start solutions, which consists of two substeps: First, solve the scattering equations in K + n by using the inverse soft homotopy (12) recursively. Then, with these solutions as start solutions, we can use the homotopy (7) to solve the scattering equations for one point in the realistic target region. As the next step, once we have all (n−3)! solutions to the scattering equations for one physically realistic point in the kinematic space, we can track these solutions to any point in the kinematic space using the homotopy (7). In the second step, the solutions of the start system can be continued to the target system much more easily, since they both live in the same physically realistic region.
The method presented above has been implemented into a C++ program. For the numerical integration of dif-ferential equations, we adopt the Runge-Kutta-Fehlberg 7(8)-th order method [57] provided by Odeint [58]; and for the numerical solution of linear equation system, we adopt the Householder QR decomposition with column pivoting provided by Eigen [59]. In obtaining the start solutions, the local accuracy is set to be 10 −15 , while in the second step, the local accuracy is set to be 10 −7 . In both steps, the Newton method is adopted to increase the precision to 10 −15 . The code is available at [https://github.com/zxrlha/sehomo].
We consider the randomly selected non-exceptional points in the phase space corresponding to 2 → n−2 scattering up to n = 13. All tests were performed on a Macintosh laptop with a 2.7 GHz processor. The results of the computation times are summarized in Table I. In the table, t n are the computation times for obtaining all ♯(n) = (n−3)! solutions, andt n ≡ t n /(n−3)! represents the average time for each solution, for solving the scattering equations with a set of prepared initial solutions in the physically realistic region of K n . That is to say, they correspond to the Step II shown in (14). Here we would also like to note that our algorithm for obtaining the start solutions (i.e. the Step I in (14)) works well. In this step, the time cost is dominated by tracking solutions from unphysical positive region to physically realistic region in K n , while solving the scattering equations in the positive region recursively is very fast. For example, it costs less than 30 minutes for n = 11 case. As a consequence of the Newton method, all solutions can be obtained with an accuracy of 10 −15 . We have also checked that all solutions are distinct each other, thus we can verify that no solution is missed.
We observed that the total time to obtain all solutions increases significantly as n increases, mainly due to a factorial increase in the number of solutions. On the other hand, the average time of obtaining one solution increases much more slowly, and it is still at O(ms) level even for n = 13. It is noteworthy that obtaining different solutions are completely independent, thus can be done in parallel.
We also found that the time costs are dominated by solving the differential equations. Therefore if higher precision on solutions are requested, only the last step, i.e. the Newton iterations should be performed within higher precision, which have only small impact on the total time cost.
In addition, due to the property of the algorithm, for two neighboring points in the phase space, clearly it will be much easier to obtain the solutions of the scattering equations from each other. Therefore, one could speed up the calculation through a book-keeping method: first the initial solutions are prepared at several typical kinematic points rather than only one point, and the closest point are adopted as the initial point when do actual calculation.
Lastly, let us make a comparison between methods in our paper and in ref. [50]. In four dimensions in the spinor-helicity formalism, the scattering equations can be decomposed into 'helicity sectors' and written in terms of two-component spinors with additional variables involved (see e.g. refs. [60][61][62][63]). In ref. [50], a method was introduced to solve the spinor-valued scattering equations proposed in ref. [61] and implemented in Mathematica. Overall, our algorithm is much faster than the one in ref. [50] for obtaining all (n−3)! solutions. Here we identify some significant differences as follows. As already pointed out, obtaining solutions is completely independent of each other in our algorithm. In contrast, in ref. [50] the solutions are obtained sequentially, and as more solutions obtained, finding the next solution becomes increasingly difficult. Consequently, we can easily obtain all solutions for high points (e.g. n = 13), while even for n = 10 it is quite challenging to solve the equations for all helicity sectors by the package in ref. [50].
Conclusion and outlook In this paper we have proposed the kinematic homotopy which connects different points in kinematic space. Such a homotopy always preserves momentum conservation and on-shell conditions. With the physical homotopy, we developed an efficient algorithm to generate all numerical solutions of the scattering equations. This opens a new window of opportunity for further explorations in various prospectives.
First of all, this powerful method allows us to solve the scattering equations with high accuracy and high efficiency in different contexts. It is interesting to investigate the properties of the scattering equations and the CHY formulas in various kinematical regions, such as collinear and multi-Regge limits. While the discussion above is limed at the tree level, our method can be simply generalized to solve the scattering equations at loop level, which have been derived from ambitwistor strings.
In practical terms, it allows one to develop a new framework to compute scattering amplitudes at tree and loop level. Once one obtains all solutions to the scat-tering equations, as a next step, it is straightforward to generate tree amplitudes or loop integrands by summing up the contributions from these solutions. For instance, since the scheme to extend the CHY formalism to loop level has been developed at least for gauge and gravity theories, this makes possible to compute the amplitudes in these theories up to the two loop order.
More interestingly, the kinematic homotopy developed in this paper has further significance beyond solving the scattering equations. It is intriguing that the kinematic homotopy may provide an avenue to study various physical quantities, such as scattering amplitudes and scattering forms [40,64,65], in the kinematic space directly.
Acknowledgments We are grateful to Claude Duhr, Song He, Fabio Maltoni, Jun-Bao Wu and Ellis Yuan for useful comments on the manuscript. We would also like to acknowledge the hospitality of Galileo Galilei Institute in Florence, ZL also acknowledges the hospitality of the CERN theory division in Geneva and ITP, CAS in Beijing. ZL is particularly grateful to Hongbao Zhang for his kind hospitality and generous support during a visit to Beijing Normal University. The work of ZL was supported by the "Fonds Spécial de Recherche" (FSR) of the UCLouvain. This work has received funding from the European Union's Horizon 2020 research and innovation programme as part of the Marie Sk lodowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104).