Analytical amplitudes from numerical solutions of the scattering equations

The CHY formalism for massless scattering provides a cohesive framework for the computation of scattering amplitudes in a variety of theories. It is especially compelling because it elucidates existing relations among theories which are seemingly unrelated in a standard Lagrangian formulation. However, it entails operations that are highly non-trivial to perform analytically, most notably solving the scattering equations. We present a new Python package ( https://github.com/GDeLaurentis/seampy ) to solve the scattering equations and to compute scattering amplitudes. Both operations are done numerically with high-precision floating-point algebra. Elimination theory is used to obtain solutions to the scattering equations for arbitrary kinematics. These solutions are then applied to a variety of CHY integrands to obtain tree amplitudes for the following theories: Yang-Mills, Einstein gravity, biadjoint scalar, Born-Infeld, non-linear sigma model, Galileon, conformal gravity and $(\text{DF})^2$. Finally, we exploit this high-precision numerical implementation to explore the singularity structure of the amplitudes and to reconstruct analytical expressions which make manifest their pole structure. Some of the expressions for conformal gravity and the $(\text{DF})^2$ gauge theory are new to the best of our knowledge.


Introduction
The scattering equations (SE) first appeared in the litterature in the context of string theory in the '70s [1][2][3] and '80s [4]. They were more recently rediscovered by Cachazo, He and Yuan (CHY) in a series of pioneering papers [5][6][7] demonstrating that the SE provide a set of algebraic equations that are key to an alternative formulation of scattering amplitudes at tree level in d dimensions. Shortly afterwards, this framework was proven to reproduce the correct results for φ 3 and Yang-Mills [8], to generalise to loop level [9,10], and to arise naturally from a worldsheet theory called ambitwistor string [11].
In this alternative QFT formulation the kinematic information of the scattering process is encoded in a set of variables describing the location of punctures on the Riemann sphere. The locations of the punctures are related to the external momenta by the SE. Tree-level amplitudes are obtained by integrating over the position of the punctures on the Riemann sphere, while removing a redundancy coming from Möbius transformations, and imposing the solution of the SE. Alternatively, this integral can be recast as a contour integral around the punctures of the Riemann sphere. The rest of the integrand (called the CHY-integrand) depends on the chosen theory and it has the nice feature of making manifest relations that are hidden in a standard Lagrangian formulation. For instance, the CHY-integrands for Yang-Mills, Einstein gravity and biadjoint scalar theory closely match the KLT relations [12,13].
The main bottleneck for the study of QFTs following this approach is the factorial growth of the number of solutions to the SE. In general, after accounting for Möbius redundancy, the the CHY formulae are supported on (n − 3)! solutions of the SE. More specifically, at three-point there are no free punctures, at four-point the SE have a single rational solution, and at five-point the there two irrational solutions. At six-point there are six irrational solutions which have been shown to be still algebraic in d = 4 [14]. Starting at seven-point in d = 4 and at six-point for general d dimensions the solutions can not be expressed in terms of radicals. At the same time tree-level amplitudes are rational functions of the external kinematics for any phase space multiplicity. Clearly some nontrivial simplification has to occur.
There exist also formulae specific to d = 4 based on the scattering equations refined by MHV degree [15][16][17]. In this case the counting is different and the number of solutions corresponds to the Eulerian numbers.
An intriguing solution found in the literature [18,19] to this factorial growth problem is to obtain the sum of residues from the integral yielding the amplitude without explicitly finding the position of the poles. This powerful approach makes the rationality of the amplitude manifest even when the punctures are irrational. However, as the analytical complexity grows with the multiplicity of the scattering process, even this approach seems to require some form of numerical or semi-numerical reconstruction in order to achieve an analytical expression for the amplitude.
In this paper we develop a purely numerical approach, followed by an analytical reconstruction with the strategy of Ref. [20]. To perform this reconstruction we need an implementation of the CHY formulae which is both sufficiently stable in singular limits and that yields amplitudes with enough numerical precision. We provide code that satisfies these criteria in a Python package which we called seampy (from "Scattering equations and amplitudes with Python").
A publicly available package to compute amplitudes within the CHY framework had already been presented in Ref. [21]. It is based on the scattering equations refined by MHV degree. However, it was not designed to provide amplitudes with the high precision needed by our reconstruction strategy. Furthermore, although the reconstructed analytical expressions we present in section 4 are specific to d = 4, our package provides numerical solutions to the SE in general d dimensions.
This article is organised as follows. In section 2 we review parts of the CHY formalism, in particular the polynomial form of the scattering equations [22], their solutions by means of elimination theory [23,24], and a variety of CHY-integrands [25,26]. The algorithm for solving the SE and the CHY-integrands are implemented in the Python package provided, which is presented in section 3. It provides high-precision floating-point solutions to the SE and numerical amplitudes. In section 4 we make use of this technology to explore the singularity structure of amplitudes and reconstruct explicit expressions in the d = 4 spinor helicity language. Finally, in section 5 we give our conclusions and outlook.

The CHY formalism
In this section we briefly review the theory underlining the CHY formalism and in the next section present its implementation in a Python library. For a more thorough introduction to the subject, with explicit step by step derivations, please consider Ref. [27] and the references therein.
Let us consider the tree-level scattering of n massless particles in d dimensions. We denote with A the set {1, . . . , n}, with k µ a (a ∈ A) the n momenta, and with z a the n special points of the Riemann sphere called punctures. The map from momentum space to the Riemann sphere, as defined in Ref. [5], is the given by where p µ (z) are d polynomials with coefficients depending on the momenta and the punctures. The contour is taken to encircle the punctures. From Eq. (2.1) it can be shown, as a consistency condition, that the following equations have to be satisfied these are the so-called scattering equations. As previously mentioned, the SE are invariant under Möbius transformations SL(2, C), that is under the following mapping Because Eq. (2.3) has effectively three free complex parameters, we can fix the position of three of the n punctures. A common choice in the literature, which we follow throughout this work, is given by Scattering amplitudes for n massless particles A n 1 are then obtained by integrating a CHY-integrand I CHY 2 over the solutions of the SE. This can be achieved either with a normal integral over delta functions, or as a contour integral over the SE. As by prescription where the Möbius measure dω and the modified product symbol are defined as 1 Color ordering is assumed for gauge theories. 2 More details on ICHY are given in section 2.2. For now it suffices to say that in general it is a function of the punctures z, the momenta k and the polarisations . (2.8) By substituting Eq. (2.7) and Eq. (2.8) back into Eq. (2.5) or Eq. (2.6) it can be shown that the amplitude A n is invariant under Möbius transformations. Note that in principle the sets {i, j, k} and {r, s, t} are independent, but in practice they are often taken to be the same for convenience sake. Clearly the requirement of Möbius invariance also imposes a restriction on the valid CHY-integrands I CHY , as we will see shortly.
We would like to use a purely algebraic approach, as it is more amenable to be implementation as computer code. To achieve this we can recast Eq. (2.5) from an integral to a summation by changing variables from the punctures z a to the scattering equations f a . This introduces a Jacobian factor, i.e. the determinant of the Jacobian matrix defined as Again, in the spirit of preserving Möbius invariance, since we have removed punctures i, j, and k from the above δ-function, we also have to remove the corresponding rows form the Jacobian. Similarly, we are not integrating over r, s and t, and therefore those columns have to be removed as well. The matrix of Eq. (2.9) with rows i, j, k and columns r, s, t removed is denoted by φ ijk rst . In the end, the relevant Jacobian for the change of variables, which is independent of the Möbius fixing choice, is given by 3 . We now write Eq. (2.5) for the scattering amplitudes as where j labels the solution of to the SE given by the set of punctures z (j) , which are themselves function of the momenta k. Note that, because of Eq. (2.10) and our choice Eq. (2.11), the Jacobian J introduces the four powers of z 1 = ∞ in the numerator. Therefore, I CHY must come with four powers of z 1 in the denominator for Eq. (2.12) to be sensible. This is a check of Mobius invariance.

Polynomial form of the SE and their solutions
We now turn to the problem of actually finding the solutions to the SE. It is easiest to consider the SE in the form found in Ref. [22], where the SE are reformulated as n − 3 polynomial equations. We can then follow Ref. [23,24] in using an elimination theory algorithm to find the solutions. The SE in polynomial form, which are equivalent to the original SE of Eq. (2.2), are given by where the sets A and S 1 are defined as and where k S and z S are defined as In the above z 1 and z n have already been set to ∞ and 0 respectively, but z 2 is still kept free. This is a system of n−3 polynomial equations (h 1≤m≤n−3 ) in n−2 variables (z 2≤i≤n−1 ). As such it can be solved by using an elimination theory algorithm. The idea underpinning elimination theory is to express the system of equations in matrix form and to introduce more variables and equations until the system is over-specified and yields a consistency condition in the form of det(M n ) = 0. Here we are going to discuss directly the general n case. A more detailed discussion can be found in the original papers of Ref. [23,24] or in Ref. [27].
In general, the aim is to obtain an equation of order (n − 3)! in the ratio z n−1 /z n−2 . The original set of 2 n−4 monomials we wish to eliminate is given by We introduce a auxiliary set which contains (n − 4)! terms. The new set of monomials is then given by which is of length (n − 3)!. Similarly, the new (n − 3)! equations are given by where H T denotes the vector of polynomial scattering equations h 1≤m≤n−3 .
This procedure ensures that the number of monomials matches the number of equations, thus allowing to express the system in matrix form.
Then, by taking partial derivatives of the entries of the extended H of Eq. (2.19) w.r.t. those of the extended V of Eq. (2.18), we could construct the (n − 3)! × (n − 3)! matrix M n whose determinant is the required equation. However, this is not necessary in practice since the matrix M n can be built recursively in a block-matrix format starting directly from the original set h 1≤m≤n−3 and their derivatives w.r.t. z 2≤i≤n−3 . We denote the derivatives with superscripts (M z = ∂ z M ) and we have when written in terms of M i−1 . After the derivative is taken z i−3 is set to zero. M n is then a function of z n−1 and z n−2 only, the required equation of order (n − 3)! in z n−1 /z n−2 is simply det(M n ) = 0, and its roots are the solutions we seek. Note that, as discussed in the introduction, it is feasible to perform this root-finding step analytically only for low phase space multiplicities. Clearly we are not at the end of the calculation yet, because we want values or expressions for the punctures themselves not for ratios. This is achieved by reintroducing one variable at a time in M . More explicitly, we first check with Eq. (2.18) the position in the vector of the variablez we want to reintroduce (say it is the j th entry) then we addz times the j th column of M to its first column, and eventually remove the j th column and the last row. This leads to a matrix of size (n − 3)! − 1 × (n − 3)! − 1 whose determinant will be a linear equation forz. There is one notable exception to this procedure, namely whenz = z 2 we set z 2 = 1 and get a linear equation for z n−2 instead.

CHY-integrands
So far we have treated the theory-independent part of Eq. (2.12). Now we consider the theory-dependent term I CHY . It can be built in a modular way from various building blocks. Here we review the definition of some of those building blocks found in Ref. [25] and in Ref. [26] which we have implemented in the Python package presented in the next section.
Starting from the building blocks that are matrices, we have the 2n×2n anti-symmetric matrix Ψ which is defined block-wise in terms of two n × n anti-symmetric matrices A and B and in terms of a third n × n matrix C. The definitions follow.
Since these are matrices we have to define an operation which converts them to a rank-one object before we can use them to construct I CHY . In the case of anti-symmetric matrices the determinant can be written as a square of a polynomial in the matrix entries. This polynomial is called the Pfaffian and it was shown to be the correct operation to perform. More specifically, since the matrix Ψ has two null vectors and its Pfaffian would be zero, it is necessary to define a reduced Pfaffian PF as where Ψ ij ij again denotes deletion of rows and columns i and j. The same reduction applies also to different arguments, such as the matrix A.
We also consider two scalar building blocks C n and W 1 . C n is a cyclic Parke-Taylor-like factor simply defined as , (2.24) and the W 1 function is defined as 4 I CHY is built from products of pairs of these building blocks. A more detailed analysis reveals that PF (Ψ), C n and W 1 come with a factor of z −2 1 , while PF (A) comes with a factor of z −1 1 . This dictates which combinations are allowed by Möbius invariance (recall that overall we need four powers of z 1 to balance out those in Eq. (2.12)). Table 1 summarises the theories that can be built out of PF (Ψ), C n , PF (A) 2 and W 1 : EG stands for Einstein Gravity, YM for Yang-Mills, BS for Biadjoint Scalar, BI for Born-Infeld, NLSM for Non Linear Sigma Model and CG for Conformal Gravity. The theories labelled with a question mark do not seem to have an agreed upon name, but they are discussed in the reference from which the W 1 function is taken. This is by no means a complete recount of all possible integrands I CHY , but it is sufficient to illustrate the framework. Also note that, as anticipated in the introduction, relations among theories through double copies are now manifest in the structure of the integrands.
Finally, remember that the CHY-integrands are not unique. For instance, a different integrand for conformal gravity is given in Ref. [28]. Table 1: Possible QFTs built out of PF (Ψ) , C n , PF (A) 2 and W 1 . A product is implied between rows and columns, eg: I CHY, EG = PF (Ψ) × PF (Ψ).

Python libraries
In this section we introduce two new packages developed in Phyton 2.7: • seampy (Scattering equations and amplitudes with Python), The former provides high-precision floating-point solutions to the scattering equations in d dimensions and a variety of numerical scattering amplitudes built from their solutions. The latter is used to manipulate and pass a high-precision phase space point as input to the numerical amplitude.
Both packages are available on the Python Package Index. The source code is available on github and the documentation on the associated github pages seampy and lips. Their installation is straightforward thanks to pip: pip install --upgrade seampy # this installs lips as well pip install --upgrade lips # but it can be installed separately The same commands can be used to update the libraries. The --upgrade option ensures that the lastest version is always used. A review of the key features of these packages is now provided. Further examples are given in section 4 and in the appendices A and B.

Solving the scattering equations
In this section we show how to easily obtain solutions for the scattering equations. All the following examples have n = 6.

Computing scattering amplitudes
First of all we can list the theories directly available for computation: To calculate an amplitude we need to generate a phase space point, as in the example for the solutions of the scattering equations: >>> oParticles = Particles(6) # arg. is multiplicity of phase space We then need to declare what quantity we want to compute. This requires us to specify a theory and a multiplicity. For example, biadjoint scalar theory (BS) amplitudes or non-linear sigma model (NLSM) amplitudes can be accessed as follows: Gauge and gravity theories also require an helicity configuration to be specified (the multiplicity is then deduced from it). Note that for gravity theories we are suppressing the repeated helicity sign since we don't have mixed cases such as dilatons. This means that in the following code snippet for conformal gravity (CG) helconf=pmpmpm stands for It is then simply a matter of evaluating any amplitude at the phase space point: >>> oBSAmp(oParticles) mpc(real='#nbr', imag='#nbr') Since most of these helicity amplitudes come with pre-factors of √ 2, we decided to normalise them in such a way that numerical coefficients in analytical expressions are rational fractions and often simply the imaginary unit. This also allows for easier comparison to other codes, which usually adopt such a normalisations. For instance, in the case of Yang-Mills amplitudes the right hand side of Eq. (2.12) is multiplied by 1/( √ 2) n−2 , so that the numerical coefficient in the Parke-Taylor expression for MHV amplitudes is i instead of ( √ 2) n−2 i, where n is the multiplicity of the process. Additional checks that don't require independent implementations of amplitudes include checking the little group scalings, mass dimensions, pole structure (more of this in section 4) or properties such as color ordering. For instance, as a sanity check, we can see that (DF) 2 is color ordered whereas conformal gravity is not. This is shown in the following snippet (we are still using the helconf=pmpmpm amplitudes declared above): >>> oNewParticles = oParticles.image("321456") # swap momenta 1 & 3 >>> abs(oCGAmp(oParticles) -oCGAmp(oNewParticles)) < 10 ** -270 True >>> abs(oDFAmp(oParticles) -oDFAmp(oNewParticles)) < 10 ** -270 False However, picking the correct cyclic permutation of the external legs leaves the (DF) 2 amplitude unchanged as well.
Finally the most stringent tests come from comparing to independent libraries. We have checked all pure gluon (Yang-Mills) tree amplitudes at 3, 4, 5, 6, and 7 point against BlackHat [29] and Yang-Mills, Einstein and conformal gravity against the code of Ref. [21]. They all match, that is their ratio differs at most by a normalisation factor fixed by convention.

Analytical reconstruction
We now consider how to recover analytical expressions for the tree-level scattering amplitudes discussed so far. There are several reasons why analytical expressions are preferable to numerical ones, such as execution speed, numerical stability and general understanding of their analytical structure. The same reconstruction technique can be applied to all the theories from Table 1. In the accompanying files we provide sample analytical amplitudes for all these theories up to six point. The the results are given both in human readable format and as expressions readable by the S@M Mathematica package [30].
In this section, we are going to explicitly discuss only the reconstruction of (DF) 2 and conformal gravity amplitudes, since they are the ones with a less well known analytical structure and therefore the most interesting to analyse. These theories are related by a double copy relation, similar to that between Yang-Mills and Einstein gravity, namely: (DF) 2 × YM ∼ CG. (DF) 2 and conformal gravity present issues with renormalisability and unitarity, since for instance (DF) 2 is built out of dimension-six operators, as implied by the name. Despite this, they are of interest for a few reasons. Namely, one type of conformal gravity arises in Berkovits-Witten twistor string [31], it is the zero-mass limit of a mass-deformed theory that reproduces Einstein gravity in the infinite-mass limit [32], and it may be useful for computing Einstein gravity amplitude in curved backgrounds for cosmological applications [33,34].
More specifically, in the following paragraphs we are going to provide: a) the first complete set of five-point (DF) 2 amplitudes (one of which we could confirm numerically with that found in Ref. [35]); b) an alternative expression to that of Ref. [31] for the fivepoint MHV conformal gravity amplitude; c) results for the leading three-particle sigularities of the six-point amplitudes in the MHV and NMHV helicity sectors. All the amplitudes we present are written in the spinor helicity language and are free from spurious singularities, unless explicitly stated. We think that, in order to obtain similar complete results at six point, it could be necessary to use spurious sigularities, which introduces a further complication in the analysis.
We make use of the high floating-point precision provided by seampy and follow the strategy introduced in Ref. [20]. Briefly summarised, we study the behaviour of amplitudes in singular limits of complex phase space to obtain the poles and their degree. We then study the amplitudes in doubly singular regions to obtain information about the structure of the denominators of the amplitude. Using this information we generate ansätze for the residues of different poles and solve linear systems for the coefficients of bases of spinor expressions in the numerators. If a reconstructed ansatz is correct, once subtracted from the numerical amplitude, it removes a singularity. We repeat the procedure until the amplitude is fully reconstructed.
Explicit examples are discussed in the following subsections.

Five-point amplitudes
(DF) 2 : five-point all-plus (explained example) In contrast to QCD amplitudes, five-point (DF) 2 amplitudes are non zero for all helicity configurations even at tree level. They are color ordered, like QCD, because their CHYintegrand contains the Parke-Taylor-like cyclic factor C n of Eq. (2.24). Therefore, the symmetry group is restricted to cyclic and anti-cylic permutations. It can be generated from two operations, which can be thought of as the rotations and reflections of a pentagon (i.e. the dihedral group D 5 ): (12345 → 23451) and (12345 → −15432) . The minus sign in the reflection comes from the partity operation applied to vector particles (J P = 1 − ). In total the group contains 10 elements (including the identity).
Following this same procedure with the rest of the spinor invariants we obtain a first look at the analytical structure of the all plus amplitude: where N is some numerator structure. Two comments are now in order. Firstly, note that the adjacent particle singularities are of second order. This reflects the fact that this theory has a quartic propagator instead of the usual quadratic one. Secondly, although in this case it is possible to obtain an expression for the numerator N , it is often not feasible to do so in this single fraction representation, especially with higher point amplitudes; and even when it is possible, the result is complicated and obscures the structure of the amplitude.
In order to obtain a compact representation, we want to write the amplitude as a sum of fractions, each of which should have a simpler denominator structure than the expression above. It is generally convenient to start by considering the double poles, since they make it difficult to numerically access the corresponding simple poles. We study doubly singular limits, that is regions of phase space where pairs of spinor invariants vanish. In practice, this can be done with the same code snippet as above, by replacing the oParticles.set function with the oParticles.set_pair one. For example, for the pair 12 , 23 we have: >>> oParticles.set_pair(" 1|2 ", 10 ** -30, " 2|3 ", 10 ** -30) By repeating the same procedure with all pairs involving 12 and recording the behaviour of the amplitude in the corresponding doubly singular limit we can generate Table 2. Since 12 is already a double pole, it is not likely for any other invariant appearing with a 2 in the table to be in the same denominator as 12 2 . Therefore, we make an ansatz where only 34 and 45 (as simple poles) appear together with 12 2 . More rigorously, we conjecture that: To check whether the above is true or not, we start by noting that the amplitude has mass dimension 5 of 1 and little group weights 6 of [−2, −2, −2, −2, −2]. Therefore, the numerator in the RHS must have mass dimension 5 and little group weights [0, 0, −1, 0, −1] in order to match the LHS. We then generate a complete set of linearly independent products of spinor invariants consistent with these constraints. In this specific case the basis contains 20 independent entries: Note that, the basis would have 290 entries if we were to generate it for the numerator of Eq. 4.2. Moreover, since we are not working in a generic phase space region but in the limit of small 12 , it turns out that 10 of the 20 basis elements only contribute to the O( 12 −1 ) part of Eq. 4.3 and thus can be ignored. We can now generate 10 random phase space points in the 12 → 1 region and solve for the coefficients of the 10 elements. The solution has only one non zero coefficient: To obtain the remaining four double poles, we can simply symmetrise the expression for the 12 double pole by applying the following cyclic permutations: Once an expression for a particular pole has been reconstructed, it can be numerically subtracted from the amplitude and the left over quantity will not contain that particular singularity anymore. Its singular limits can then be studied, ansätze made and reconstructions performed until all the poles have been successfully obtained and the amplitude fully reconstruced.
The final result for the all plus (DF) 2 amplitude follows. On the left hand side we give the amplitude written using the symmetries discussed above. This is the format used throughout the rest of the article. For the sake of clarity, below we reproduce on the right hand side the same expression with the meaning of the symmetries made explicit. 5 Natural units are assumed. In practice the mass dimension can be numerically obtained by re-scaling all the momenta. 6 Little group transformations modify spinors, while leaving four momenta unchanged. Thus little group scalings can be numerically obtained by re-scaling the spinors. For more details see Ref. [36].
A (DF) 2 (1 + , 2 + , 3 + , 4 + , 5 + ) = This MHV amplitude is the only one we could already find in the litterature, specifically in Ref. [35], where it was written in terms of Mandelstam invariants. The expression we provide is more concise, makes its symmetry explicit and is free from spurious singularities.
We have numerically checked that the two expressions agree. The one we found follows. (DF) 2 : five-point MHV (non-adjacent) The following is the last independent five-point amplitude. All others can be obtained by permutations and/or conjugation of the amplitudes presented here.

Conformal gravity: five-point MHV
An all-multiplicities expression for MHV conformal gravity amplitudes exists thanks to work by Berkovits and Witten [31]. Here we present an expression specific to five point which makes manifest the absence of terms with pairs of double poles. In order to convey the increse in complexity that a six-point amplitude entails here we present an expression for the three-particle double poles as well as for the simple poles of non-adjacent three-particle singularities in a six-point MHV (DF) 2 amplitude. In the above expression N would contain several thousand terms. It is therefore crucial to identify appropriate ways to perform a partial fraction decomposition, since smaller denominators would in turn imply smaller numerators and thus easier systems of linear equations to generate and solve. However, further studies will be necessary to check whether such a decomposition requires the introduction of spurious singularities, like for NMHV amplitudes in Yang-Mills, and if so what form these spurious poles would take.

Conformal gravity: NMHV (partial)
To conclude, we present an expression for the three-particle double poles in the six-point NMHV conformal gravity amplitude. To the best of our knowledge this is the first analytical result, albeit a partial one, for NMHV conformal gravity amplitudes. Here N would contain even more terms then in the six-point (DF) 2 example. Similar expressions where the symmetries of the poles are made manifest are also possible in Einstein gravity amplitudes, for example the following represent the three-particle simple poles in the six-point NMHV sector. However, this symmetric approach, which is also free from spurious singularities, makes it highly non trivial to obtain the rest of the amplitude (i.e. the numerator N ). Indeed, the compact expressions that we are aware of come from BCFW recursions and have a quite different structure: We have reproduced this result already known in the literature by applying our analytical reconstruction strategy to a single BCFW factorisation channel at a time, which is significantly simpler than the full amplitude 7 . Compared to the previous partial result, we note that this representation manifestly does not contain two-particle Mandelstam invariants, but introduces many spurious singularities and hides the symmetries which were manifest in the above partial result. The strategy of studying a factorisation channel at a time could prove fruitful also in the case of conformal gravity and (DF) 2 amplitudes, but the quartic propagator introduces a significant complication in the BCFW recursion.
In fact, the usual A L A R /p 2 factorisation is broken by the presence of higher order poles in the Laurent expansion in the shift parameter. We attempted to achieve such a factorisation by means of a Taylor expansion of the numerator A L A R around the pole. However, this involves taking a derivative with respect to the shift parameter which in turns requires the amplitudes to be well defined in the neighbourhood of the factorisation point. This seems to be equivalent to the factorisation formula (Eq. 2.18) given in Ref. [32], where the derivative is implicit in the fact that we have to take the zero mass limit of expressions like (A L (m 2 )−A L (0))/m 2 . This would also explain why our approach fails: the amplitudes we use are well defined only exactly at the factorisation point, where the legs are on-shell and massless.
However, we do have the six-point amplitude through the CHY formula and there is no need to generate it recursively from lower point amplitudes. At the same time, we expect single factorisation channels to have an easier analytical structure than the full amplitude. This suggests to still look at the amplitude via the residue theorem: 7 In this case a 21] shift was used.
We can then study one term in the sum in the RHS at a time. Note that the simultaneous need to generate singular phase space limits and to numerically extract the residue from a Laurent expansion in some cases requires to increase the working numerical of precision.
The numerator N , having mass dimension of 46, is unfortunately still too complicated to be determined. We see that the conformal gravity residue has more poles and poles of higher order compared to Einstein gravity one, as well as some spurious singularities of order higher than one. Furthermore, note that for this shift the contour integral vanishes for Einstein gravity but not for conformal gravity. Therefore, in the latter case we would have to include a boundary term coming from the residue at infinity. Some of the other possible shifts have the advantage of vanishing on the contour, but the structure of the residues remains similarly complicated. Further work will be required to see whether a reasonably compact analytical expression can be obtained for these residues.

Conclusion and outlook
In this article we have briefly reviewed the CHY formalism for massless tree-level scattering, and more specifically the problem of solving the scattering equations and applying their solutions to CHY-integrands.
In order to overcome the analytical complexity of the computation, we have developed a Python package (seampy) which allows to numerically solve the scattering equations and to computate tree amplitudes with high floating-point precision for the following theories: Yang-Mills, Einstein gravity, biadjoint scalar, Born-Infeld, non-linear sigma model, Galileon, conformal gravity and (DF) 2 .
Finally, we have discussed how to recover analytical expression in the spinor helicity language from numerical evaluations. In particular, we have presented the first complete set of five-point (DF) 2 amplitudes, a new form for the five-point MHV conformal gravity amplitude and a discussion with partial results for six-point amplitudes in both (DF) 2 and conformal gravity.
In the accompanying files we have provided sample analytical amplitudes for all mentioned theories up to six point. The the results are given both in human readable format and as expressions readable by the S@M Mathematica package.
Let us remark the fact that despite not all the solutions to the scattering equations are rational (except at three and four point), and in some cases not even be expressible in terms of radicals (beyond six point), the tree-level amplitudes built from them are purely rational functions. This is made clear by reconstructing explicit rational analytical expressions from numerical evaluations. The expressions we obtain are usually compact, with a clear symmetry structure when available and free from spurious singularities, unless explicitly stated.
We have observed that complexity increases significantly from five-point to six-point amplitudes and given explicit examples. In the previous section we discussed ways to look at simpler building blocks rather than the full amplitude all at once, such as a modified BCFW recursion for the quartic propagators of conformal gravity and (DF) 2 or a more naive application of the residue theorem. These approaches seem promising, since they can still be carried out numerically while resulting in simpler structures to which apply the analytical reconstruction. However, more remains to be done to make this feasible in practice for the more complicated theories.
Finally, going forward it might be interesting to use this numerical approach to the CHY formalism together with the analytical reconstruction tools to look at other interesting quantities such as double copy structures, BCJ numerators, amplitudes with mixed particle content, and loop-level amplitudes.

Appendices A lips (phase space generator)
In this first appendix we present in more details the lips Python package. It is an objectoriented high-precision floating-point phase space generator. It is not the focus of this work, but it is necessary in order to pass a sufficiently precise phase space point to the scattering equations solver function or numerical amplitude object.
The lips phase space generator is built on two layers. The lower one, called Particle, describes the kinematics of a single particle. Though setters and getters, it provides selfupdating numerical tensors for the left and right spinors, four vectors and rank two spinors. This means that if, say, the value of the four momentum is changed, then the values of the spinor attributes are immediately recalculated to reflect the change. We can see the naming conventions in the following code snippet: >>> oParticle = Particle() >>> oParticle.l_sp_u # left spinor with index up (λα) >>> oParticle.r_sp_d # right spinor with index down (λ α ) >>> oParticle.four_mom # four momentum with index up (P µ ) >>> oParticle.r2_sp # rank two spinor (Pα α ) By default the Particle object is initialised with random complex momenta. However, this can be overruled by specifying the optional paramenter real_momentum=True. A custom value for any of these attributes can also be passed. For instance, we can set the momentum to be along the x axis: >>> oParticle.four_mom = [1, 1, 0, 0] The second layer is a list subclass, called Particles. It is a base-one list of Particle objects with several methods attached associated to it. The reason why the list is rebased to start from 1 instead of 0 is simply to match the notation in the amplitudes community. As we have observed, it is initialised as follows: >>> oParticles = Particles(6) # argument is the multiplicity It also accepts an optional parameter, now called real_momenta, which is by default set to False, and which gets automatically passed down to all the Particle objects in the Particles list, thus generating a complex or real phase space point.
Furthermore, as discussed in conjunction with the analytical reconstruction, the Particles phase space can be manipulated to generate specific configurations. For instance, we can generate phase space point with vanishing angle bracket 12 by calling: