BCJ numerators from differential operator of multidimensional residue

In previous works, we devised a differential operator for evaluating typical integrals appearing in the Cachazo–He–Yuan (CHY) forms and in this paper we further streamline this method. We observe that at tree level, the number of free parameters controlling the differential operator depends solely on the number of external lines, after solving the constraints arising from the scattering equations. This allows us to construct a reduction matrix that relates the parameters of a higher-order differential operator to those of a lower-order one. The reduction matrix is theory-independent and can be obtained by solving a set of explicitly given linear conditions. The repeated application of such reduction matrices eventually transforms a given tree-level CHY-like integral to a prepared form. We also provide analytic expressions for the parameters associated with any such prepared form at tree level. We finally give a compact expression for the multidimensional residue for any CHY-like integral in terms of the reduction matrices. We adopt a dual basis projector which leads to the CHY-like representation for the non-local Bern–Carrasco–Johansson (BCJ) numerators at tree level in Yang–Mills theory. These BCJ numerators are efficiently computed by the improved method involving the reduction matrix.

As discussed in [1,15,16], for the partial amplitudes that satisfy the BCJ relation, the BCJ numerators can often be obtained by acting with the famous Kawai-Lewellen-Tye (KLT) matrix, which is proposed in [17], on the partial amplitudes.An alternative method for computing the BCJ numerators, proposed in [18], is set up in the pure spinor formalism of string theory [19].In this approach, the worldsheet integral for an n-point string amplitude is decomposed into (n − 3)! compact integrals [20][21][22] and the BCJ numerators constructed in [18] are closely related to the coefficients of these compact ones.The decomposition often requires some help from the algebraic-geometry-based techniques, such as the Gröbner basis [23].There has also been progress in constructing the BCJ numerators from the Lagrangian that respects the color-kinematic duality in [24,25].
In this paper, we propose another method for constructing the BCJ numerators based on the Cachazo-He-Yuan (CHY) formula [26][27][28] and the polynomial scattering equations [29,30].The CHY form, which is a compact representation of scattering amplitudes, in general takes the following form at tree level for n external particles, where PT(12 • • • n) denotes the famous Parke-Taylor factor for this particular color ordering, and f i 's denote the scattering equations which are universal for the theories that have a string origin.Notice here we have already taken care of the SL(2, C) conformal symmetry by gauge fixing σ 1 → 0, σ n−1 → 1, σ n → ∞.The reduced Pfaffian Pf ′ (Ψ) is the only theory-specific part of the CHY form.Its explicit expression for a variety of theories can be found in the original papers [26][27][28].
The scattering equations f i are given by, while a simple transformation found in [29,30] takes the original equations to the polynomial ones, where 2 and notice that only (n − 3) equations are independent due to the SL(2, C) symmetry.The polynomial form of the scattering equations is what we actually use in this paper.The Jacobian of the said transformation is simply given by a Vandermonde determinant.Notice that instead of taking σ n−1 → 1, we have replaced it with an auxiliary variable σ 0 which makes these polynomials homogeneous and obviously at σ 0 → 1 the original expressions are recovered.
The CHY formula, as a beautiful representation of amplitudes, has computational advantages when the BCJ numerators and their double-copy are concerned, as demonstrated in [12,27,[31][32][33][34][35].This is a natural result of the string origin of the CHY formula and the multidimensional complex integral appearing in a CHY form is in parallel with the high-energy limit of an integral defined on the modular space of the string worldsheet [20,21,[36][37][38][39][40][41][42][43].
We find that the BCJ numerators can be collected via a decomposition of the CHY form into a basis of integrals, which conceptually is similar to the aforementioned decomposition that string worldsheet integrals undergo.The CHY forms can be evaluated using Integration Rules and the cross-ratio method at tree and loop levels [44][45][46][47].Recalling that the CHY form mathematically amounts to a multidimensional residue around the common zeros of the polynomial scattering equations, 1 the differential operator for evaluating such multidimensional residues proposed in our previous papers [48,49] comes in handy for evaluating and manipulating the CHY forms.Another method for computing multidimensional residues involving the Groebner basis or the H-basis is discussed in [50,51].A useful Mathematica package for computing such residues is given in [52].
We first summarize the conjecture for this differential operator here.Let g 1 , g 2 , ..., g k be homogeneous polynomials in complex variables z 1 , z 2 , ... z k .If their common zeros lie on a single isolated point p, for a holomorphic function R(z i ) in a neighborhood of p, there exists a differential operator D that computes the residue of R at p as follows, where the general ansatz for D reads, Here ∂ i = ∂ ∂z i , i = 1, . . ., k and r i 's are non-negative integers satisfying the Frobenius equation k i=1 r i = k i=1 deg(g i ) − k.The coefficients a r 1 ,r 2 ,...,r k are constants independent of z i 's, determined uniquely by two sets of constraints arising from: 1. the local duality theorem [53,54] and 2. the intersection number of the divisors D i = (g i ).
Detailed discussions on these constraints can be found in [48] and this conjecture, although first proposed for the CHY form, is widely applicable to any such residues.In the case of the tree-level CHY formula, in particular, we have (n − 2) integration variables σ i including the auxiliary one σ 0 and also (n − 2) polynomial factors in the denominator.Among them (n − 3) factors are the polynomial scattering equations (1.2) while the last one which we usually denote as h 0 is theory-specific.When h 0 = σ i , we call the CHY form of the prepared form and scrutinize such cases in [49], where an analytic solution of the corresponding differential operator is obtained.
The decomposition of the differential form in the CHY integral is translated to the decomposition of its corresponding differential operator, and the latter is carried out systematically using what we will define in later sections as the reduction matrix.A crucial observation in our approach is that the number of independent coefficients appearing in such a differential operator is always (n − 3)!, regardless of the order of the operator.As will become clear in detailed discussions, this observation allows us to develop a method that relates higher-order differential operators with lower-order ones, through the reduction matrices.Eventually all the differential operators are efficiently decomposed into a basis in which every operator is of the same order.This method is completely theory-independent and applies to all that admit a CHY representation of amplitudes.For Yang-Mills in particular, our approach beautifully preserves the gauge invariance of a n-point amplitude in (n − 2) external legs manifest in every intermediate step. 2t loop levels, we expect our method to apply straightforwardly to those theories in which the loop-level CHY forms are already discovered; these loop-level CHY forms are scrutinized in [35,[55][56][57][58][59][60][61][62][63].Our method has also a neat byproduct that is a highly efficient evaluation of a general CHY form, both at and beyond tree level.As shall be discussed in detail, the factorized form of the CHY integrand can always be rendered a prepared form by a sequence of reduction matrices, while the prepared form is analytically evaluated at tree level and one loop in our previous work [49].This paper is organized as follows.In Section 2, we review the concept of dual basis and it allows us to present a conceptually straightforward procedure for finding the BCJ numerators.In Section 3, we study in detail the systematic method that realizes the construction of BCJ numerators using the differential operator we have just reviewed.We also include in this section discussions and comments on the number of independent coefficients being always (n − 3)!.In Section 4 we further illustrate this proposed method with 4-and 5-point YM tree amplitudes.For the 5-point case we focus on two characteristic terms in this section, while the full results are given in Appendix A. In Appendix B we present the explicit expressions for the reduction matrices appearing in the 6-point amplitudes.

Review of BCJ Numerator
A number of theories with color-kinematic duality, such as Yang-Mills theory, nonlinear sigma model (NLSM) and bi-adjoint scalar theory, are shown to respect both the Kleiss-Kuijf (KK) relation [64] and the BCJ relations [1][2][3][4][5][6][7][8][9][10][11].The fundamental BCJ relations are given by, As studied in [1,15,16], exploiting the KK and the BCJ relations, any partial amplitude in the bi-adjoint scalar theory can be decomposed into the minimal basis whose dimension is (n − 3)!.The bi-adjoint scalar partial amplitudes can be packed in the matrix m αβ where the indices denote the color-orderings.This matrix is also the inverse of the KLT matrix [17], with which a general partial amplitude in any theory with the KK and the BCJ relations can be written as Here N 1,β,n−1,n are functions of kinematic variables only and are called the BCJ numerators in the minimal basis, although strictly speaking, they may be rational functions, instead of polynomials.The number of non-vanishing BCJ numerators is (n − 3)!.

BCJ Numerators from CHY Form
The integrand of a general CHY form (1.1) can be decomposed into a basis of differential forms.As mentioned above, the dimension of this basis is (n − 3)! and we are allowed to choose the basis elements to be the following, where PT(1βn − 1) is the Park-Taylor factor for the order (1, .
And β is an element of the permutation group Suppose the coefficient of each basis element is C 1βn−1n and then the CHY form for the color-ordered tree amplitude can be rewritten as, PT(1αn − 1) (2.4) For notational brevity, let us denote the CHY-like integral, that is an integral around the solutions of the scattering equations, as an inner product.Recall that the bi-adjoint scalar amplitude m αβ also has a CHY representation [31] and hence it can be denoted as an inner product as follows, Similarly, in this notation (2.2) reads, which holds for any color-ordering (1, α, n − 1, n).
Now we adopt a dual basis projector PT 1αn−1n of the Parke-Taylor factor PT(1, α, n− 1, n), with which we show that the coefficients C 1βn−1n are precisely the BCJ numerators.This projector, by definition, satisfies the condition where the dual basis projector can also be obtained through the KLT matrix S αβ as follows, That is to say, the projectors are always linear combinations of the Parke-Taylor factors of different color orderings.Since (2.6) is true for any color ordering, it also holds for their linear combinations.Therefore, we can safely replace the Parke-Taylor PT(1αn − 1n) on both sides of (2.6) with its dual projector PT 1αn−1n and this leads to, (2.9) As a demonstrating example, we consider the four-point NLSM amplitude at tree level, whose CHY form is given by, 134 PT(1234) (2.10) One the one hand, let the color ordering (1234) be the basis of the one-dimensional bi-adjoint matrix m (1234)(1234) = Replacing the Parke-Taylor in (2.10) with this projector, we obtain, (2.12) This agrees with the previous result.

Ambiguity of BCJ Numerator
So far we have worked in the minimal basis whose dimension is (n − 3)! for n external lines; another often-used basis, the Del Duca-Dixon-Maltoni (DDM) basis [65], is of (n − 2)! dimensions.In the latter, the CHY form of the tree amplitude can be written as [35] A n (1α where the functions N ′ 1β ′ n denote the BCJ numerators in the DDM basis, which is not to be confused with N 1βn−1n previously discussed in the minimal basis.The number !. Yet due to the amplitude relations of the double-colored scalar field, these differential forms are not all independent.Besides, the BCJ numerators are not unique and their ambiguities arise from these amplitude relations.The BCJ numerators constructed in the minimal basis can always be turned into those in the DDM basis and vice versa, by adding or removing pieces proportional to these relations.
In the language of the CHY form, the BCJ-like amplitude relations [35,66] are given as the following, where the summation is taken over all possible permutations keeping the relative orders, namely, the order-preserving (OP) permutations.ξ i denotes the position of the leg i and denotes the ideal generated by the scattering equations Here a n−1 = n − 1 and the other a i and b j are drawn from the set {2, 3, • • • , n − 2}.
The number of the independent relations here are (n − 3) × (n − 3)!.Hence, the total number of independent differential forms is cut from (n − 2)! down to (n − 3)!, the dimension of the minimal basis.
The coefficients in each BCJ relation above form a vector E r , r ∈ [2, n − 2], whose components are indexed by γ that is an element of the order-preserving shufflings.The components of them are The module of E r over rational functions of s ij can fully characterize the ambiguity of the BCJ numerators.Let N 0 be a (n − 2)!-dimensional vector whose components are labeled by the color-orderings of the form (1β ′ n), β ′ ∈ S n−2 .For an ordering belonging to {(1βn − 1n)|β ∈ S n−3 }, the respective component is simply N 1βn−1n , while for those orderings outside this set, their respective components are all zero.Then any BCJ numerators are of the following form, where R i 's are arbitrary rational functions of s ij .These rational functions are of the same mass dimension, of which the value is theory-specific.
Let us once again take the 4-point NLSM amplitude as an example.In this case the color-orderings in the DDM basis are {(1234), (1324)} while in the minimal basis we only have the first element.Then we have where we have used the on-shell condition in N 0 .By choosing R 2 = −s 13 , we recover the result in [34] that is

BCJ Basis from Differential Operator
As demonstrated above, with the dual basis projector, the BCJ numerators can be procured through a standardized procedure.However, in practice, computing the inner product of the projector and the Pfaffian, before decomposing the latter, can be quite cumbersome for higher multiplicities.Hence, instead of directly evaluating such inner products, it is natural to decompose the Pfaffian first into the minimal BCJ basis PT(1βn − 1n) (β ∈ S n−3 ), which trivializes the inner products.In this section, we propose a systematic method for such decompositions with the help of the differential operator reviewed in Section 1.
Consider the following differential form, where is the homogenized polynomial scattering equation of degree j and σ 0 is the auxiliary variable as introduced in (1.2).The last factor in the denominator h 0 and the numerator P (σ) are both polynomials in σ's.Recall the conjecture mentioned in Section 1 that, for such a (n−2)-form Γ we can associated with it a differential operator D Γ , of which the action computes the residue corresponding to Γ, that is Henceforward, at the level of CHY-like integrals (integration of a differential form around the solutions of the scattering equations), each vector of the minimal basis, namely PT(1βn − 1n), is equivalent to a differential operator D β .
Likewise the Pfaffian can also be replaced by differential operators.Schematically, the expression of the (reduced) Pfaffian appearing in a CHY-integral, although theoryspecific, is always of the following form, where the numerators N h 0 are nothing but bilinears of the external momenta k i and, where possible, the helicities of external particles ǫ i as well.The denominators are factorized polynomials of the integration variables σ i .After homogenization of all the factors in the denominators with the auxiliary variable σ 0 , each term on the right-hand side is translated to a differential operator D h 0 and the residue corresponding to the Pfaffian is simply, The decomposition of the Pfaffian into the minimal basis is thus translated into the decomposition of each differential operator D h 0 associated with a given term in the Pfaffian into a basis formed by those operators D β associated with the minimal vectors.Throughout this section, instead of decomposing the functions directly using the scattering equations, we will work towards replacing the Pfaffian with a linear combination of the differential operators corresponding to the BCJ vectors.
Recall that the gauge fixing of the SL(2, C) symmetry, i.e (σ 1 , σ n−1 , σ n ) → (0, 1, ∞), changes the degrees of the denominators both in the Parke-Taylor factors and in the Pfaffian, resulting in that their respective differential operators are of different orders.
As we are about to discuss in detail, we can construct a basis from the lowest-order differential operators that trivializes the decomposition, and more importantly, relate all the higher-order differential operators to these lowest-order ones by what we would like to coin as the reduction matrices as long as the denominators are factorized, which is naturally true in many theories.Hence, all the differential operators we encounter can be treated on an equal footing, although they are of different orders at the beginning.The factorized forms of the denominators, as we will see shortly, also allow us to express right away the a-coefficients appearing in a differential operator with deg(h 0 ) > 1 in terms of the reduction matrices and the a's in another differential operator with deg(h 0 ) = 1.The latter is of the prepared form defined in [49], for which the solutions of a's are derived in the same paper.Moreover, for Yang-Mills in particular, with the help of the gauge-invariant reduced Pfaffian Pf ′ Ψ 1n (where the 1st and the n-th external lines are fixed) worked out in [67], the results expressed in terms of these reduction matrices are always manifestly gauge invariant in the (n − 2) unfixed external lines.

Canonical Coefficients in Differential Operator
The differential operator D Γ associated with a given form (3.1), as reviewed in Section 1, is of the following form, where n−2 ∂ s 0 0 and {s} ord(D) denotes the solutions to the Frobenius equation for all the theories that have the CHY representations, we can always solve the local duality conditions corresponding to these polynomials first, without any knowledge of the specific theory which only comes into play through the factor h 0 .These conditions read, where q j (σ) scans over all the monomials in σ's of the degree deg(q j ) = ord(D Γ )−j.For a given number of external particles, the number of the a-coefficients and the number of local duality conditions both grow as the degree of h 0 increases, and it might not be obvious what to expect of the solutions of (3.6) at the first glance.However, a curious observation is that the number of independent a-coefficients after solving the equations (3.6) is always (n − 3)!, independent of the degree of h 0 .3 Although we are in principle free to choose whichever (n − 3)! a-coefficients in a given differential operator as our parameters, for the method we are to present shortly, there is one choice that is particularly convenient, and we would like to call these acoefficients chosen as the parameters the canonical coefficients.Suppose the degree of the polynomial h 0 is d h 0 and the order of the relevant differential operator is, We define the canonical coefficients as the ones of the following form, All the other a-coefficients are related to these canonical ones through scattering equations and can be obtained straightforwardly.In the rest of the section, we focus on the canonical a's only.
The minimal basis formed by the Parke-Taylor factors PT(1, β, n − 1, n) , although commonly used, is not a preferable choice for our method.Instead, we would rather construct a new basis from the Parke-Taylor factors such that this new basis satisfies the following conditions: 1. the denominators of the elements are all of the same degree , which is the lowest degree possible for the Parke-Taylor factors and leads to the order of the corresponding differential operator being m 0 = (n−1)(n−4)/2; 2. this basis makes the vector of the canonical a-coefficients a unit vector of dimension (n − 3)!.That is, for a vector b j (j = 1, 2, • • • , (n − 3)!) in this new basis, the canonical a-coefficients appearing in the corresponding differential operator D b j form the unit vector as follows, ) . (3.9) For this reason, we call this basis the unit basis and it is obviously equivalent to the minimal basis, up to a linear transformation.Here we only present the definition of this new basis.The construction of this basis is relatively straightforward which we will illustrate with examples later.
In the unit basis, the decomposition of a given differential operator D (m 0 ) of order m 0 becomes trivial.With all the non-canonical a-coefficients expressed in terms of the canonical ones, the operator D (m 0 ) can be written as the following, where the non-canonical a-coefficients are expressed in terms of the canonical ones as the following, and therefore The coefficients c s 2 •••s n−3 s 0 γ in the summation are obtained by solving the local duality constraints related to the scattering equations (3.6) and can be functions of the external momenta only.Note that these coefficients are universal for all differential operators of the same order, for a given n.Hence we have, where M = (n − 3)! and e i is the M-dimensional unit vector with the i-th component being 1.In the second line we have used Thus we have decomposed a general operator D (m 0 ) in the unit basis D b j .These canonical coefficients obtained through the decomposition can be viewed as the "BCJ numerators" in this unit basis, and are equivalent to those in the minimal basis mentioned in the previous section, up to linear transformations.The decomposition of higher-order operators requires the help of the reduction matrix that we will discuss later in this section.

Discussion on the Dimension of BCJ Basis
Before we proceed to the rest of our proposal, let us pause for a moment and make a few comments on the canonical coefficients we have just defined.
The number of the canonical coefficients in any differential operator with n external particles, as we have observed, is always (n − 3)!, independent of the order of the operator.There are three other numbers of the same value (n − 3)!.
1.The number of independent partial amplitudes at tree level in a theory that respects the KK and the BCJ relations.[68,69], where

The dimension of the twisted cohomology
Here X is an open set defined as the space C n−3 with the discontinuity points removed, while these discontinuity points are those where the worldsheet coordinates coincide with each other.The branches of the multi-valued function |σ ij | α ′ s ij , which appears in the worldsheet integrals of string amplitudes, can be conveniently characterized by the open set X.
3. The number of solutions of the tree-level scattering equations, as verified in [26,29,30].
Together with the number of canonical coefficients, these four quantities of the same value (n − 3)!, summarized in Figure 1, are universal for all the theories in which the amplitudes respect the KK and the BCJ relations.And here we try to understand this curious coincidence of the four quantities.
It is usually true that a theory in which the amplitudes satisfy the KK and the BCJ relations has a "string origin" in the following sense.There exists a theory defined on the string worldsheet in which the tree amplitudes in general take the following form, where h 0 is a monomial of σ ij with degree n−2, which is a result of conformal invariance, and N h 0 are functions of external momenta and helicities.When taking the limit On the string theory side, or the "stringy theory" side should the theory be not one of the known string theories, the Integration-by-Parts (IBP) relations demand that the addition of a total derivative term ∇ ω ξ(σ i ) to the integrand in (3.15) does not alter the worldsheet integral.That is, we do not need to consider all the differential forms that can possibly enter (3.15).Rather, it suffices to take into consideration only those that are not related by total derivatives.The collection of these IBP-inequivalent differential forms is precisely the twisted cohomology H (n−3) (X, ∇ ω ).
The basis of this twisted cohomology is carefully analyzed in literature [20][21][22].A convenient choice of basis is Note that the kernel and the image of ∇ ω , which define the twisted cohomology, are independent of the parameter α ′ .Therefore the twisted cohomology remains the same when going from the string side to its field theory limit.Moreover, notice that this basis of differential forms is identical to the minimal BCJ basis on the field theory side, and therefore the respective numerators on both sides should also admit a one-to-one correspondence.
On the string side, the integrand is cast into the aforementioned basis using the IBP relations whereas on the field theory side, where instead of worldsheet integrals we have the CHY integral, the integrand is cast into the BCJ basis using the scattering equations.We have observed in examples, 4 that the IBP relations used on the string side can be connected to the scattering equations on the field theory side.Since the BCJ basis on the field theory side can be viewed as the equivalence classes of the scattering equations, the identical dimensions of the BCJ basis and the twisted cohomology seems to suggest that the ideal of those IBP relations used for this purpose can be related to the ideal of the scattering equations.
The relation between the dimension of the twisted cohomology and the number of solutions of scattering equations can be understood from a different angle.The tree-level scattering equations have (n − 3)! solutions and on each solution, the twisted cohomology is simply following differential form where σ (0) i is the position of the solution.Namely, the dimension of the twisted cohomology on a single solution of the scattering equations is 1 while the total dimension of the twisted cohomology is then the number of solutions, given that the differential form (3.17) is of the dlog form.The dlog form can be guaranteed.Although the degree of the denominator h 0 in the integrand (3.15) appears to be (n − 2) after the gauge fixing of conformal symmetries, the integrand can always undergo a partial fraction decomposition which results in a sum of terms whose denominator are one degree lower.
Last but not the least, for any given order of the differential operator, the number of the canonical coefficients, which are solely determined by the scattering equations, remains the same and appears completely unaware of the factor h 0 .This indicates that, if the degree of h 0 is fixed, there are only (n − 3)! independent differential operators 4 Consider the following term in the 4-point amplitude for example.The integral in the string-like theory and the CHY form are respectively The relation which is used to reduce this integrand into the basis integrand PT(1234 are respectively: the IBP-induced equation ∇ ω σ2−1 = 0 on the string side; the scattering equation = 0 on the CHY side.In the limit α ′ → ∞, the IBP-induced equation becomes the scattering equation.In higher multiplicity amplitudes, the situation is similar.Similar ideas are discussed in [41].We will discuss more systematically in future project.D h 0 .This is equivalent to that there are only (n − 3)! independent differential forms in a CHY integral, when the degree of h 0 is fixed.
The degree of h 0 can always be turned into (n − 3) by partial fraction decompositions, after all, and from this derives the legitimacy to construct the basis of independent differential forms with fixed-degree h 0 's only.But to keep gauge invariance nicely manifest when decomposing the differential forms is a different matter.This shall be accomplished via the reduction method we are about to propose, where we do need to work with h 0 of different degrees.

Reduction Matrix
As we have seen already, the aforementioned unit basis of differential operators D b j is advantageous for decomposing any given differential operator D (m 0 ) of the same order.Now we move on to the decomposition of differential operators of higher orders into the same unit basis.
First and foremost, let us make sense of the idea of decomposing a higher-order operator into a basis spanned by lower-order operators only.Consider a differential operator D (m) Γ of order m and recall that it is associated with a (n − 2)-form of the form (3.1).Suppose that the last factor in the denominator h 0 is factorized as h 0 = h ′ 0 (σ)q(σ) where both h ′ (σ) and q(σ) are homogeneous polynomials in σ's, and that the degree of the numerator deg(P ) m − deg(q).Then the residue at the origin computed by D (m) Γ can be equated to one computed by a lower-order operator D(m−dq) where d q = deg(q).That is, at the level of residues, for which such operators are designed, we have . Note that the factor h 0 we encounter are always completely factorized with every factor of degree one.On the other hand, the numerator P (σ) is the remainder of the Vandermonde determinant, arising from the polynomial scattering equations, with part of it canceled by the Parke-Taylor associated with the color ordering.By a simple counting of degrees, one easily observes that the degree of the numerator is (n − 1)(n − 4)/2 for n external lines. 5As a result, we can always choose q(σ) such that m − d q = m 0 and D(m 0 ) Γ is readily decomposed as discussed above.
Now we investigate in detail the relation between the higher-and lower-order op-erators.Consider the meromorphic forms Γ and Γ ′ , the latter taking the form of the former only with h 0 substituted by h ′ 0 , and their corresponding differential operators are D .When P (σ) is an arbitrary polynomial of degree deg(P ) m − d q , the self-consistency of such operators requires, Plugging in the solutions of the respective non-canonical a-coefficients on both sides and expressing both differential operators in terms of their canonical coefficients only, the equation above becomes, where d = deg(h 0 ) and the operators D is defined in (3.12).For this equation to hold for an arbitrary P (σ), the coefficients of the surviving derivatives . . .
where we name the matrix M n,m q(σ) the reduction matrix whose elements are just functions of s ij .Note that the reduction matrix depends only on the factor q(σ) and the orders of the differential operators while it knows nothing about the specific expression of the factor h ′ 0 .Moreover, its dependence on the operator orders is quite trivial, as we will soon observe in examples.This reduction operation can be performed order by order, as long as the polynomial h 0 can be factorized into a product of degree-one polynomials.As a result, the canonical coefficients in D (m) Γ can eventually be related to those in an operator of order m−(d−1).
where each q (r) is a σ ij factor of the original h 0 .
So far we have established the relations between the canonical coefficients belonging to two operators D (m) Γ and D (m ′ ) Γ ′ .These relations hold for an arbitrary polynomial numerator and thus it is always legal to relate the canonical coefficients of two differential operators, as long as h 0 and h ′ 0 differ by only a factor.The operator , although of a lower order, should not be confused with D(m 0 ) Γ that gives rise to the original residue associated with Γ, which is the main goal of this section.
As we are to discuss in detail shortly, the reduction method is employed in two stages of handling an operator D (m) Γ with m > m 0 .In the first stage we construct the differential operator D(m 0 ) Γ equivalent to the original operator and the coefficients of the former can be expressed conveniently with the correct reduction matrix.At this point, we have already completed the decomposition while the canonical coefficients to be determined are those in an operator of order m 0 .By virtue of the factorized forms of the denominators in CHY forms, (3.21) applies beautifully and in the end, we only need the canonical coefficients corresponding to a prepared form, of which an analytical solution is derived in [49]. 6o construct this operator D(m 0 ) Γ , let us write down first the ansatz for it as follows, where the coefficients are denoted as b s2 •••s n−2 s0 to signify that this operator D(m 0 ) Γ is of a somewhat different nature than those whose coefficients are denoted as a's.To be precise, these b-coefficients are not directly obtained by solving the local duality and the intersection number constraints; rather, they are related to the usual a-coefficients as we are to present shortly.

Recall that we have
Γ in the sense that the actions with them result in the same residue for Γ when deg(P ) = m 0 , that is, on the left-hand side, we obtain where in the second equal sign we have used the fact that the degree of P (σ) is n − 3 and its derivatives of higher or lower orders do not survive.Likewise, the right-hand side gives,  In practice, as we will illustrate with examples, the highest order of the differential operator that we encounter is n − 2 for n particles.Recalling that m 0 is corresponding to the differential operator with d h 0 = n − 3, (3.26) simplifies to where the right-hand side can always be written as linear combinations of the canonical coefficients {a γ(0)•••γ(n−4)(n−3) }.Let q 1 (σ) be a factor of h 0 in Γ and deg(q where M n,m 0 +1 , with M n,m 0 +1 q 1 denoting the reduction matrix relating the two operators of orders (m 0 + 1) and m 0 respectively and c m 0 is the constant integer obtained in (3.27) that depends only on m 0 for canonical coefficients.Notice that {a γ(0)γ(1)•••γ(n−4)(n−4) } are the canonical coefficients in the operator D (m 0 ) Γ .In general, the degree of the numerator appearing in a n-point CHY expression can only be as low as m 0 , and consequently, we are not allowed to reduce the order of differential operators any further than that.This echoes with our choice for the basis for the decomposition.Nonetheless, it is always legal to adopt (3.21) to work out the explicit expressions of these canonical coefficients quickly.

Gauge Invariance of Yang-Mills BCJ Numerators
In the last section, we have laid out a general description of our method for computing the BCJ numerators in the minimal basis and also discussed in detail the key ingredient of this method that is the reduction matrix.In this section, we take the Yang-Mills 4and 5-point tree amplitudes to illustrate how the story plays out.

Yang-Mills 4-Point Amplitude
For a color-ordered tree amplitudes with four external particles in YM, its corresponding CHY-form reads, where we have taken care of the SL(2, C) symmetry by gauge fixing (σ 1 → 0, σ 3 → 1, σ 4 → ∞) and are left with only one integration variable σ 2 and its respective scattering equation f 2 = 0.The reduced Pfaffian for Yang-Mills can be repackaged in the following way as proposed in [67], where we have adopted the similar notations used in [67] as follows, The external leg 1 and leg 4 are fixed in this way of writing in the sense that they only enter through ǫ 1 and ǫ 4 in the F i 1 •••ir ǫ 1 ,ǫ 4 functions.For a term in (4.2), if it contains F i , it vanishes as ǫ i → k i ; i.e. a term in (4.2) is manifestly gauge invariant in the external leg i if the quantity F i is present in one of the traces appearing in this term.Hence, the terms painted black in (4.2) are already manifestly gauge invariant in both leg 2 and leg 3 at the level of the reduced Pfaffian and throughout our computation this manifest gauge invariance is preserved.As for the red terms, each of them contains at least one C ii , which vanishes as ǫ i → k i on the scattering equations.In this sense, they are also manifestly gauge invariant in the unfixed legs.The gauge invariance in the two fixed external legs, leg 1 and leg 4 in our case, is not manifest in a single term in (4.2), but can be verified when all the terms are summed up.
Just like in the case of NLSM, for four external lines, the dimension of the minimal basis is 1 and the basis is simply PT(1234).The unit basis in this case is but PT(1234) multiplied by a normalization factor, namely, b 1 = s 12 (s 12 + s 13 ) s 13 PT(1234).(4.5) The order of the corresponding differential operator is 0, i.e this operator is the identity.But let us still formally denote this identity operator as After gauge fixing and homogenization, each denominator in (4.2) can only be a character drawn from the set {σ 2 , σ 2 − σ 0 , σ 2 (σ 2 − σ 0 ), (σ 2 − σ 0 ) 2 }.For the last two characters, the reduction matrix M σ 2 −σ 0 renders them one degree lower and the same as the first two.When the denominator is σ 2 , the integral is of the prepared form and the corresponding a 00 can be read off.The solution of a 00 in the prepared form with σ 2 being the denominator is denoted as A σ 2 below.When the denominator is σ 2 − σ 0 , we may replace it with σ 0 (σ 2 − σ 0 ) and employ the reduction matrix M σ 2 −σ 0 yet again, making the integral another prepared form with the remaining factor σ 0 , for which the solution of a 00 is denoted as A σ 0 .Explicitly we list the values of these solutions and the reduction matrices, Thus each term in (4.2) is replaced by its corresponding operator as follows, We have denoted as the replacement of the original term in the Pfaffian by its corresponding evaluation in terms of the reductive matrices and the solutions of prepared forms.As a result, the reduced Pfaffian (4.2) evaluated at the solution of the scattering equations reads, Pf ′ (Ψ 14  14 where we have used implicitly the amplitude relations to cast the result into the minimal basis. Following the discussion in Section 2.1, the BCJ numerator in the minimal basis can be extracted from the Pfaffian by taking inner product with the dual basis projector PT 1234 , which simply picks out the coefficient of PT(1234) in (4.9), namely, Notice that the unphysical pole s 1,23 appearing in the first line above drops out in the second line, after algebraic manipulations.This coefficient, up to amplitude relations, is the sum of all the A σ appearing in (4.8) weighted by the relevant reduction matrices and hence, the solutions of the canonical coefficients we encounter in the evaluation of the Pfaffian can be viewed as "the BCJ numerator" in this unit basis. 7o check our result N 1234 , we compare it with the BCJ numerator in the DDM basis {PT(1234), PT(1324)}.The BCJ relations discussed in Section 2.2 render the first two terms in (4.10) conveniently in the DDM basis as the following, The last four terms in (4.10) can also be molded into the DDM basis by plugging in the explicit definitions for the F -functions as introduced in (4.3) and we obtain It is easy to recognize that the last two lines above, namely the terms without any denominator contribute to the coefficient of PT(1234) while the first two lines, namely all the terms with a factor s 12 s 13 contribute to the coefficient of PT(1324).Unfortunately, in the DDM basis, we have lost the manifest gauge invariance in leg 2 and leg 3.
There is another bonus feature of the results we have obtained.Since the BCJ numerators have the ambiguity we have discussed Section 2.2 and can be rearranged, we are allowed split the term s 12 ǫ 1 • ǫ 4 ǫ 2 • ǫ 3 in (4.12) into two halves and use the BCJ relations to place them in the coefficients of PT(1234) and PT(1324) respectively.Together with (4.11), we have come to another set of BCJ numerators, slightly different from the aforementioned results.It is then easy to verify that these numerators respect the "crossing symmetry" under the swapping of leg 2 and leg 3. 8

Yang-Mills 5-Point Amplitude
In this subsection, we present another example, the five-point tree amplitudes in Yang-Mills theory.We use the same method to deal with this example and the complication compared with the four-point case is more technical than conceptual.In this example, it will become more clear how to deal with different structures appearing in the Pfaffian such that the manifest gauge invariance in the (n−2) unfixed external lines is protected throughout the computation.
Like always, we choose the gauge (σ 1 → 0, σ 4 → 1, σ 5 → ∞) and the CHY-form for the color-ordered five-point tree amplitude in Yang-Mills reads, The reduced Pfaffian written with the notations introduced in (4.3) takes the following form, where we have fixed the external line 1 and line 5 in Notice that every term in the first line contains F 2 , F 3 , F 4 in one of the (• • • ) brackets and therefore they are manifestly gauge invariant in all the unfixed external legs.As for the rest painted red, they all contain at least one C ii and the C ii factor only is manifestly gauge invariant in the i-th leg on the scattering equations.That is, on the scattering equations, every term in (4.14) is manifestly gauge invariant in the three unfixed legs.The gauge invariance in the two fixed legs, although not manifest in each term, is nevertheless protected when all the terms in (4.14) are summed up.
The minimal basis for five points is chosen as {PT(12345), PT(13245)} and it is straightforward to construct the unit basis that is a linear combination of the minimal basis.We merely need to solve for the canonical coefficients defined as (3.8) corresponding to the elements of the minimal basis respectively and find the transformation denoted as B 5 which makes the vector of the canonical coefficients a unit one.Then the unit basis is simply, The unit basis vectors take the following explicit forms, The Parke-Taylor factors, although they seemingly have denominators of degree 3 after gauge fixing, can be split into a sum with lower-degree denominators.For instance, and likewise for PT(13245).Hence, each term of b 1 and b 2 has only degree-2 denominators and their corresponding differential operators are consequently the second-order ones.Then it is easy to verify that the canonical coefficients, (a 011 , a 101 ) in particular, in the differential operators corresponding to b 1 and b 2 form unit vectors (1, 0) and (0, 1), while the non-canonical coefficients are expressed in terms of the canonical one by solving the equations (3.6).More explicitly, we write down the expressions for these differential operators as follows, D b 1 = 1, 0, − s 12 s 14 , − s 12 (s 12 +s 14 +s 24 ) 2s 14 (s 12 +s 23 +s 24 ) , − s 2 12 +(s 14 +s 23 +s 24 )s 12 +s 14 (s 13 +s 23 ) 2s 14 (s 12 +s 14 +s 24 ) , − s 12 +s 14 +s 24 2(s 12 +s 13 +s 23 ) b 2 = 0, 1, − s 13 s 14 , − s 13 (s 12 +s 23 +s 24 ) 2s 14 (s 12 +s 14 +s 24 ) , s 12 (s 14 −s 13 )+s 14 s 23 −s 13 s 24 2s 14 (s 12 +s 23 +s 24 ) , s 12 +s 23 +s 24 2(s 12 +s 13 +s 23 ) Now we are ready to compute the differential operators associated with the terms in (4.14) and decompose them into the unit basis (4.18).Without loss of generality, we only discuss the following two characteristic terms in detail while the rest can all be treated in like manner, where we have substituted the definition of C 44 while leaving the F -functions intact.
The prefactor multiplying the reduced Pfaffian (4.14) is given by, where J(σ) is the Vandermonde matrix like in the 4-point case and this quantity is a monimial of degree 2 in σ ij 's in the gauge of choice.
Consider the term I 1 first.Its denominator is of degree 3 and the corresponding differential operator is therefore of order 3, which we denote as D (3) σ 2 (σ 3 −σ 0 ) 2 with the denominator homogenized.The coefficients present in this differential operator are the following, {a 012 , a 102 , a 111 , a 201 , a 021 , a 003 , a 210 , a 120 , a 300 , a 030 } , ( where the first two coefficients colored in red are the canonical coefficients on which we focus throughout the computation while the rest can all be written as linear combinations of the canonical ones.For brevity we will not display the explicit expressions for the non-canonical coefficients in this section, remembering they are obtained by solving (3.6).The four black coefficients in the middle do have a contribution to the final residue while the last four in gray do not at all, due to the fact that the degree of the numerator N 5 is no more than 2.
The differential operator D σ 2 (σ 3 −σ 0 ) 2 , as before, has to be identified with a secondorder operator, which is easily achieved by using the reduction matrix M 5,3 (σ 3 −σ 2 ) associated with the factor (σ 3 − σ 2 ). 9 Denoting the resulted second-order operator as D(2) σ 2 (σ 3 −σ 0 ) , the coefficients of D(2) σ 2 (σ 3 −σ 0 ) are related to those in (4.22) through (3.27), and therefore the canonical coefficients in D(2) σ 2 (σ 3 −σ 0 ) are simply {2a 012 , 2a 102 }.Recall that the canonical coefficients in a usual second-order differential operator of the form (3.5) are {a 011 , a 101 } and the reduction matrix M 5,3 (σ 3 −σ 2 ) connects the two sets of the canonical coefficients as follows, 2M 5,3 where the subscript of each vector denotes the denominator their respective differential operator is associated with.Notice that the first vector belongs to the second-order D(2) σ 2 (σ 3 −σ 0 ) , the intermediate one belongs to the original third-order D σ 2 (σ 3 −σ 0 ) 2 while last one another second-order operator corresponding to the same denominator.The reduction matrix reads, Let us stress once again that the matrix M 5,3 σ 3 −σ 2 does not depend on the explicit form of the remaining factor in the denominator that is σ 2 (σ 3 − σ 0 ) in this case.The canonical coefficients {a 011 , a 101 } can be further rewritten as where the right-hand side is the canonical coefficients associated with a prepared form with the denominator σ 2 .Their values are easily read off from the analytical solution derived for prepared forms in our previous work [49], which reads Consequently, the differential operator that computes residue associated with I 1 reads, where the operators D ijk are defined in (3.12) while the coefficients of the derivatives can be read off from (4.18).The residue associated with this term I 1 is obtained as the following, Recall that the differential operator (4.27) is already decomposed into the unit basis (4.18), namely, The integrand I 1 we start with is then decomposed into the unit basis (4.16) with the same coefficients.To cast this term into the minimal BCJ basis, we only need to include the transformation B 5 and we obtain, where the inner products need not to be calculated by brute force, rather, by the definition of dual basis, they can simply be read off from (4.29) as the coefficients of the respective Parke-Taylor factors.
Now we move on to the second term I 2 .This integrand I 2 is manifestly gauge invariant with the unfixed legs only as a whole, while the three terms of it are not.To preserve the gauge invariance, we first employ the reduction matrices for each term in the C 44 factor and sum the three matrices up right away, which prevents later algebraic manipulations from spoiling the gauge invariance.Namely, we start with the following replacement where like before the wavy arrow denotes the replacement of a factor in the denominator by its respective reduction matrix.The last term M 5, 3  1 is present more as a formality and obviously the identity.The first two reduction matrices above take the following forms, M 5,3 Summing up these three reduction matrices, we obtain the quantity denoted as M 5,3 of the following form, where in the first equal sign the subscript (i) of the canonical-coefficient vector denotes the i-th term in (4.33) it corresponds to.In the last equal sign we have decomposed the operator associated with I 2 into the unit basis.The explicit expression for this operator is given by using the reduction matrix M 5,2 σ 2 −σ 3 once more to replace the factor (σ 2 − σ 3 ), leaving us with a monomial σ 2 in the denominator and thus the integrand of the prepared form.This reduction matrix reads, M All other terms in (4.14) can be dealt with in the same fashion.Since the two terms we have discussed in detail are representative enough, the remaining ones possess no conceptual novelty and we present the final expression of them in Appendix A. To verify this result, one can always cast the expressions into the DDM basis, as we have discussed in Section 2.2 and shown in the 4-point case, by collecting the coefficients of the bilinears, such as In doing so, it is easy to observe that those non-local poles disappear.

Conclusion and Outlooks
In this paper, we have developed a systematic method to construct the BCJ numerators, starting from the CHY forms of scattering amplitudes and using the differential operator proposed in our previous work [48].This method is based on the key observation that the number of canonical coefficients in such a differential operator is always (n − 3)! for a n-point CHY form, independent of the order of the operator.In the process of solving for the canonical coefficients, we have built the reduction matrices to simplify the differential operator and improve the efficiency of computation.The reduction matrices are also universal for all the theories.In the end, as we have demonstrated, we always arrive at the prepared forms for which the coefficients are solved analytically in [49].The values of these canonical a-coefficients are just the BCJ numerators in the unit basis we have constructed, and these canonical coefficients make the evaluation of the CHY forms a quick exercise.
Both the BCJ numerators and the amplitudes obtained this way enjoy the manifest gauge invariance in (n − 2) out of the n external legs, and their final expressions are always of factorized forms.
It is hopeful to formulate a closed form for the reduction matrix in future works, since a given reduction matrix M σ ij is determined only by the factor σ ij and the polynomial scattering equations.The reduction matrix behaves linearly as 1 σ ij in practice.Moreover, the polynomial scattering equations have nice combinatoric structures to exploit, which may allow us to produce analytical results for a general reduction matrix.As discussed in [51], the polynomial scattering equations is a Macaulay H-basis.This property may be helpful to prove the observation of the number (n − 3)! of canonical a-basis.
In this paper, the concept of dual basis is adopted to extract the non-local BCJ numerators in the minimal basis, in which the non-local propagators can be removed using the BCJ relations.It is another interesting direction to look for the algebra that extracts the local BCJ numerators directly, in like manners as the dual basis does the non-local ones.Such an algebra should also be theory-independent, just like the dual basis projectors.A possible candidate may be the Hopf algebra [71] or some improvised version of it.
Besides the algebraic properties, it is also important to investigate the geometric aspects of the BCJ numerators.The geometry of the scalar amplitudes, namely the cubic structures, is scrutinized in the language of associahedron in [69,72,73].There can be ways to dress these geometries with numerators encoding the kinematic information, such that the dressed geometries capture the properties of other theories, such as gauge theory and gravity [74][75][76].
Our method can be generalized to higher loops directly.Even if the generalized CHY forms for higher-loop amplitudes remain unknown, the dimension of the looplevel BCJ basis can still be calculated from the scattering equations.Another future direction we would like to take is to carry our method over to string theory.String amplitudes, written in the forms of worldsheet integrals, have a number of features in common with the CHY-like integrals.It is reasonable to hope that there exist similar differential operators and even reduction matrices which can help us evaluate those worldsheet integrals efficiently while preserving the factorized form.
A Five-Point (n − 2)-leg Gauge Invariant Form In Section 4.2 we have displayed the process of evaluating two terms in (4.14), using the reduction matrices.In this appendix, we present the evaluation of the remaining terms in (4.14).It turns out that when the terms containing the same F i 1 •••ir ǫ 1 ,ǫ 5 are grouped together, the poles which are not present in the cubic structures drop out.Thus in this way we show the final expressions of these terms.Moreover, in further simplifying the expressions, the following Jacobi-like identities come in handy, The terms containing the factor F ǫ 1 ,ǫ 5 read, where the transformation matrix B T 5 is already taken care of in the final expression.The terms containing the factor F i ǫ 1 ,ǫ 5 read,   (A.9)

B 6-Point Reduction Matrix
For higher-multiplicity amplitudes, just like their 4-and 5-point counterparts, only the reduction matrices associated with σ ij and the canonical coefficients of prepared forms are the necessary ingredients in the construction of the BCJ numerators and the evaluation of a given CHY form.In this appendix, we present the explicit expressions for the reduction matrices and the canonical coefficients we may encounter in the case of 6 external lines.These quantities are theory-independent.The prepared forms we may end up with are always of the form A 6 σ i , where i = 0, 2, 3, 4. It is easy to observe that the degree of h 0 appearing in the 6-point CHY form can be at most 4, resulting in the possible reduction matrices being the following, M 6,4 σ i −σ j , M 6,3 σ i −σ j , M 6,2 σ i −σ j , M 6,2 σ i .
The three independent polynomial scattering equations read The canonical coefficients forming the vector A σ i are obtained by solving the local duality constraints and the intersection number constraint for the ideal {h 1 , h 2 , h 3 , h 0 }, where h 0 = σ i .Such cases are studied in detail in [49] and the analytical solutions can be directly read off from the tableaux representations as the following, As for the reduction matrices, those associated with the factor σ 0 can be immediately read off from the definition of the reduction matrix, namely, M 6,2 σ 0 = I , M 6,3 σ 0 = 1 2 I , M 6,4 σ 0 = 3 × 2) The factors multiplying the identity matrix all come from taking the derivatives.
The other reduction matrices are slightly more complicated.We take M 6,2 σ 2 as an example and the rest can all be obtained the same way.The equations relating the the two sets of the canonical coefficients associated with M  where f = 5 i=1 (−1) δ i,5 s 1i s i6 and L 6 σ 2 is The reduction matrices of the form M 6,2 σ i −σ j can obtained in the same fashion and those of the forms M 6,3 σ i −σ j and M 6,4 σ i −σ j can be related to the former, that is, M 6,2 σ i −σ j = (L 6 σ i −σ j R 6 ) −1 , M 6,3 σ i −σ j = 1 2 M 6,2 σ i −σ j , M 6,4 σ i −σ j = 3

2
Ambiguity of BCJ Numerator 3 BCJ Basis from Differential Operator 3.1 Canonical Coefficients in Differential Operator 3.2 Discussion on the Dimension of BCJ Basis 3.3 Reduction Matrix 4 Gauge Invariance of Yang-Mills BCJ Numerators 4.1 Yang-Mills 4-Point Amplitude 4.2 Yang-Mills 5-Point Amplitude 5 Conclusion and Outlooks A Five-Point (n − 2)-leg Gauge Invariant Form B 6-Point Reduction Matrix 1 Introduction

Figure 1 :
Figure 1: Four numbers of the value (n − 3)! at tree level deg(P ) must be the same on both sides, giving rise to a set of linear relations between the two sets of canonical a-coefficients.Schematically, solving these relations we obtain, γ(1)•••γ(n−4)(d−1)

. 25 )
Comparing (3.24) and (3.25), we arrive at of residue computing.Their respective coefficients are related by a transformation matrix whose elements are only integers arising from simple combinatoric factors and differentiations.Note that those coefficients as 2 •••s n−2 s 0 with s 2 + • • • s n−2 > m 0 , which are not related to any b s2 •••s n−2 s0, do not contribute to the residue anyway.The non-canonical a-coefficients can always be rewritten in terms of the canonical ones that can be further expressed in terms of the reduction matrices and the a-coefficients belonging to the lower-order operators.That is, D(m 0 ) Γ can indeed be related through the relevant reduction matrices to D (m 0 ) Γ ′ , which is a standard operator of the form (3.5) associated with Γ ′ = q(σ)Γ where deg(q) = m − m 0 .Moreover, with all the coefficients in D(m 0 ) Γ expressed in terms of the canonical a-coefficients in D
Notice that the coefficients appearing in the these relations are not canonical and these non-canonical coefficients can be written in terms of the canonical ones, which is realized is a 6 × 10 matrix while R 6 a 10 × 6 one.The matrix R 6 is determined solely by the scattering equations and therefore universal for a given number of external legs, which reads− 2s 12 s 46 [i,j s k,h],l ≡ s ij s khl − s hj s kil .The reduction matrix is then M 6,2 σ 2 = (L 6 σ 2 R 6 ) −1 . (B.8) [2,6 s 3,5],6 s 36 s 326 0 s 26 s 256 36 s 346 + 1 − s 46 s 246 s 36 s 346 − s 46 s 256 2s 36 s 346 s s