A Monte Carlo Approach to the 4D Scattering Equations

The scattering equation formalism is a general framework for calculation of amplitudes in theories of massless particles. We provide a detailed introduction to the 4D scattering equation framework accessible to non-experts, outline current difficulties solving the equations numerically, and explain how to overcome them with a Monte Carlo algorithm. With this submission we include {\tt treeamps4dJAF}, the first publicly available {\sc Mathematica} package for calculating amplitudes by solving the scattering equations, supporting MHV analytical and N$^{k-2}$MHV numerical computations. The package provides a powerful and flexible computational tool for calculating tree-level amplitudes in super Yang Mills theories, Einstein supergravity and conformal supergravity. We tabulate sets of numerical solutions up to 9 points in all MHV sectors and 12 points in the NHMV sector which can be used for fast evaluation of amplitudes.


Introduction
The calculation of scattering amplitudes in quantum field theories is both an essential ingredient in interpreting data from high energy physics experiment, and a powerful tool for probing the deepest questions of theoretical physics. The development of new mathematical techniques for calculating amplitudes has lead to important advances in both of these areas. In the 1980s the introduction of spinor helicity notation gave way to new simplified computations of four dimensional amplitudes which previously had seemed intractable, for example [1][2][3]. In recent years there have been many advances in the set of technical tools that exist for calculating perturbative scattering amplitudes. Key tree-level techniques include recursive methods for calculation of higher point amplitudes from lower point inputs, known as BCFW recursion [4,5], and the formulation of field theory amplitudes in terms of string worldsheet calculations, known as twistor-string theory [6][7][8].
The simplicity of the Parke-Taylor form of the Yang-Mills MHV amplitude [3] inspired a description in terms of 2D current algebra [9]. This idea was then generalised to N = 4 super Yang-Mills amplitudes [6][7][8], and was also found to calculate N = 4 conformal supergravity amplitudes [10]. A similar worldsheet formula for N = 8 supergravity was found in [11]. The spectra of these string models contains only field theory degrees of freedom. Cachazo, He and Yuan extended the formula of [11] to a framework for calculating scattering of particles in arbitrary dimensions for a wide variety of theories in terms of a unified set of scattering equations [12][13][14]. We will refer to these equations as 'general d scattering equations'. In the CHY framework, tree-level amplitudes for different theories of massless particles are supported on the solutions of these scattering equations, which were first discovered in the context of ordinary string theory in [15,16]. Worldsheet expressions of this form are related in a deep way to recursive formulae arising from BCFW [17][18][19][20].
The power of the spinor helicity formalism in four dimensions along with the generality of the scattering equation formalism were then combined to produce 4D ambitwistor string theory [21,22]. This model gives a worldsheet description for the tree-level S-matrices of Yang-Mills theories and Einstein gravity which are valid for any number of supersymmetries and are supported on refined scattering equations which are graded by helicity degree. We will refer to these equations as '4D scattering equations'. In recent work we extended this formalism to include N = 4 conformal supergravity [23].
The scattering equations are understood in most detail at tree level, and this paper will consider only tree level amplitudes. There are extensions to the general d scattering equations to calculate loop integrands [24][25][26][27]. The 4D scattering equations currently only support tree level computations, although some preliminary one loop expressions exist [18].
Formulae supported on scattering equations have been successful in representing and calculating abstract theoretical properties of amplitudes such as soft limits [28][29][30][31][32], collinear limits [33] and relations between the S-matrices of different theories [14,34]. Direct eval-uation of amplitudes by solving the scattering equations is difficult, and some approaches solve the integrations by different methods [35,36]. The equations have (n − 3)! solutions at n points, and it is likely that finding all of these solutions analytically for generic kinematics is not possible. Calculating amplitudes in this framework then becomes a primarily numerical problem, which has been addressed only for the general d equations. These equations can be reduced to a simplified polynomial form [37], which Mathematica's inbuilt NSolve algorithm can solve numerically up to 9 points on a standard laptop. CHY provide an algorithm for finding individual solutions at higher points [12], but there are difficulties finding all solutions in this way. Solutions which start out distinct at the start of the algorithm generically degenerate and the full set of solutions is not recovered at the end; details are discussed in this paper. The 4D equations do not currently have an equivalent simplified polynomial form and depend on a larger set of variables than the general d equations, and as such NSolve cannot solve them above 7 points. We overcome these difficulties with a Monte-Carlo solution finding algorithm which finds all solutions.
There exist published Mathematica packages for evaluation of tree level amplitudes using BCFW [38,39], but to date there have been no equivalent packages for calculating amplitudes using the scattering equations either analytically or numerically. [40] provide an algorithm for calculation of amplitudes numerically without directly solving the scattering equations, but no explicit implementation is given. Available Mathematica packages focus on Yang-Mills theory, and we provide the first explicit publicly available algorithms for calculating Einstein supergravity and N = 4 conformal supergravity amplitudes at tree level.
The key purpose of this paper is to present the Mathematica package treeamps4dJAF, which explicitly calculates tree level amplitudes by solving the scattering equations. The analytical framework required for understanding the details of the 4D scattering equation formalism is covered, and we provide details of how to find full sets of solutions for a given set of numerical momenta and MHV degree by Monte Carlo algorithm. The package supports analytical computations in the MHV sector, and numerical computations in general MHV sectors for a wide variety of theories.
We hope that this paper, along with its associated package, will provide a computational tool that can serve for a number of different purposes. It can be seen as an lead-in to the scattering equation formalism those new to the field, as well as a general introduction to calculating amplitudes in Mathematica. It can be used as a reference to check if new techniques for calculating amplitudes agree with previously known results for the theories currently supported, and a play box for discovering further theories that can be described in the 4D scattering equation formalism. We intend to update the package with integrands for new theories once they have been discovered.
The structure of the paper is as follows. We start with a review of the general d and 4D scattering equations and the integrands for currently supported theories in section 2. Section 3 gives an overview of the analytical tools necessary to calculate amplitudes from the 4D scattering equation integral. Numerical methods for solving the equations outside of the MHV sector by Monte Carlo algorithm and for extracting component amplitudes from a superamplitude are detailed in section 4. We conclude and discuss further directions in section 5. Appendix A introduces the Mathematica package treeamps4dJAF which gives a concrete implementation of all of the algorithms discussed in the paper and is included with the paper. The package's key functions are detailed in this appendix along with usage examples. Finally appendix B gives detailed proofs of various properties of the 4D scattering equations, including those used in section 3.

Review and Conventions
In this section we review relevant details of the general d and 4D scattering equation formalism. For an n point N k−2 MHV amplitude we define N to be the set of all of the particles, i.e. N := {1, ...n}. We will refer to the following equations as the 'general d scattering equations' [13] j∈N j =i where k i are a set of n-point null momenta obeying momentum conservation, and s i are points on the Riemann sphere. In the general d formalism, tree-level amplitudes can then be calculated as integrals of some integrand f (s) over delta functions enforcing these equations, This type of integral arises from the calculation of a correlation function on a string worldsheet, and hence we will refer to the s as worldsheet variables. The equations have the standard SL(2) symmetry arising from global conformal transformations on the worldsheet, which reduces the number of integration variables by three. To balance this, three delta functions are removed and a corresponding Jacobian factor added, as denoted by the notation i∈N . Any valid integrand f must be covariant under this SL(2) symmetry, along with any other symmetries inherited from the corresponding amplitude. The number of integrations matches the number of delta functions, and hence this integral simply instructs us to sum the integrand times the relevant Jacobian factor over the (n − 3)! solutions to the scattering equations. In this way, calculating tree-level scattering for massless particles in many theories has been reduced to solving a set of algebraic equations. Further details can be found in [12,13].
In four dimensions we have additional structure due to the factorization of the Lorentz group SO(3,1) = SL(2) L ×SL(2) R , and the corresponding MHV classification which splits amplitudes up into different sectors depending on the number of negative helicity particles scattered. These additional structures can be described in terms of the spinor-helicity formalism, reviewed in [41,42]. To clarify the notation, we define angle and square brackets such that |i ∈ SL(2) L and |i] ∈ SL(2) R describing the external momenta. Then k i = |i] i| := |i] ⊗ i| is the explicitly null external momentum on leg i. We form antisymmetric SL(2) invariant brackets as ij := det(|i |j ) and [ij] := det(|i] |j]).
We work with supersymmetry in an on-shell superspace via the Grassmann variables η A i , where A is an R-symmetry index and η A i are related to the supermomenta by q A i = |i η A i . Amplitudes in supersymmetric theories are most concisely described as a supersymmetry covariant superamplitude. This corresponds to scattering on-shell supermultiplets of the theory, and the superamplitude is a Grassmann expansion in terms of the η A i . The superamplitude will be a Grassmann function with some well-specified weight in η variables N G . Amplitudes where N G is not a multiple of N are zero, and from this we define the Grassmann degree k G := N G N of the amplitude. Component amplitudes for individual states can be extracted by integrating against N G relevant η variables. A review of these techniques can be found in [41,42], and the extension to non-maximal supersymmetry is discussed in detail in [43].
As an example consider N = 4 super Yang-Mills with supermultiplet Φ ± = g − η 1 η 2 η 3 η 4 + η I η J η K ψ − IJK + η I η J φ IJ + η I ψ + I + g + . The MHV superamplitude A(+ + −−) has k G = 2, and for maximal supersymmetry the positive and negative supermultiplet are the same and hence we can choose which legs are + and −. For non-maximal SUSY the choice of assignment of superfields to each leg results in different superamplitudes. To extract a mixed fermion-gluon four point component amplitude, we read off the relevant η from the supermultiplet, and integrate the superamplitude against these Grassmann variables; The 4D scattering equations are a refinement of the general d scattering equations which depend on the MHV degree as well as the number of particles [21], and calculate supersymmetry covariant expressions. They make use of spinor-helicity notation and are written in terms of the angle and square bracket spinors rather than momentum dot products. This introduces an extra variable t i into the system for each particle which is related to the little group scaling of the spinor variables. We will think of the t i as additional worldsheet variables, promoting each worldsheet coordinate from a point on the Riemann sphere to a point in C 2 . We will use different coordinates for the worldsheet variables as σ i = 1 Combining all of worldsheet variables together as a matrix σ ∈ C 2×n , we will work with minors of this matrix defined as The external particles are now grouped into two different sets L and R with L R = N to respect the MHV degree of the amplitude. For Yang-Mills theory and Einstein gravity, the left set will correspond to the set of negative helicity particles (or more generally, those particles that sit in the negative helicity super-multiplet of the theory), and the right set will correspond to the positive helicity particles (resp. positive supermultiplet). In these cases the MHV degree k of the amplitude will be the size of the left set, k := |L|. More generally we can define an integrand in terms of some abstract left set, and indeed for example conformal supergravity amplitudes can have negative helicity particles in the left or the right sets, which we explain shortly.
We will refer to the following equations as the '4D scattering equations' [21]. They are specified by an integer n and a left set L.
As in the general d case, we can then write amplitudes for many different theories as integrals of some integrand f over these scattering equation delta functions; We will also analyse the integrands for some specific theories. Integrands are currently known for Einstein gravity and Yang-Mills theory with any valid number of supersymmetries [21], and for graviton multiplets in N = 4 conformal supergravity [23]. The supersymmetry structure is encoded by a set of fermionic delta functions of the form δ N η l − r∈R ηr (lr) , which can be thought of as scattering equations for the Grassmann on-shell superspace variables. In this paper we include the fermionic delta functions in the integrand, rather than as part of the scattering equations. The integrands for these theories are where det is an instruction to remove one row and column from the matrices before taking the determinant. H is the Hodges matrix andH the dual Hodges matrix, defined as In conformal supergravity the left set L does not correspond to the set of negative helicity superfields, which we denote Φ − (with Φ + the positive helicity superfields; Φ − Φ + = N ). The left set L of the scattering equations has |L| = k G , where k G is the Grassmann degree of the superamplitude considered, and the choice of L does not affect the amplitude. The factors H, F,H andF in conformal supergravity are generalisations of the gravitational inverse soft factor outside of the MHV sector. The H factors are the diagonal elements of the Hodges matrices, and the F are given bỹ (2.7) As shown in [28], there are A(n, k) = n−3 k−2 solutions to the n-point N k−2 MHV scattering equations. n k are the Eulerian numbers [44], tabulated in Figure 1. The inductive proof of the number of solutions finds that A(n, k) can be broken down by taking either a left set particle soft to reduce to an amplitude of the form A(n − 1, k − 1) or a right set particle soft to reduce to A(n − 1, k). Each of these solutions then has a given multiplicity such that

Calculating Amplitudes with the 4D Scattering Equations
In this section we describe the analytical results necessary for solving the scattering equations to calculate amplitudes given a certain integrand. These techniques will then be used explicitly to find amplitudes numerically and analytically in the package treeamps4dJAF.
Detailed proofs are given in appendix B. As shown in equation (2.4), at n points there are 2n 4D scattering equations depending on 2n worldsheet σ variables. There is a GL(2) symmetry on the worldsheet which acts as inhomogenously on the left and right set of worldsheet coordinates for G ∈ GL(2) as Under this GL(2) any minor of the form (lr) remains invariant, and hence the scattering equations are invariant. Fixing this gauge symmetry leaves 2n − 4 remaining degrees of freedom. Generally in this work we will restrict to gauge transformations specified by two particle labels i, j which fix σ i = ( 1 0 ) and σ j = ( 0 1 ), and refer to this operation as 'gaugefixing particles i and j'. As the GL(2) transformations act inhomogeneously on the left and right set, we can only fix either two left or two right particles this way. The details of gauge-fixing and the symmetries of the equations are explained in appendix B.2.
The system now appears to be over specified, and we must remove four equations. Two spinor equations i and j either from the left set or from the right set can be reduced to a momentum conserving delta function on support of the other scattering equations, as proved in appendix B.3. We then refer to having 'deleted particles i and j', and the remaining equations are a well-specified set of 2n − 4 equations in 2n − 4 variables.
The number of integrations in equation (2) is the same as the number of delta functions, and hence the integrations instruct us to sum over all of the solutions of the scattering equations. Deleting equations l, l ∈ L we can calculate the Jacobian of the remaining equations to solve the delta function integrals as follows .
where J n L is the Jacobian of the scattering equations with respect to the sigma variables, and the superscript l, l refers to removing four rows and columns corresponding to the two particle labels from the matrix. Details of the Jacobian to the scattering equations are explained in appendix B.5. It is now a well-formulated problem to solve the scattering equations and sum a theory dependent integrand f over the full set of solutions to produce an amplitude.
Calculating MHV amplitudes is generally more simple than calculating N k−2 MHV amplitudes, and we find this simplicity to be reflected in the structure of the 4D scattering equations. Analytical solutions to the 4D scattering equations are not known for general MHV degree, but in the MHV sector we can construct analytical solutions. First consider the case where the left set is L = {1, 2}. Then the MHV equations become The most obvious choice of equations to remove in this case is the two left set equations which become the overall momentum conservation delta function as detailed in appendix B.3, and we recover a Jacobian of 12 2 . Similarly we gauge-fix particles 1 and 2 to the identity in the Grassmannian, and arrive at the following form for each delta function of the right-set equations, which are solved by a Schouten identity; The full MHV solution along with its minors and the associated expression for the Jacobian of the delta functions is 1r . (3.5) The Jacobian as calculated this way is in agreement with the calculation from appendix B.5, where we also consider the Jacobian for higher MHV degree. Finding analytical solutions for generic kinematics outside of the MHV sector is currently an unsolved problem, apart from at 6 points NMHV where the scattering equations have four solutions. These solutions are found for the general d equations for d = 4 in [45]. Abel's theorem states that there is no algebraic solution in terms of nth roots to a general polynomial equation of degree five or higher with arbitrary coefficients [46]. To find analytical N k−2 MHV solutions above 6 points NMHV some underlying structure would have to exist within the coefficients of the equations, otherwise general analytical solutions are excluded by Abel's theorem. Full sets of solutions can be found analytically for some specific choices of momenta [47].
Due to these difficulties in finding analytical solutions, calculating amplitudes by solving the scattering equations outside of the MHV sector is primarily a numerical problem, which we address in section 4.
One key strength of the scattering equation formalism is that once a full set of solutions are known, amplitudes can be calculated in any theory at relatively small computational cost. The relevant integrand is chosen, and the only necessary operation is to sum over the solutions. Solutions to the scattering equations are graded only by MHV degree and not by a specific choice of left set. This implies that there must exist some transformation on the worldsheet which can map an integrand supported on scattering equations for one left set into an integrand for a different left set of the same length. For example such a mapping will allow calculation of all 6 point NMHV gluon amplitudes in Yang Mills theory with only one solution to the scattering equations; eg.
The following is an explicit co-ordinate transformation on the worldsheet which swaps two particles between the left and right sets of the scattering equations. 1 We single out two legs l 0 ∈ L and r 0 ∈ R which which will be swapped.
Under this transformation we find that the scattering equation integrand transforms as where L has l 0 swapped with r 0 . We now have an explicit transformation for how to calculate a new integrand f for a swap of the choice of left sets for the scattering equations. Details of the transformation under this mapping are given in appendix B.4. Repeatedly applying this transformation can be used to reassign any left set.

Numerical Methods
The 4D scattering equations can be thought of as 2n − 4 equations with variables in the Grassmannian Gr(2, n), as explained in appendix B.2. The equations are parametrised by a 1 We thank Paul Heslop for suggesting this transformation set of spinors obeying momentum conservation. A solution to the scattering equations at n points N k−2 MHV will be a mapping from external data to n−3 k−2 points in Gr(2, n). In the MHV case we give this mapping analytically in section 3, but finding analytical solutions for k > 2 is complicated due to the combinatorially increasing number of solutions. A well-specified problem is to provide explicit numerical momenta, which will usually be randomly sampled, and to then solve the resulting equations numerically. CHY provide an inverse-soft type algorithm for finding individual numerical solutions to the general d equations [12], but there are difficulties in constructing the full set of solutions in this way which we discuss in section 4.1. We provide an explicit algorithm which given a set of numerical momenta samples random numerical points in Gr(2, n) to find solutions to the equations stochastically. Algorithms of this type are known as Monte Carlo algorithms. Monte Carlo methods in high energy physics are well studied [48,49] and their application to solving non-linear algebraic equations is straightforward.
With enough computing power and time, any non-linear system of equations can be solved by Monte Carlo algorithm. The two key questions to address are when to stop the algorithm, and what distribution to sample the initial guess points from. The scattering equations are well-suited for solution in this way because the number of solutions is known, which gives a clear stopping condition. We address the sampling question in section 4.2. Finding a set of solutions this way is stochastic and can take a long time, with time complexity now distributed as a random variable which depends on n and k. The expectation of the time complexity increases as n increases and as k moves towards n 2 . One advantage of the 4D formalism that makes it better suited for solution by Monte Carlo algorithm is that the (n − 3)! solutions are broken down into Eulerian numbers of solutions, tabulated in Figure 1. This means that the algorithm can stop after finding a smaller number of solutions than in general dimensions.
Once a full set of numerical solutions along with the corresponding momenta and left set are known for a given number of points and MHV degree, they can then be used to calculate amplitudes in different theories for a selection of different external states by substituting different integrands into the sum over solutions. Solutions up to 12 points NMHV and 9 points in all MHV sectors are currently accessible to the algorithms of treeamps4dJAF, and we tabulate full solution sets with rational external data for these cases in the accompanying data file SolutionLookupTable.csv.
Tree-level amplitudes are all rational functions of external momenta, and hence for rational numerical external data they will be a rational number. The solutions to the scattering equations are in general not rational numbers, but given a set of rational kinematics we can calculate to very high precision at relatively low computational cost via deterministic algorithm once all solutions are known. It is then possible rationalize to the closest rational number to give exact numerical results for the amplitude. We provide support for this type of calculation in treeamps4dJAF.

Difficulties with CHY's Inverse Soft Algorithm
One proposed algorithm to find numerical solutions to the general d equations is that of CHY, which takes one of the momenta soft with parameter to reduce the equations from n points down to n − 1 points [12]. The soft parameter is then reintroduced, and the soft equation at O( ) is solved for each of the (n − 4)! solutions to the n − 1 point equations. The solutions then have a multiplicity of n − 3, and these points are input back into the system with moving slowly up from 0 to 1. As this algorithm involves slowly bringing the soft parameter back to the full n point system, it is referred to as an inverse soft algorithm. (n − 3)! solutions to the n point equations will be found in this way, but there is no guarantee that all of these solutions will be distinct, and hence they will not necessarily cover the full solution space.
The inverse soft algorithm is based on an inductive argument for counting the number total number of solutions, and we can extend it to the 4D case using the analogous 4D solution counting argument, which is reviewed in section 2. We take one soft parameter for a left set particle 1 ∈ L so that |1 → |1 , and a further parameter˜ for a right set particle n ∈ R so that |n] →˜ |n]. The 4D scattering equations become (4.1) We can see that worldsheet variable for the particle in the soft limit decouples, and the equations reduce to n − 1 point N k−2 MHV equations when = 1,˜ → 0 and n − 1 point N k−3 MHV equations when˜ = 1, → 0. Evaluated on the solution to the lower point equations, the remaining equation for the particle that decoupled gives the multiplicity for each solution, as shown in equation (2.8).
As detailed in [12] this method is sufficient to produce individual solutions for a specific n and k, but we find difficulties when trying to construct all of the solutions in this way. In four dimensions we find that two of the different solutions constructed from a lower point amplitude can converge to the same higher point solution, as shown at 6 points NMHV in Figure 2. Hence the maximum number of solutions this algorithm can find is n−3 k−2 , and generically it does not find all of the solutions.
MHV and MHV solutions are known analytically, and the first non MHV case is at 6 points NHMV, with four solutions. The soft limit gives a 5 point MHV amplitude, and the˜ soft limit produces a 5 point MHV amplitude. The soft limit equations both have two solutions. Figure 2 describes the norm of the MHV solutions as they evolve from = 0 up to = 1 in blue, and the MHV solutions from˜ = 0 up to˜ = 1 in yellow.
We define a matrix norm on the solutions by taking the standard norm on C n , after flattening the 2 × 4 worldsheet matrix with gauge-fixed rows deleted down to C 8 . Lines crossing on the plot in this description is not sufficient for the solutions to be equal, and indeed where the lines cross for ∈ [0, 1] the solutions are not equal. At = 1 the solutions converge. Interestingly when we continue the algorithm to run for slightly larger than the 1 we see that the solutions separate again. In the case shown in Figure 2, both of the solutions from 5 points MHV collide separately with two individual solutions from 5 points MHV . In general for different randomly selected numerical momenta at 6 points NMHV we find zero, one or two pairs of solutions converging. The number of pairs that converge appears to be random based on the selection of random momenta. Solution crashing is a ubiquitous phenomenon; it occurs for nearly all choices and the difficulty is rather to find cases where the solutions do not crash than to find cases where they do.
These difficulties with the inverse soft algorithm inspire developing new methods, and we solve this problem using a Monte Carlo algorithm.

Monte Carlo Algorithm
The Monte Carlo equation solving algorithm in treeamps4dJAF is implemented via NSolveMonteCarlo. Many initial random points are sampled from a distribution described below, and chosen as the initial conditions for a FindRoot calculation. These initial calls to FindRoot run for many iterations, and stop after only one digit of precision is met. Most initial guess points will be far from a solution, and will not converge to 1 digit of precision by the specified number of iterations. Those which do not converge are discarded, and the ones that do converge go back into FindRoot up to a higher precision. These points are now solutions, which are compared with a list of all currently found solutions and duplicates are discarded. The algorithm stops when a suitable stopping condition is met, which may be after a specified amount of time or number of iterations, or when enough solutions are found. A pseudo-code for this algorithm is

Return the solutions end function
NSolveMonteCarlo is tailored to the 4D scattering equations in the function NSolveScatteringEquations4D. This specifies the stopping condition to be when all of the Eulerian numbers of solutions are found, and selects an appropriate distribution to sample the random points from. Figure 3 gives a statistical analysis of the time complexity of the algorithm 2 and of the distribution of solution points for all currently accessible n and k. Figure 4 expands on this analysis for 6 points NMHV, giving a histogram of the timings. Timings are positively skewed with a similar shape for other n and k. Based on this, the algorithm can currently handle at most around 500 total solutions. Accessing the next cases would require around 1000 solutions, at 10 points N 2 MHV and 13 points NMHV.   of the solution points is gauge dependent, and as we always gauge fix particles 1 and 2 it is not possible to make a direct comparison between data points with k = k and with k = n − k , even though these cases at first glance should be symmetrical.
We perform a statistical analysis at 6 points NMHV, as guided by our analytical understanding in the MHV sector. A solution at 6 points NMHV can be considered as a matrix in C 2×4 after removing gauge-fixed columns. We project this down to a vector of real numbers in R 16 , and hence collect 48 real numbers for each set of numerical momenta which we refer to as 'solution points'. We then run the algorithm for statistically many sets of random numerical momenta, and collect all of the solution points together into one dataset. A histogram of this data is plotted in Figure 5. The data are best fitted by a Cauchy distribution, for which the probability density function has the form (4. 2) The data are symmetrically distributed around 0, and hence we see that the location parameter x 0 = 0. This leaves only one remaining parameter γ, which is tabulated for some cases in Figure 3. Apriori we would expect γ be a function of n, k and the size of the uniform distribution from which the random momenta are sampled. We can use some intuition from the analytical result in the MHV sector to reduce this down to just γ = γ(n, k). Solution points in the MHV sector have a form such as Re ( 12 1r ), as shown in section 3. Hence the size of the uniform distribution from which the random events are sampled will divide out statistically, and will not affect the distribution of the solution points. This independence is guaranteed for higher MHV degrees as solutions has mass dimension 0, and spinors have mass dimension 1 2 . Performing a statistical analysis of numerical solutions to the scattering equations outside of the MHV sector also verifies these properties, for example as shown in Figure 3 where different values for γ are tabulated. Our data in Figure 3 show that γ is insensitive to k also. We find this independence of γ on n and k to be intriguing and counter-intuitive.
From this analysis, it seems clear that we should sample the initial points from the relevant Cauchy distribution. However, very many initial points must be sampled for each iteration of the algorithm, and sampling points from Cauchy distribution is significantly slower than sampling from a uniform distribution. We find in practice that it is more efficient to approximate this Cauchy distribution with a uniform distribution. To approximate a normal distribution with a uniform distribution, one could choose a symmetric range with a width which is a small multiple of the standard deviation, eg. U ([−2σ, 2σ]). This poses a problem with the Cauchy distribution as it has infinite standard deviation. We use instead the median absolute deviation (MAD), which is roughly equivalent to the γ parameter of the Cauchy distribution, and is tabulated for different n and k in Figure 3. To solve the equations for a given n and k, NSolveScatteringEquations4D uses the tabulated MAD values from Figure 3 to sample initial points from U([−2 MAD, 2 MAD]).
When finding solutions by Monte Carlo algorithm, a large percentage (∼ 3 4 ) of the solutions tend to be found comparatively quickly, and it can take a long time to find the remaining solutions. It is possible that sampling a certain percentage of the solutions from the approximated uniform distribution and then sampling the remaining solutions from the relevant Cauchy distribution could help to overcome this problem. We leave this for future work.

Efficient Component Amplitude Extraction
Extracting individual component amplitudes from a super amplitude written as an expansion in Grassmann parameters is a well specified and well understood operation, as explained in section 2. Explicitly evaluating these calculations on the computer is not as simple as the operation itself might suggest. The most naïve application is to expand out the Grassmann super amplitude in terms of each of its factors, and to then throw away any of the terms which are equal to zero either due to η 2 = 0 or because do not match the integration measure. This algorithm is prohibitive in terms of memory usage in the computer. For example a product of m factors each with a sum of m terms of Grassmann numbers will result in a total of m m terms when expanding out naively, nearly all of which are zero. This motivates finding a more efficient algorithm.
In the scattering equation formalism, the Grassmann delta functions are always written in the form l∈L δ N η l − r∈R ηr (lr) . We can see that integrals of this expression over a subset of the η variables of dimension kN is a special case of a Grassmann integral of the form I G := d m η m i=1 ( m j=1 A ij η j ); in our real use case the matrix A has some components given by worldsheet minors, and the rest zero. Grassmann integrals of this type can be evaluated in an especially neat analytical form, and we find that I G = det(A). Numerical computation of determinants is generally implemented with an algorithm of time complexity around O(m 3 ), and hence we can extract Grassmann components in a very efficient way using this formula. Grassmann integration of functions of this form then consists simply of assigning the correct values to the matrix A and calculating its determinant. Any product of Grassmann delta functions takes this form, and specifically the fermionic delta functions of the scattering equations, as well as many other standard representations for superamplitudes.
We can improve the RAM and time efficiency of this algorithm further for numerical computations by storing the matrix as a sparse array, where all elements are assumed to be zero unless they are specified (in our usage case as worldsheet minors). We can then hold the evaluation of the determinant until explicit numerical values are substituted in for the worldsheet minors, and then calculate the determinant of a numerical matrix. Calculation in this way is very efficient, saving the need to store large intermediate analytical expressions in the RAM.

Conclusions
In this paper we provide the necessary tools for solving the 4D scattering equations analytically in the MHV sector and numerically for higher MHV degree. We outline a Monte Carlo algorithm for solving the scattering equations numerically, and an algorithm for extracting component amplitudes from Grassmann delta functions efficiently. This allows computation of tree-level amplitudes in Yang-Mills theory and Einstein gravity with any number of super symmetries, and N = 4 conformal supergravity. Amplitudes in these theories can be calculated explicitly using the accompanying Mathematica package treeamps4dJAF, which implements the algorithms described in this paper. When integrands are conjectured for new theories it will now be straightforward to test and evaluate them using the package. We include some discussion of the key functions from the package in this paper, and in the supporting files we provide the Mathematica package along with example code and full documentation, and a lookup table with solutions up to 12 points NMHV and 9 points in all MHV sectors.
A natural question to ask is whether these techniques are applicable to the general d scattering equations. It should be simple to map a full set of solutions at n points from the 4D scattering equations to the general d scattering equations for d = 4, which would allow calculating amplitudes for the wider variety of theories supported by these equations in four dimensions. Solving the general d equations directly using this algorithm is more difficult as they are not refined by helicity degree and all (n − 3)! solutions must be found.
Mathematica's inbuilt algorithms can already solve the general d equations numerically up to 9 points with a deterministic algorithm, so it would be necessary to improve on the efficiency of the Monte Carlo algorithm if it is to become a viable method for this problem.
treeamps4dJAF provides only a basic implementation of the Monte Carlo equation solving algorithm, and there are a number of different ways that it could be improved. NSolveMonteCarlo does not support parallel computation of different calls to FindRoot, and does not use an optimal algorithm for deciding whether solutions are duplicate. The sampling method used where all initial points are taken from an approximated uniform distribution could be improved on. An updated algorithm with these changes in Mathematica should have significantly better behaved time complexity. Re-writing the core solution finding algorithm in a different programming language better suited to low-level numerical computation such as C++ could also increase efficiency dramatically. These changes should allow solution finding in the 4D case for higher n and k, and it is possible that they could allow solving the general d equations by Monte Carlo algorithm to become viable.
It is intriguing to find that the Cauchy distribution arises in the statistical analysis of the solutions to the 4D scattering equations, and especially to find that the parameter of the distribution is insensitive to changing n and k. It would be interesting to explore how the functional form of this distribution is related to the physics of the equations; the statistical analysis of this paper is dependent on choice of GL(2) gauge, but it would be possible to do a similar analysis of GL(2) invariants.
A clear limitation of the 4D scattering equations is that they currently only support tree-level calculations, and an obvious future direction is to investigate calculation of looplevel integrands using these methods. Loop-level scattering equations exist in general dimensions [50,51], and future work could implement numerical algorithms to evaluate loop-level integrands as sums over solutions to the general d equations.
The 4D scattering equations are not restricted only to calculating amplitudes, they also cover form factors [52] and possibly further structures in quantum field theory. They can also be used to calculate more physically relevant standard model amplitudes [53]. The tools we have provided in this paper should be directly relevant to calculation of form factors, and could play a role in finding new types of structures that can be supported on scattering equations.

Acknowledgements
I would like to thank Simon Badger, Jake Bourjaily, Paul Heslop and Dan Rutter for useful discussions, Tim Whitbread for help with numerical algorithms, Themis Botsas for guidance on statistical methods, Theresa Abl for bug testing the package, the community at mathematica.stackexchange.com for help with many Mathematica related questions, and especially Arthur Lipstein, without whose excellent supervision I wouldn't have been able to complete this work. I have been funded for this work by EPSRC PhD scholarship EP/L504762/1.

A. The Mathematica Package treeamps4dJAF
treeamps4dJAF provides a set of computational tools in Mathematica for analytical and numerical calculation of amplitudes at tree-level. The package's most high-end functions are tailored for calculations in the 4D scattering equation formalism using the techniques derived in this paper. For example the package provides functionality for calculating MHV amplitudes analytically with specified external states in the supermultiplets of Yang-Mills and Einstein gravity theories, and graviton multiplets in conformal supergravity. Moving out of the MHV sector the package also provides functionality for solving the scattering equations numerically and using these solutions to calculate amplitudes. At a more lowlevel analysis of the code, many functions and a framework are provided for analytical computations in general in the spinner helicity formalism for amplitudes, including for example functions for dealing with antisymmetric brackets and the Schouten identity, and for evaluating expressions in terms of momentum twistors numerically. Each function is detailed in the documentation of the package, and can be accessed using Information via ?Amplitude or ??Amplitude. Here we provide an overview of the package's key functions.
The packages reserves, clears and protects a number of default heads for usage with symbolic calculations. These can be accessed using GetHeads and SetHeads, and they are explained in Figure 6.

Name
Default Head Mathematical Notation "Angle" ang

A.1 Calculating Amplitudes
The most high-level functions of treeamps4dJAF are used for calculation and readable display of analytical and numerical amplitudes.
Amplitude is the package's main function. It is used to calculate tree-level analytical MHV amplitudes and numerical amplitudes in all MHV sectors by solving the 4D scattering equations. Analytical output of Amplitude will be a function of the external momenta in terms angle and square brackets and Grassmann η variables; see Figure 6 for details. The input for Amplitude is overloaded in a number of different ways. Analytical amplitudes are specified in the form; • Amplitude[theory, N, states] • Amplitude[theory, N, n, L] Numerical amplitudes can take any of the following forms; • Amplitude[theory, N, momenta, states] • Amplitude[theory, N, momenta, n, L] • Amplitude[theory, N, momenta, L solutions , solutions, states] • Amplitude[theory, N, momenta, L solutions , solutions, L] n and L specifications can be used to calculate Grassmann superamplitudes, and providing external states will perform the Grassmann integrations to give a component amplitude. If a set of numerical momenta matching the state specification are entered, Amplitude will attempt to solve the relevant scattering equations, and return either a numerical value or a superamplitude Grassmann expansion with numerical coefficients. On entering a set of numerical momenta along with left set L solutions and a full set of solutions solutions which match the number of points and MHV degree of the specified states, Amplitude will sum the solutions over the relevant integrand and return the numerical amplitude.

In[1]:=
Amplitude["YM",0,"g-","g-","g+","g+"]//AmplitudesDisplay ListTheories[] returns a list of all theories with number of supersymmetries currently supported by treeamps4dJAF. The package's core functions for working with the 4D scattering equations are ScatteringEquations4D[n, L, gaugefix, delete] returns the scattering equations at n points with left set L, with worldsheet variables specified by gaugefix fixed to the identity, and with equations given by delete removed. ScatteringEquations4D[momenta, L, gaugefix, delete] returns the equations as a function of worldsheet minors only, for the specified numerical momenta. If gaugefix and delete are not entered, ScatteringEquations4D returns the scattering equations at n points with left set L in terms of worldsheet minors.
SolveScatteringEquations4D[n, L, gaugefix ] returns the analytical MHV and MHV solutions to the scattering equations derived in section 3 at n points with left set L and particles specified by gaugefix gauge-fixed to the identity.  [1,2] ang [2,3] , s[4, 1] → ang [1,2] ang [1,4] , s[4, 2] → ang [1,2] ang [2,4]  Integrand returns a theory-dependent tree-level integrand which can be summed over solutions to the scattering equations to return an amplitude. The output of Integrand will be an analytical expression which is a function of the worldsheet and external momenta, and Grassmann η variables. The input for Integrand is overloaded in different ways, for example; • Integrand[theory, N, states] • Integrand[theory, N, n, L] The specifications for the states work the same as for Amplitude, which inherits its specifications from Integrand. ang [1,2]squ [3,4] cir [1,2]cir [3,4] In [2]:= Integrand["CG", 4, "h-", "h-", "h+", "h+"] ApplyScatteringSolutions[expression, momenta, L, solutions] applies numerical momenta to expression, and then sums expression and the relevant Jacobian over the given solutions to the scattering equations, for the number of points specified by the length of momenta and left set L.

A.5 File Input-Output for Sets of Solutions
Solving the 4D scattering equations using NSolveScatteringEquations4D is stochastic, and can take a long time for higher points and MHV degrees. treeamps4dJAF provides functions for storing and reading set of Solutions to and from a CSV file, along with a lookup table of solutions in SolutionLookupTable.csv which can be used without solving the equations.
ReadSolution[filename, n, L, solutionnumber ] looks for a solution specified by n and L in filename.csv, and returns the relevant data if found.

B. Analytical Details of the 4D Scattering Equations
In this appendix we provide detailed proofs of the analytical results from section 3, along with some further general results for the 4D scattering equations.

B.1 Recovering the General d Scattering Equations for d = 4
Solutions to the 4D equations are grouped into sets for the different N k−2 MHV sectors. The general d equations depend only on n and hence do not encode this grouping. As the 4D scattering equations contain more information than the general d scattering equations we have 4D ⇒ general d (for d = 4), which we prove in this section. A different argument is given in [37]. The proof also provides an explicit method for finding solutions to the general d equations using those from the 4D specific case. Integrands for the general d equations for d = 4 can also be mapped to integrands for the 4D equations [54].
First we prove a lemma which holds for the general d equations. Define a world-sheet dependent momentum P (s) := j∈N k i s−s i . Then we prove that P (s) 2 = 0 ⇔ the general d equations. Note that there are no second order poles in P (s) 2 as all of the external momenta k i are null.
where we use partial fractions and relabel indices to arrive at the result. P (s) is now written explicitly as a sum of its poles, and hence can only be zero if all of its residues are zero. Then we conclude that P (s) 2  where we use the same steps as in the previous calculation, and in the second last equality we use the 4D scattering equations and k i = |i] i|. Given that |λ(s)] λ(s)| is explicitly constructed as a null vector, we then see that the 4D scattering equations imply that P (s) 2 = 0, and hence by the lemma they imply the general d scattering equations. This proof also provides an explicit mapping from a solution to the 4D scattering equations to a solution to the general d equations. We map a point in the solution space of the 4D equations to a point in the n-punctured Riemann sphere by writing each column as σ i = t −1 in agreement with equation (49) of [45], up to choice of SL(2) fixing.

B.2 Symmetries, Little Group Scaling and Grassmanians
The scattering equations have a GL(2) symmetry which can be realised in different ways in terms of a worldsheet redefinition, or a worldsheet redefinition with a corresponding little group rescaling. The worldsheet GL(2) symmmetry in (3.1) is a combination of the standard SL(2) symmetry of global conformal transformations in the string worldsheet s variables, and a GL(1) transformation corresponding to a rescaling of the worldsheet t variables. Any function f (σ) which is integrated against the scattering equation delta functions must transform as f (σ) → f (σ)(det G) n−2k under (3.1) to balance out the transformation of the measure. All of the integrands for the theories considered in section 2 satisfy this transformation law, as enforced by their underlying 4D ambitwistor string models [21]. Before considering joint worldsheet and little group transformations, we first analyse the little group scaling of amplitudes supported on the scattering equations. First consider a general amplitude A n,L with some arbitrary integrand f (σ, |i , |i]), Perform a different little group scaling for each particle, such that |i → α i |i , |i] → α −1 i |i]. Then we see that where in the last equation we choose β = α det G to keep the equations invariant. The measure and delta functions combined transform such that We have seen how the little group transformation of the amplitude A comes from little group covariance of the integrand f . Any f which integrates to an amplitude must transform covariantly under the little group, and hence we need that f transforms as a scaling transformation for any combined little group transformation and worldsheet rescaling. We can then conclude that for any f which describes an amplitude there must exist some x and y real numbers such that Given f which transforms in this way, we can choose the little group scaling α to such that (det G) 3n−2k+x α 2n−4k+y = 1, and hence this transformation is a symmetry for any amplitude supported on the scattering equations. This GL(2) invariance is the GL(2) invariance of the Grassmannian Gr(2, n), and hence in this sense we can think of the solutions to the scattering equations as living in Gr(2, n).

B.3 Deleting Equations and Momentum Conservation
In this appendix we demonstrate how to remove four equations to give a momentum conserving delta function. There are three possible cases; we either consider the equations for two particles in the left set, or for two in the right set, or for one in each. First consider the case of two particles in the left set. Without loss of generality, label these particles to be 1 and 2. Defining the delta functions for these particles to be ∆ 1,2 we see that ∆ We would now say that we have 'deleted equations 1 and 2', and the remaining 2n − 4 equations now give a well-specified system. Note the Jacobian 12 2 for this calculation. The calculation for deleting two equations in the right set goes by the same steps, and labelling the two particles in the right set to be r 1 and r 2 we find that There is one remaining choice; we could delete one equation from the left set and one equation from the right set. Choosing equations in this way does not produce a momentum-conservation delta function, and hence does not result in a solvable system for spinors satisfying momentum conservation. Label the left set particle as 1 ∈ L, and the right set particle as n ∈ R. Then following the same analysis as above, we arrive at (B.16) These equations looks deceptively similar to momentum conservation. We prove here that they are in fact not the same, and give a constraint corresponding to non-generic kinematics. Solving the first equation as |1] = 1 n1 n−1 i=2 |i] in and substituting into the second, we find that these delta functions imply that where between the first two equalities we split the second sum into two terms, relabel the indices use a Schouten identity. Hence to keep consistency with these equations we must set either |n or the Mandelstam invariant to zero, and neither of these choices correspond to generic kinematics. From this we see that it is not possible to delete one equation from each set.

B.4 Permutations and Choice of Left Set
A given n point N k−2 MHV amplitude will be supported on scattering equations with a specified left set L. In this section we show how solutions to the scattering equations for one left set can be used to calculate amplitudes with the same MHV degree that are supported on different left set L . This mapping of different left sets is important computationally as it will allow us to solve the equations only once in each MHV sector and then calculate all amplitudes of this MHV degree for the specified numerical momenta. The explicit worldsheet transformation swapping particles l 0 ∈ L and r 0 ∈ R in between the sets is given in equation (3.6). Under this transformation, we find that the scattering equations transform as where we used a Schouten identity in the worldsheet variables for the l and r equations. The second terms in the l and r equations are zero on support of the l 0 and r 0 equations, and hence we see that the scattering equations remain the same up to changing the particles l 0 and r 0 between the left and right set. The delta functions pick up an overall factor of (l 0 r 0 ) −4 . Calculating the transformation of the measure, we find that the scattering equation integral transforms as in equation (3.7).
If we now also perform a permutation of the external data on the same legs l 0 and r 0 and look at the transformation properties of the integrands, we can understand how the amplitude does not depend on a choice of left set for maximally supersymmetric theories and for N = 4 conformal supergravity. This transformation can also be used to understand permutation invariance under swapping between the right and left sets for N = 8 supergravity. Showing permutation invariance of gravity amplitudes under a swapping two legs which carry the same helicity superfield is straightforward, and simply requires renaming the worldsheet variables on the permuted legs.
Let us now look at the form of the integrand for Yang-Mills theories. Under this transformation, we find that d 2×n σ GL (2) δ 2×n (SE n L ) 2r 0 4 , which is exactly the factor required to modify the Park-Taylor formula for particles 1 and 2 negative helicity gluons to have particles 2 and r 0 with negative helicities. It is interesting to note that this structure extends outside of the MHV sector at the level of the integrand in the scattering equation formalism.

B.5 The Jacobian of the Scattering Equations
In this appendix we detail some properties of the Jacobian of the 4D scattering equations. For a general worldsheet integral over the scattering equation delta functions, we can solve the integrals and write the expression as a sum over the solutions to the scattering equations. To do this we need an explicit expression for the Jacobian J n L (σ) as follows. We use the notation for matrix A that A ij has rows and columns i and j removed. Let l, l ∈ L and delete equations l and l to arrive at d 2×n σ GL (2) δ 2×n (SE n L ) f (σ) = δ 4 (P ) σ sol ∈solutions f (σ sol ) ll −2 det(J n ll L (σ sol )) .

(B.19)
Note that as shown in appendix B.3, the two rows/columns deleted must either both be in the left set or both in the right set, and there is an extra associated factor eg. ll 2 such that the full determinant of the Jacobian is ll −2 det(J n ll L ) for two left set particles deleted and gauge-fixed.
At n points with left set L we find the Jacobian to be where the matrix is written in a block form with blocks of sizes of the left and right set, and each element of these matrices is broken down into a 2 × 2 matrix which is a tensor product of a spinor with a worldsheet vector. This matrix has determinant 0, which is simple to check analytically for example in Mathematica. This is insured by the fact that gauge-fixing removes 4 of the sigma variables, and hence we have to remove four of the rows of the matrix to produce a well specified system. Similarly we must remove four of the columns of J n L , which is equivalent to deleting two equations as shown in appendix B.3. In the MHV sector this matrix is block diagonal and we can calculate the determinant in terms worldsheet minors directly. Taking the left set to the particles 1 and 2 and also gauge-fixing and deleting these particles we find that det J n 12 where in the last equation we substitute in the MHV solution from section 3. It is also possible to analytically evaluate the Jacobian outside of the MHV sector. Assume 1, 2 ∈ L and gauge-fix and delete particles 1 and 2, and use the formula for the determinant of a block matrix to find that det J n 12 where r, r ∈ R and l, l ∈ L/{1, 2}. As ∂Er ∂σ r is block diagonal, it is comparatively simple to calculate its determinant and inverse. Using equation (B.22) we find the determinant to be det ∂E r ∂σ r = r∈R l<l ∈L ll (ll ) (rl) 2 (rl ) 2 . (B.25) To invert this matrix we see that in the r, r indices it is simply δ rr , leaving the calculation of the inverses of the 2 × 2 blocks. We use the result for a 2 × 2 matrix M that M where the determinant is taken over l, l ∈ L/{1, 2}, combined with the tensor product of spinor and worldsheet indices. In the NMHV sector we see that l = l = 3 and we can calculate the determinant of the remaining 2 × 2 matrix using equation (B.22). As an example at 6 points NMHV we find that det(J 6