An Algebraic Approach to the Scattering Equations

We employ the so-called companion matrix method from computational algebraic geometry, tailored for zero-dimensional ideals, to study the scattering equations. The method renders the CHY-integrand of scattering amplitudes computable using simple linear algebra and is amenable to an algorithmic approach. Certain identities in the amplitudes as well as rationality of the final integrand become immediate in this formalism.


Introduction
In last couple of years, amazing progress has been made by Cachazo, He and Yuan [CHY] in a series of papers [1][2][3][4][5], where tree-level amplitudes of a host of quantum field theories can be calculated using solutions of a set of algebraic equations. These are called the scattering equations and appear in the literature in a variety of contexts [6][7][8][9][10][11][12][13][14].
The mysterious relationship between the CHY approach and the standard QFT paradigm has been explained from different points of view. In [15], using the BCFW on-shell recursion relation [16,17] the validity of the CHY construction for φ 3 theory and Yang-Mills theories has been proven. A broader understanding is achieved using ambitwistor string theory [18][19][20][21][22][23][24][25][26][27][28], where using different world-sheet fields, different integrands in the CHY approach for different theories -which we will call CHY-integrands, a function of the coordinates z i in a Riemann surface -have been derived alongside with the natural appearance of scattering equations. A nice point of ambitwistor approach is that it provides the natural framework for loop scattering equations as studied in [21,24], which lead to a breakthrough in [28]. A third understanding is given in [29], where inspired by the field theory limit of string theory, a dual model has been introduced, based on which a direct connection between the CHY approach and the standard Feynman diagram method has been established not only at the tree-level in [30,31], but also at the one-loop level (at least for φ 3 theory) in [32] (see also [33]) 1 .
Although conceptually the CHY approach is remarkable and very useful for many theoretical studies of properties of scattering amplitudes, when applying to real evaluation, one faces the problem of solving scattering equations, which has (n − 3)! solutions in general. Furthermore, when n ≥ 6, one encounters polynomials of degree exceeding five, rendering analytic solutions in radicals hopeless. Nevertheless, while the solutions can be very complicated, when putting them back into the CHY integrand and summing up, one obtains simple rational functions. These observations have led people to wonder if there is a better way to evaluate the CHY-integrand without explicitly solving the scattering equations. In [35], using classical formulas of Vieta, which relate the sums of roots of polynomials to the coefficients of these polynomials, analytic expression can be obtained without solving roots explicitly. More general algorithms are given by two works. In one approach [36], using known results for scalar φ 3 theory, one can iteratively decompose the 4-regular graph determined by the corresponding CHY-integrand to building blocks related to φ 3 theory, thus finishing the evaluation. In another approach [30,31], by careful analysis of pole structures, the authors wrote down a mapping rule, so that from the related CHY-integrand, one can read out contributions of corresponding Feynman diagrams.
Both approaches are powerful and have avoided the need of solving the scattering equations explicitly. Furthermore, based on these perspectives, especially the mapping rule, one can use Feynman diagrams to construct the CHY-integrand. These results produce a very interesting phenomenon: two different CHY-integrands can produce the same result. For example, there are two very different CHY-integrands for scalar φ 4 theory: one is given in [5], while another one is given in [30,31]. We are naturally led to wonder how to explain the equivalence of different CHY-integrands.
In fact, as a rational function of coordinates z i on a Riemann surface, the equivalence can occur on three different levels.
of course rather trivial and in order to proceed to the other two levels of equivalences, we need to change our viewpoint to algebraic geometry, i.e., to transform the scattering equations to a set of polynomials of (n − 3) variables, defining an ideals I; 2. The difference of two CHY-integrands can be written as J(I 1 − I 2 ) = P Q where both P, Q are polynomials and J is the Jacobian we will review shortly. If P belongs to the ideal I, then for each solution of the scattering equations P = 0, thus at the second level we say that I 1 is equivalent to I 2 ; 3. However, in practice, most of the time something more complicated happens and we find that though P does not belong to the ideal I and J(I 1 − I 2 ) = 0 when and only when we sum over all solutions. If this happens, we say that I 1 is equivalent to I 2 at the third level. It is clear that this is the most involved situation, and indeed, in practice this is the most frequently encountered.
Motivated by the above considerations and bearing in mind that indeed the most conducive perspective on studying the scattering equations is through the language of algebraic varieties and polynomial ideals [37,38], we turn to this method of attack. The above problem thus translates to finding the sum over the rational function P Q evaluated at the roots of a zero-dimensional ideal I, and testing whether the sum is zero. Luckily, there is a theorem in commutative algebra, due to Stickelberger, which addresses the situation [39]. We will discuss the theorem and the associated algorithm in illustrative detail. It turns out that this method not only checks the equivalence at the third level, but also evaluates the integration without solving the scattering equations. In this sense, it is in the spirit of the methods in [36] and [30,31]. Although it is sometimes less efficient compared to these two methods, it does provide a very different angle to approach the problem and could have very advantageous repercussions.
The structure of the paper is as follows. We begin with a brief review of the tree-level scattering equations in §2, before laying down the foundations of the theory of zero-dimensional ideals in §3, especially that of companion matrices. We then illustrate the technique with ample computational examples in §4, before concluding with remarks in §5.

Review of Tree-level Scattering Equations
In this section, we offer a brief review of tree-level scattering equations and the reader is referred to [1][2][3]5] for details. The scattering equations are given by where s ab = (k a + k b ) 2 = 2k a · k b , and k a with a = 1, 2, ..., n are n massless momenta for n-external particles and z i are complex variables living on CP 1 with n punctures. Although there are n equations, only (n − 3) of them are linear independent after using the momentum conservation and massless conditions which translate to the following three relations which are, in fact, the consequence of the SL(2, C) symmetry on the CP 1 . Because of this, we can insert only (n − 3) delta-function. To make sure the result does not depend on which three equations have been removed, we make following combination and define the measure 2 with z ij = z i − z j . With the above, the general tree-level amplitude is given by where dω = dzrdzsdzt zrszstztr comes after we use the Möbius SL(2, C) symmetry to fix the location of three of the variables z r , z s , z t by the Faddeev-Popov method. Different QFTs give different forms of the CHYintegrand F(z). Invariance under the Möbius transformation requires F(z) to have proper transformation behaviors, i.e., under z = az+b cz+d , we have To simplify expression (2.4) further, we integrate out the delta-functions to arrive at the key expression where three arbitrary indices i, j, k correspond to three removed scattering equations while three arbitrary indices r, s, t correspond to the above mentioned three fixed locations. The sum is over the solution set of the scattering equations, which is generically a discrete set of points. Furthermore, in the above, the Jacobian matrix Φ is calculated as (a for rows and b for column) and |Φ| rst ijk is the determinant of Φ after removing the i-th, j-th and k-th rows and r-th, s-th and t-th columns.

Specific Examples:
Now we list some examples in the literatures [2,3] (more can be found in [5]). According to the CHY formula, the integrand unifying scalars(b = 0), gluons(b = 1) and gravitons (b = 2) is given by where the sum is over permutations on n elements by the symmetric group S n , up to cyclic ordering of Z n , Ψ is a 2n × 2n antisymmetric matrix defined by Ψ = (where t is the transpose of the matrix), with A, B, C being n × n matrices with components and Pf Ψ is the reduced Pfaffian (square-root of the determinant) of Ψ defined by where 1 ≤ i, j ≤ n and Ψ ij ij is the matrix Ψ removing rows i, j and columns i, j. We recall that the Pfaffian of a 2n × 2n antisymmetric matrix can be computed as where sgn(σ) is the signature of σ ∈ S 2n . Importantly, Pf Ψ ij ij is non-zero on the solutions of scattering equations, while Pf Ψ is independent of the choice of i, j.
Specifically, we have that • For color-ordered bi-adjoint scalar φ 3 theory, (2.12) • For color-ordered Yang-Mills theory with ordering {1, 2, ..., n}, Pf Ψ . (2.13) • For gravity, (2.14) Having presented the above examples, let us go back to (2.6). As is clear from the expression, the right hand side is a rational function in the complex variables z i . To employ methods developed in algebraic geometry, we need to associate solutions to a zero-dimensional algebraic variety defined by some polynomials. In other words, we should rewrite E a defined in (2.1) to an equivalent polynomial system. This has been done in a beautiful paper [37], where it has been shown that scattering equations are equivalent to following set of polynomials where the sum is over all n! (n−m)!m! subsets S of A = {1, 2, ..., n} with exactly m elements and k S = b∈S k b and z S = b∈S z b . The algebraic geometry, notably the affine Calabi-Yau properties of (2.15), has been investigated in [38].
A very useful observation made in [15,37] is that If all k 2 S = 0, then values of z a are all distinct. The set (2.15) has not fixed gauge. One of the choice of gauge will be to set, as is standard with points on CP 1 , the three points z 1 = ∞, z 2 = 1 and z n = 0. Under this choice, the set of polynomial is reduced to In summary,h defines a zero-dimensional ideal in the polynomial ring in n − 3 variables. Then, using the standard Bézout's theorem, the number of points in this ideal (solutions of the scattering equation) Instead of computing the amplitude with formula (2.6) by summing over all solutions of scattering equations, we will show in next section that, using the companion-matrix method, we can compute the amplitude A n = sol P Q as the trace of certain matrix composed of so-called companion matrices T z i without the explicit solutions of scattering equations.

The Mathematical Framework
As mentioned in the introduction, it is expedient to consider the problem within the framework of ideal theory. Our problem is thus the following.

Problem:
Let I = f i be a zero-dimensional ideal in R = C[x 1 , . . . , x n ] generated by f i=1,2,...,k (x 1 , . . . , x n ) ∈ R and let r(x 1 , . . . , x n ) be an arbitrary rational function in the fraction field of R. Because dim C I = 0, I = N j=1 {z j } is a discrete set of, say N , points. We wish to evaluate where each summand is an evaluation of p at one of the discrete set of zeros z j . In particular we wish to test whether this sum is 0. This is the level 3 equivalence mentioned in the introduction.
Of course, the idea is to solve this without explicitly finding the roots z j . This can be done using the technique of companion matrices [40] (cf. also [41]). Suppose a Gröbner basis for I has been found for some appropriate monomial ordering and B is an associated monomial basis for I, which can be seen as a vector space of dimension d. Then the multiplication map by the coordinate variable x i is an endomorphism of quotient rings. In the basis B of monomials, this is a d × d matrix and is called a companion matrix. Clearly, {T i } all mutually commute and thus can be simultaneously diagonalized.
We have the following [39]: The complex roots z i of I are the vectors of simultaneous eigenvalues of the companion matrices T i=1,...,n , i.e., the corresponding zero dimensional variety consists of the points: We point out that the original statement of the theorem is in terms of annihilators in algebraic number theory and is perhaps a little abstruse. Fortunately, the computational algebraic-geometry community has rephrased this into the readily usable form of companion matrices [40,42]. In particular, we have the following important consequence: where the evaluation of the rational function r on the matrices T i is without ambiguity since they mutually commute.
We remark that because r is rational, whenever the companion matrices appear in the denominator, they are to be understood as the inverse matrix.

Warmup
Before proceeding to examples in our context, we present two simple exercises to demonstrate our algorithm. Computations can be made in Macaulay2 [42] or Singular [43], or the latter's interface with Mathematica [44]. Let We know, of course, that there are 5 roots Now we consider two functions, where one is polynomial and another, rational: p(x, y, z) = 3x 3 y + xyz, Q(x, y, z) = 3x 3 y + xyz It is easy to find, after summing over the solutions, that We now show how the companion matrices work without finding the roots (3.3) explicitly.
In the lex ordering of x ≺ y ≺ z, the Gröbner basis and the monomial basis are, respectively, Therefore, the sum over the roots of p is and we have nice agreement with (3.4). Before going to examples of scattering equations, let us give some remarks. First, the theorem in its original form is for polynomial test functions r, while functions we will meet in scattering equations are rational functions, i.e., the form P Q with both P, Q are polynomials. Luckily, the theorem and corollary can be generalized trivially since we can diagonalize companion matrices simultaneously because the next remark. Now, there is a second part of the theorem which states that the companion matrices can be simultaneously diagonalized if and only if the ideal I is a radical ideal. That is, there are no multiple roots. However, as shown in [15], if all k 2 S = 0, the solutions of z i will all be different, so we indeed have a radical ideal and find simultaneous eigenvalues readily.
Third, since there are (n − 3)! solutions, the size of T i will be in general d = (n − 3)! which will become very large with n. Although with this counting, the efficiency of the method may be arguable, it does make the following property manifest: after summing over all solutions, the final result must be rational functions of k, .
-8 -In the following, we will use several examples to demonstrate the companion matrix method. The n = 4 case is simple. The companion matrix is 1-dimensional, equaling to the single solution of scattering equations. We compute the amplitudes in scalar φ 3 , Yang-Mills and gravity theories to show the validity of the method. For n = 5, we first study the amplitude of scalar φ 3 theory, and show that the amplitude-level identity can be understood by the fact that the trace of matrix is a linear mapping, and use it the explain a 7-point identity proposed in [31]. For the amplitudes of Yang-Mills and gravity theories, we will show that the companion matrix method indeed produce the correct amplitudes.
For n = 6, the scalar φ 3 theory will be shown to detect the pole structures so that the amplitude can be constructed by setting appropriate kinematics. The Next-MHV gluon amplitude is also presented as an example to show the validation of companion matrix method in a more difficult situation. Finally, for n = 7 amplitudes of scalar φ 3 theory, we demonstrate that, when companion matrices are computed in the diagonal form, the diagonal elements of the integrand matrix (which we recall to be an (n − 3)! × (n − 3)! matrix for n-points) have one-to-one mapping to the integrand computed at the (n − 3)! solutions of the scattering equations, so they are not only equivalent at the amplitude level, but also at the level of each solution as indicated by Stickelberger's theorem.

Four-Point Amplitudes
The n = 4 case is trivial. There is only 4 − 3 = 1 variable left, so the companion matrix is just a complex number. Let us remove three scattering equations E 1 , E 2 , E 4 and gauge-fix three points z 1 = ∞, z 2 = 1 and z 4 = 0. The remaining one scattering equation is We can define the ideal It is a linear function, so the Gröbner basis and monomial basis are trivially The polynomial reduction of z 3 B = {z 3 } over Gröbner basis of ideal I gives the remainder { s 34 s 23 +s 34 }. Thus in the quotient ring, the companion matrix is given by We now proceed to the three cases of concern.

Scalar φ 3 Theory
For the 4-point amplitude in scalar φ 3 theory, we wish to compute (recall that the three points z 1 , z 2 , z 4 have been gauge fixed) where we have used the simplification so that the factor 1/z 2 23 z 2 34 cancels the numerator of 1/|Φ| 124 124 . We see that the final expression is summed over the (discrete) solution set of the scattering equation which is rather trivial here. The summand is a rational function in the free variable z 3 which we define as P/Q; of course, P = 1 here and Q will be used later.
Finally, using the simple expression for the companion matrix T z 3 from (4.3), we have after some identities between Mandelstam variables have been used. This is indeed the same answer as the standard known result as given in the introduction.

Yang-Mills Theory
For 4-point amplitude in Yang-Mills theory, we want to compute (under gauge-fixing To avoid the divergence when taking the limit z 1 → ∞, one of the removed rows(columns) in Ψ should be 1, otherwise some terms in Pf Ψ 8×8 would lead to infinity. Let us then choose the reduced Pfaffian as The large z 1 dependence of Pf Ψ 8×8 is then 1/z 2 1 , and together with the factor from the scalar part, we obtain a finite integrand when taking the z 1 → ∞ limit. Explicitly, the new matrix Ψ ≡ Ψ 12 12 is a 6 × 6 matrix, whose Pfaffian is given by The reduced Pfaffian Pf Ψ 8×8 in this case is a rational function with denominator z 2 3 (z 3 −1). Together with the factor 1/z 23 z 34 = 1/z 3 (z 3 − 1), they cancel the numerator of 1/|Φ| 124 124 , leaving a z 3 in the denominator of integrand.
Therefore, it is immediate that the numerator of the integrand comes entirely from the numerator of the reduced Pfaffian: The denominator of the integrand, on the other hand, is where Q is the denominator of integrand for the scalar φ 3 theory from (4.4). In summary, by computing Tr(P Q −1 | z 3 →Tz 3 ), we arrive at The pole structures are similar to the scalar φ 3 theory, while the terms without poles come from the gluon four-vertex. Of course, by momentum conservation and the property i k i = 0, we can further write the above result as a function of all independent kinematics, for example by using identities j k 4 = − j k 3 − j k 2 − j k 1 and 4 k 4 = 0. This result agrees with the one computed directly by Feynman diagrams.

Gravity
For the 4-point amplitude in gravity, we want to compute under the gauge-fixing z 1 = ∞, z 2 = 1, z 4 = 0. As in Yang-Mills theory, we choose the reduced Pfaffian as and as above, we know that the squared reduced Pfaffian (Pf Ψ 8×8 ) 2 is a rational function with denominator z 4 3 (z 3 − 1) 2 . This cancels the numerator of 1/|Φ| 124 124 , leaving a z 2 3 in the denominator of integrand, so that the numerator of integrand equals to the square of numerator of reduced Pfaffian: while the denominator of integrand is Combining all together, we have where Q is the denominator of integrand for scalar φ 3 theory from (4.4) and the expressions for P YM and Q YM are given in (4.11) and (4.12). In the present case of n = 4, there is only one solution for scattering equations, and the companion matrix is really 1-dimensional in (4.3), so although in general By BCJ relation [45], we can rewrite this to the familiar one 20) in agreement with the known result by KLT relation [46][47][48][49][50].

Five-Point Amplitudes
For n = 5 amplitudes, there are five scattering equations, but only two of them are independent. Under the gauge-fixing z 1 = ∞, z 2 = 1, z 5 = 0, the Dolan-Goddard's formula [15] gives: We can solve these two equations to get two solutions: We note that the companion matrices actually formally "live" in the ideal I itself by satisfying scattering equations, i.,e., and likewise, s 45 T z 3 + s 35 T z 4 + s 25 T z 3 T z 4 = 0 2×2 . This is, of course, a general property by construction since the companion matrices are constructed as multiplication (on a particular basis), so that substituting into the defining polynomials would vanish in the quotient ring. The situation is very much analogous to the classical result of Cayley-Hamilton that a matrix satisfies its own characteristic polynomial. It is worth to emphasize this discussion as

COROLLARY 2
The companion matrices satisfy the defining polynomials of the given ideal.
The above corollary shows some kind of equivalence between solutions of scattering equations and companion matrices of monomial basis over the Gröbner basis of scattering equations. With these companion matrices, we now proceed to compute the trace of the integrands to obtain the amplitude for different theories.

Scalar φ 3 Theory
The 5-point amplitude of scalar φ 3 theory is given by where we have used that and as above, defined the appropriate P and Q, which are, explicitly, Now, we wish to compute the trace of the matrix P Q −1 upon substituting z 3 and z 4 by their associated companion matrices, instead of summing over all the complicated solutions of the scattering equations. In other words, we should replace the variables z 3 , z 4 as T z 3 , T z 4 in the integrand, i.e., P = P | z 3 →Tz 3 ,z 4 →Tz 4 , Q = Q| z 3 →Tz 3 ,z 4 →Tz 4 (Hereafter we will always use P , Q to denote the matrices after replacing z i to T z i ). The product of variables z 3 , z 4 changes to the product of matrices T z 3 , T z 4 , and since the companion matrices are commutable, their order does not matter in here. Then we should compute the inverse of matrix Q , and the final result is given by Tr(P Q −1 ).
Recalling that the physical poles appearing in the color-ordered amplitude are s 12 , s 23 , s 34 , s 45 , s 15 , we can define them as the independent Mandelstam variables, and rewrite all the other Mandelstam variables in P, Q, T z 3 , T z 4 by using following identities:  which agrees with the known result [30,31]. Let us further consider an example, corresponding to the two-cycles 3 {(1, 2, 3, 4, 5), (1, 3, 5, 2, 4)}, in the language of [30,31,36]. Using the CHY-integrand defined by above two-cycles, we have which is represented by the so-called pentacle diagram (shown in Figure 1)from the view of integration rules. Using the mapping rule given in [30,31], the answer is known to be zero. By directly computing the trace, we indeed find that Tr(P 1 Q −1 1 ) = 0 and confirms this result. In fact, for this example, although CHY-integrands of A 5 and A 5 are different, after simplification, their difference appears only in the numerator, i.e., where they share the same denominator Q, but different numerators After putting back the companion matrices, we find that Realizing that the polynomials have the simple relation we obtain the identity amongst these amplitudes as Tr(P 2 Q −1 ) = Tr((P 3 + P 4 − P 1 )Q −1 ) = Tr(P 3 Q −1 ) + Tr(P 4 Q −1 ) + Tr(P 1 Q −1 ) Above example demonstrates an idea how to find relations among different amplitudes. Starting from different CHY-integrands, we can equalize their denominators by multiplying proper polynomial both at the denominator and the numerator. After that, the relations among different amplitudes can be understood from the relations among different numerators. Let us demonstrate above idea by another example, i.e., the 7-point amplitude-level identity given by eq.(3.7) of [31], viz., amplitude obtained from the CHY-integrand Under gauge-fixing z 1 = ∞, z 2 = 1, z 7 = 0 and excluding the 1-st, 2-nd and 7-th scattering equations, the Jacobian is Thus we can immediately get the numerator of integrand after inserting the above three terms. The first term gives while the other two terms give Note that while the trace Tr((P 1 − P 2 − P 3 )Q −1 ) is zero. Note also the following decomposition so that we can write P 1 − P 2 − P 3 = P 4 + P 5 , with respectively with Tr(P 4 Q −1 ) = 0, Tr(P 5 Q −1 ) = 0. We thus conclude that strictly speaking, the amplitude-level identity between (4.35) and (4.36) is up to some CHY-integrands which have vanishing amplitude. More explicitly, the identity (4.35)=(4.36)+(4.45) holds exactly at the integrand-level, while (4.45) has vanishing final result, so that (4.35)=(4.36) holds at the amplitude-level. This provides the amplitude-level identity an explanation from the basic linearity of the trace.
We now follow the standard computation procedure: 1. Write down the expressions for |Φ| 125 125 and Pf Ψ 10×10 , and work out P YM (z 3 , z 4 ), Q YM (z 3 , z 4 ); 2. Replace the variables z i 's by companion matrices T z i , as P = P | z i →Tz i , Q = Q| z i →Tz i ; 3. Compute the inverse of Q and the trace Tr(P Q −1 ).
The result for un-specified helicities is quite lengthy. For illustration, let us consider the 5-point amplitude with helicity A YM 5 (g − 1 , g − 2 , g + 3 , g + 4 , g + 5 ). The polarization vector is defined as and we choose the reference momenta as r 1 = r 2 = k 3 , r 3 = r 4 = r 5 = k 2 . Thus settled, the only surviving products of polarization vectors are − (k 1 )· + (k 4 ) and − (k 1 )· + (k 5 ). After imposing momentum conservation for ± (k i ) · k j to reduce the ambiguity, we can simplify the 8 × 8 matrix Ψ ≡ Ψ 12 12 as This greatly simplifies the result of reduced Pfaffian, which reads, after our gauge-fixing, where we recall again that i,j = i j , κ i,j = i k j . The factor of scalar part 1/(z 12 z 23 z 34 z 45 z 51 ) after gauge fixing is , and the Jacobian |Φ| 125 125 is the same as in the scalar theory, where Q YM is a polynomial of z 3 , z 4 and Mandelstam variables, and it is also the denominator of integrand. The numerator of 1/|Φ| 125 125 cancels the denominator of Pf Ψ and that of scalar part, leaving a factor z 3 (z 3 − 1)(z 4 − 1) in the numerator. Combined with the numerator N Ψ of Pf Ψ 10×10 , they contribute to Then it is straightforward to apply the replacements P YM (T z 3 , T z 4 ) = P YM (z 3 , z 4 )| z i →Tz i , Q YM (T z 3 , T z 4 ) = Q YM (z 3 , z 4 )| z i →Tz i , and compute the trace Tr(P Y M (Q YM ) −1 ). To make the computation more efficient, we can firstly apply the polynomial reduction of P YM (z 3 , z 4 ), Q YM (z 3 , z 4 ) over GB(I). The remainders R(P YM ), R(Q YM ) are polynomials of z 4 only, since the monomial basis is {1, z 4 }. Then we can proceed by replacing z 4 → T z 4 for the remainders, and compute the corresponding trace. This gives the same result as with the original P YM , Q YM , but the computation would be much faster. With Mathematica, we obtain  Let us consider the gravity amplitude A G 5 (1 −− , 2 −− , 3 ++ , 4 ++ , 5 ++ ), so that we can use the same reduced Pfaffian Pf Ψ 10×10 as in the Yang-Mills case. Here, we do not have the factor of scalar part, but the square of the factor of the reduced Pfaffian. The numerator of 1/|Φ| 125 125 cancels the squared denominator of reduced Pfaffian z 2 3 z 2 4 (z 4 − 1) 2 (z 3 − z 4 ) 2 , leaving a factor of (z 3 − 1) 2 in the numerator. Hence, we have Q G = Q YM , and P G = (z 3 − 1) 2 N 2 Ψ with N Ψ given in (4.51). Thus, all the ingredients have been computed in the Yang-Mills situation above, and we only need to work out the trace Tr(P G (Q G ) −1 ), which gives a lengthy result: where we can see that all poles s i,j , i, j = 1, . . . , 5 appearing therein, indicating the colorless structure of gravity amplitude. This complicated expression can be simplified by non-trivially imposing momentum conservation and Schouten identities. Applying the algorithm described in the appendix of [53] , for instance, we can simplify which agrees perfectly with the result given by KLT relation [46][47][48][49][50].

Gravity and n-point KLT Relations
More generally, for n-point amplitude, under the usual gauge-fixing z 1 = ∞, z 2 = 1, z n = 0, we wish to compute In order to write down the reduced Pfaffian, we need to compute the Pfaffian of a (2n − 2) × (2n − 2) matrix, which is quite complicated. Direct computation using the above formula is obviously very difficult, just like the direct computation of gravity amplitude by Feynman diagram. So we would like to follow the KLT formalism, and compute the gravity amplitude as square of Yang-Mills amplitudes.
An important property of the reduced Pfaffian is that, it can be expanded [3] as where α, β are permutations of labels 2, 3, . . . , n − 2, and S[α|β] is the S-kernel. The appearance of A YM n is a consequence of certain integrand summing over all (n − 3)! solutions of scattering equations in the original derivation, and in the companion matrix method, it corresponds to the trace of that integrand when changing variables to companion matrices. In any event, it is a constant, and can be dragged out of the trace.
Using this expression, we can expand one Pf Ψ in the gravity amplitude, Pf Ψ 2n×2n z 1α 2 z α 2 α 3 · · · z α n−2 ,n−1 z n−1,n z n1 . (4.60) The trace Tr(P (T z i )Q −1 (T z i )) for the set α gives A YM n (1, α, n−1, n), and the summation over permutations of α can be taken out of the trace, and we thereby arrive at the KLT relation 4 .

Six-Point Amplitudes
We proceed onto six-point amplitudes, i.e., n = 6. Using the standard gauge-fixing z 1 = ∞, z 2 = 1, z 6 = 0, Dolan-Goddard's polynomial form [37] of the scattering equations is given by We can thus define the ideal I = f 1 , f 2 , f 3 in the polynomial ring C[z 3 , z 4 , z 5 ]. The degree of ideal I is 6, so according to Bézout's theorem, it has 6 solutions, though it is not possible to obtain analytic expressions for these solutions, as already seen in the 5-point cases. Let us then consider the companion matrix method.
We generate the Gröbner basis for I in Lexicographic ordering z 3 ≺ z 4 ≺ z 5 . Analytically, the explicit expression of GB(I) is rather complicated, especially in the presence of so many parameters s ij in the ring. By varying the exponents to some high power, the polynomial reduction of the monomials z a 3 3 z a 4 4 z a 5 5 (with a i from 0 to some finite number, say 20) over GB(I) gives the monomial basis The polynomial reduction of z 3 B, z 4 B and z 5 B over GB(I) gives the companion matrices T z 3 , T z 4 , T z 5 , which are 6 × 6 matrices. Again, we need to compute P = P | z 3 →Tz 3 ,z 4 →Tz 4 ,z 5 →Tz 5 Q = Q| z 3 →Tz 3 ,z 4 →Tz 4 ,z 5 →Tz 5 , and the final amplitude is given by A 6 = Tr(P Q −1 ), without summing over all solutions of scattering equations.
Since the operations we need are multiplication of matrices, taking inverse or trace of matrices, so in principle it can be done analytically. However, the symbolic manipulation for n = 6 case is quite complicated, especially when taking the inverse of matrix Q and simplifying the tedious trace result in Mathematica, so we introduce random numeric kinematics -i.e., by Monte Carlo assignments of the parametres s ij -to get the final result. One will see that, as is customary with coefficient fields in polynomial rings, trying a few large prime numbers would suffice very quickly.

Scalar φ 3 theory
We can write the amplitude as

Prime Kinematic Strategy:
The idea is the following. Since we know for scalar φ 3 theory, the final result of A 6 should take the form  9 3 = 84 terms, which we denote as S ı , ı = 1, 2, . . . 84, and the amplitude is expanded as A 6 = 84 ı=1 c ı S ı , where c ı is either 0 or 1. To each physical pole we now randomly assign a prime number, i.e., we are working with the much simpler polynomial ring C[z] instead of C(s) [z]. In this case, the computation of Tr(P Q −1 ) is trivial within seconds, and the result as well as S ı 's are all numbers. Next, we shall find the solutions 84 ı=1 c ı S ı = Tr(P Q −1 ) for c ı 's. However, doing this by brute-force is impossible since there are 84 c ı 's and each one can take 0 or 1, so one would go through all 2 84 possibilities, which is far beyond any computational ability.
We therefore adopt the following strategy: instead of setting all coefficients to numbers, we can assign all physical poles to prime numbers except one pole. For example, we would leave s 345 , to detect first the coefficients of S ı 's which contains the pole 1 s 345 . Keeping one symbolic variable s 345 would extend the computation time of Tr(P Q −1 ) up to minutes, but it is still very manageable, while keeping two or more symbolic variables would make the computation of Tr(P Q −1 ) in Mathematica very hard for a laptop.
Let us see the above strategy in action. Setting the kinematics (coefficient variables) as, e.g., is also a part of A 6 . With the same procedure, we find that for physical pole s 123 , Putting all the above together, we therefore conclude that This prime-numeric method can be applied to all the cases of n = 6 amplitudes of scalar φ 3 theory.

Yang-Mills theory
For Yang-Mills theory, when n = 6, we meet the first "not so simple" gluon amplitude, i.e., the next-MHV amplitude, so it is worthwhile to verify the companion matrix method with this non-trivial example. To illustrate, let us consider the split helicity amplitude A YM 6 (g − 1 , g − 2 , g − 3 , g + 4 , g + 5 , g + 6 ), and choose the reference momenta as r 1 = r 2 = r 3 = k 4 , r 4 = r 5 = r 6 = k 3 , so that only 1,5 , 1,6 , 2,5 , 2,6 are non-zero. The object we want to compute is Here both the Jacobian |Φ| 126 126 and reduced Pfaffian Pf Ψ 12×12 are very complicated, so it is almost impossible to compute it analytically. As in the scalar φ 3 example, we can follow the semi-analytic procedure, and set the physical poles as some prime numbers, while keeping Ψ 12×12 (all i,j = i j , κ i,j = i k j and k i k j in Ψ) analytic. In this case, the ideal and Gröbner basis are just algebraic systems of polynomials with integer coefficients, while the elements of companion matrices are rational numbers. So the computation is very fast.
The Jacobian under the chosen gauge-fixing is where D Φ is polynomial in z 3 , z 4 , z 5 . The reduced Pfaffian together with the factor of scalar part give under the chosen gauge-fixing for some polynomial numerator N Φ . So we have Note that N Ψ originates from the Pfaffian of a 10×10 antisymmetric matrix, where by definition, each term in the Pfaffian is a product of five elements in the matrix. So each term in N Ψ is a product of five elements selected from i,j , κ i,j , k i k j , combined with a monomial of z 3 , z 4 , z 5 , for example, 2z 3 3 z 4 4 2,6 κ 1,2 κ 3,1 κ 4,1 κ 5,1 . Finally we can take the replacement P YM = P YM | z i →Tz i , Q YM = Q YM | z i →Tz i and compute the trace Tr(P YM (Q YM ) −1 ). It is given as Rewriting them as spinor products and applying the simplification algorithm for spinor expression, we get a one-page long result, which remarkably agrees with the known answers [54,55].

Seven-Point Amplitudes
The companion matrices T z i are simultaneously diagonalizable, and according to Stickelberger's theorem, the complex roots z i of ideal I are the vectors of simultaneous eigenvalues of the companion matrices T z i . Thus when they are evaluated in the diagonal form, the matrices P = P | z i →Tz i , Q = Q| z i →Tz i , P Q −1 are also diagonal, and it builds the one-to-one mapping between diagonal elements of (n − 3)! × (n − 3)! matrix P Q −1 and the integrand P/Q evaluated at the (n − 3)! complex solutions of scattering equations.
To demonstrate this, let us go through a 7-point example of scalar φ 3 theory. As usual, let us gauge fixing z 1 = ∞, z 2 = 1, z 7 = 0, and the amplitude is given by It can be checked directly that, each set of diagonal elements With the arithmetic result, it is possible to determine the terms appearing in amplitude by setting appropriate kinematics. In fact, in this example, we know that the result should be the sum where naively the summation is over all possible products of 4 physical poles s ı i , i.e., 14 4 = 1001 terms with c ı , ı = 1, . . . 1001, being either zero or one. By choosing 1001 different group of kinematics for physical poles, we get 1001 linear equations of (4.94), and solving them gives the c ı .
Indeed, the number of terms grows very fast with n in (4.94). The number of independent poles is n = (n−1)(n−2) 2 − 1 for massless theory, while the number of possible terms in the expansion is n n−3 . For n = 8, the number is 15504, and for n = 9 the number is 296010. So it is not very doable when n is large. However, for φ 3 theory, the number of color-ordered diagram is much smaller, and the counting is . So for n = 7, the possible terms appearing in (4.94) is 42 (an auspicious number). For n = 8, the number is 132, and for n = 9, the number is 429, etc. If we restrict to the 42 possible terms in (4.94), then it is possible to compute the coefficients c ı by choosing one set of kinematics.
One can let each physical pole be assigned a random prime number, and compute Tr(P Q −1 ) and then let Mathematica go through all 2 42 possibilities of c ı 's to find the summation 42 ı=1 cı sı 1 sı 2 sı 3 sı 4 = Tr(P Q −1 ). If the prime numbers in kinematic variables are distributed randomly in a very large scale, e.g., primes between 2 to 10000, then usually we can find one unique solution for c ı in the spirit of Egyptian fractions. This enables us to do one computation and fix all coefficients.
For example, let us compute With the kinematics shown in (4.90), we find the unique decomposition and agrees with the result given by CHY mapping rules [30,31].
There is a way to directly determine whether a certain term 1 sı 1 sı 2 sı 3 sı 4 is present in the result or not, by setting the kinematics s ı 1 = a, s ı 2 = a 2 , s ı 3 = a 4 , s ı 4 = a 8 , and others random primes not equaling to a. If this term exists, then the denominator has a factor a 15 . Again in the A 7 example, if we instead set is not a term in A 7 , while the 5 13 factor indicates that possible terms involving 1 s 12 s 67 s 123 must exist, which provides further information for detecting other existing terms. By this way, we can check all possible terms by setting kinematics for each one.
The number of solutions for scattering equations grows as (n − 3)!, while the companion matrix grows as (n − 3)! × (n − 3)!. When n = 8, we need to invert the matrix Q 120×120 , and at n = 9, the matrix Q 720×720 , etc. This sets the limitation on the computation of higher n.

Conclusions and Outlook
In this paper, motivated by the explanation of equivalence of different integrands in the CHY setup, we propose a new method using companion matrices, borrowed from the study of zero-dimensional ideals in computational algebraic geometry, to evaluate the integrand. One advantage of the method is that the rationality of final integral is obvious. Thus although our method may not be as efficient as the one proposed in [30,31,36], it does give a new angle to study the important problem of scattering amplitudes.
As shown in the plethora of examples, when the number of external legs grows, the analytic expression of companion matrix becomes harder. In fact, when n ≥ 6, the best way to do it is by assigning the kinematic variables to random prime numbers in order to reconstruct the analytic result. The salient feature of our method is that it is purely linear-algebraic, involving nothing more than finding the inverse and trace of matrices. The linearity of the trace, for example, was demonstrated to immediately lead to non-trivial identities in the amplitudes. Now, since the physical problem is very symmetric as can be seen by the polynomials given in (2.15) and (2.16), one is confronted with an immediate mathematical challenge. If we could analytically find, say by induction, the Gröbner basis and subsequent monomial basis for the polynomial form of the scattering equations in some appropriate lexicographic ordering, then one would find a recursive way to construct the companion matrix explicitly, much like the recursive construction of tree-level amplitude by using BCFW deformation [16,17]. Working out this construction is hard but worthwhile, as it would give explicit analytic results for the amplitudes and provide a deeper understanding of the CHY formalism.