Expansion of EYM Amplitudes in Gauge Invariant Vector Space

Motivated by the problem of expanding single-trace tree-level amplitude of Einstein-Yang-Mills theory to the BCJ basis of Yang-Mills amplitudes, we present an alternative expansion formula in the gauge invariant vector space. Starting from a generic vector space consisting of polynomials of momenta and polarization vectors, we define a new sub-space as gauge invariant vector space by imposing constraints of gauge invariant conditions. To characterize this sub-space, we compute its dimension and construct an explicit gauge invariant basis from it. We propose an expansion formula in the gauge invariant basis with expansion coefficients being linear combinations of Yang-Mills amplitude, manifesting the gauge invariance of both expansion basis and coefficients. With help of quivers, we compute the expansion coefficients via differential operators and demonstrate the general expansion algorithm by several examples.

1 Introduction recursive formula, at the end the EYM amplitude would be expanded to the basis of Yang-Mills amplitudes with legs 1 and n being fixed. Coefficients of each Yang-Mills amplitude is a linear combination of C ha (h), which are polynomial functions of polarization vectors and momenta whose precise definition can be found in literatures [18].
While the expansion of EYM amplitude in KK basis of Yang-Mills amplitudes has been solved completely, since KK basis is not the minimal basis of color-ordered Yang-Mills amplitudes, a question naturally arises: what would happen when expanding an EYM amplitude to the minimal basis, i.e., the BCJ basis of Yang-Mills amplitudes? In a first thought, it seems that this question has already been solved by the generalized KLT relation (1.1). However in (1.1) the momentum kernel S[σ|σ] and A R are difficult to compute and we also need to sum over all S n−3 permutations. Hence the generalized KLT relation dose not work well in practical computation. One could also start with expression (1.2) and reformulate KK basis to BCJ basis by BCJ relations. However, computation of several examples is suffice to suggest that the algebraic manipulations are rather complicated. The resulting expansion coefficients are rather cumbersome without any hints of systematic and compact reorganization, because there are too many equivalent expressions. In paper [21], a new method is proposed by introducing the differential operators into this problem. The differential operator is originally applied to the research of relating amplitudes of different theories [22], and later on a series work shows how to apply differential operators to the expansion of EYM amplitude to KK basis [21,23,24]. Furthermore, attempts of generalizing differential operators in the expansion of EYM amplitude into BCJ basis have been carried out for some simple cases where EYM amplitudes contain one, two or three gravitons. However a systematic method for generic EYM amplitude with n gluons and m gravitons is still in demand.
In this paper, we are trying to fulfill this request by providing a systematic method to compute the expansion coefficients of EYM amplitude with m graviton in the BCJ basis. Besides the use of differential operators, we would also introduce the principle of gauge invariance. Since Yang-Mills amplitudes of BCJ basis are linearly independent, if we can write an EYM amplitude as Yang-Mills amplitudes in BCJ basis, the gauge invariance of polarization tensors of gravitons would be transformed partially into the gauge invariance of expansion coefficients, encoding in the polarization vectors. Hence the gauge invariance sets strong constraints on the form of the expansion coefficients. In fact, the gauge invariance principle has already played important role in the study of scattering amplitude. It is expected that the gauge invariance could completely determine the amplitudes of certain field theories [25][26][27], and further exploration can be found in various aspects [18,22,28,29]. Especially demonstrated in [18], it is the constraints of gauge invariance that make a compact formula available for expansion of EYM amplitude in KK basis. However the potential application of gauge invariance is still not fully exploited. In this paper, we would like to take a different view on the understanding of gauge invariance. Just as what have been done for the symmetries in amplitudes of N = 4 super-Yang-Mills theory, since the principle of gauge invariance is a strong constraint for field theory, we would like to make it manifest in the level of scattering amplitudes.
With the new understanding of gauge invariance, in this paper we will show how to expand general EYM amplitude into BCJ basis of Yang-Mills amplitudes systematically. Organization for this paper is as follows. In §2, we review some backgrounds. In §3, we introduce the gauge invariant vector space living in a general vector space consisting polynomials of Lorentz contractions of momenta and polarization vectors. We compute the dimension of gauge invariant space, characterize the explicit form of vectors, and finally construct the gauge invariant basis. In §4, we define gauge invariant vectors and differential operators in quiver representations, which is the description of mathematical structures of these vectors and operators. With help of quivers, we implement a systematic algorithm to compute expansion coefficients. In §5, we illustrate our method by several explicit examples, the EYM amplitudes with up to four gravitons in the purpose of clarifying some subtleties. In §6, we conclude our discussion and point out some problems to be solved in future. Detailed proofs of some propositions as well as some explicit BCJ coefficients in BCJ relations are presented in appendices.

The expansion of EYM amplitudes to Yang-Mills amplitudes in BCJ basis
In this section, we review some background knowledge which is useful in the later discussion of expanding EYM amplitude to BCJ basis of Yang-Mills amplitudes. Firstly, as reviewed in [21], an arbitrary colorordered Yang-Mills amplitude can be expanded to BCJ basis with three particles being fixed in certain positions relating to the color-ordering, as A n (1, β 1 , ..., β r , 2, α 1 , ..., α n−r−3 , n) = (2.1) The expansion coefficients, namely BCJ coefficients, are firstly conjectured in [3] and later proven in [4], with the expression Secondly, we review the differential operators which are originally introduced in [22]. An important differential operator is the insertion operator defined by Physically it stands for inserting a graviton k between i and j when i, j are adjacent gluons in a trace. If i, j are not adjacent, for instance T ik(i+2) , we can write it as and its physical meaning is also clear 2 . Another important operator is the gauge invariant differential operator defined as . (2.5) It has a physical meaning of imposing gauge conditions, i.e., changing a → k a . For an arbitrary polynomial of polarization vectors and momenta, if it vanishes under operator G a , we can conclude it is gauge invariant for gluon a. Gauge invariant operators are commutative, i.e., [G a , G b ] = 0, so the result of a sequential operators does not depend on the ordering, and we can denote a sequential gauge invariant operator as The insertion operator and gauge invariant operator satisfy the following commutative relation, with T ij := ∂ ( i · j ) , and it is valid after applying to any functions of polarization vectors and momenta 3 . Finally let us present a general discussion on the expansion of EYM amplitude to BCJ basis. For particles with spin, the corresponding Lorentz representations are carried out by polarizations, e.g., polarization vector˜ µ i for gluon and polarization tensor µν h i for graviton. When expanding EYM amplitude to BCJ basis, the polarization tensor of graviton is factorized into two parts µν is inherent by the polarization vector of gluon in Yang-Mills basis, while the other part ν h i is absorbed into expansion coefficients. More explicitly, the expansion coefficients are rational function of momenta k µ κ , κ = 1, . . . , n, h 1 , . . . , h m and polarization vectors µ hκ , κ = 1, . . . , m. A crucial difference between expanding to KK basis and BCJ basis is that, the BCJ basis is truly an algebraic independent basis and the corresponding expansion coefficients must be gauge invariant, i.e., In the former formulation (2.8), independent Yang-Mills amplitudes are taken to be expansion basis, and each coefficient as a function of momenta and polarization vectors hκ should satisfy conditions of gauge invariance for all hκ with κ = h 1 , h 2 , . . . , h m . In the latter formulation (2.9), b gauge-inv 's are the expansion basis and the expansion coefficients become a linear summation of A YM 's times rational functions of momenta. This alternative form of expansion has appeared in paper [21], where in order to distinguish two different kinds of basis we named the later one b gauge-inv as building blocks 4 .

Building up expansion basis in gauge invariant vector space
As mentioned, in the expansion of EYM amplitude, the gauge invariant coefficients c gauge-inv as well as expansion basis b gauge-inv are crucial. They are polynomials of polarizations and momenta which vanish under gauge conditions. In this section we would like to start from a most general vector space and localize a gauge invariant sub-space from it. The expansion basis we are looking for is living in this sub-space.

Gauge invariant vector space and its dimension
Let us start from the most general polynomial h, constructed by Lorentz products of n momenta k 1 , k 2 , . . . , k n and m polarizations 1 , . . . , m with m ≤ n. By Lorentz invariance and multi-linearity of i , this polynomial must be the form where for each monomial the degree of is m and each i , i = 1, . . . , m appears once and only once, while the coefficients α's are rational functions of Mandelstam variables of momenta. If we take all monomials as vector basis 5 , we can build up a vector space V n,{ 1 ,..., m} over the filed of rational functions of Mandelstam variables, where any polynomial h n,m belongs to this vector space.
In order to find out the gauge invariant vector space from V n,{ 1 ,..., m} , let us impose gauge invariant condition on h n,m . This can be achieved by applying differential operators G i 's to (3.1), i.e., Such operator establishes a linear mapping between different vector spaces as where in the resulting vector space the polarization t has been replaced by k t and no longer appear, which is denoted by t . This linear map is surjective 6 by noticing the reduction of B[V ], i.e., We can successively apply different gauge invariant operators G i 's, i = 1, . . . , m and establish a mapping chain of vector spaces. Since all G i 's are commutative, the result dose not depend on the ordering of successive applying, and we can denote the mapping chain as The superscript labels the removed polarizations i in the vector space. Note that different ordering of applying G i 's produces different mapping chains which at the end leads to the same vector space, so (3.4) in fact represents a collection of mapping chains. The kernel of linear map These monomials are not linearly independent. There are relations by momentum conservation i ki = 0 and transverse condition i · ki = 0. Furthermore, we consider only the parity even case, i.e., without total antisymmetric tensor µ 1 ...µ D . 6 The property of surjectivity is the cornerstone in our discussion, however for vector space of polynomials without ( · k) m it no longer holds.
Physically it means that the vectors of kernel are gauge invariant for i-th particle. Using the fact that the linear map is surjective (3.4), by fundamental theorem of linear map, we get Then the dimension of kernel can be computed by the difference of dimensions of vector space as When applying more than one G i 's, this relation can be generalized to For example let us consider the simplest case s = 0, Vector space V n,0 is the field of rational functions of Mandelstam variables, so the basis is just 1 and dim V n,0 = 1. For vector space with only one polarization, the kernel Ker G 1 [V n,1 ] represents all vectors vanishing under gauge invariant operator. This is the gauge invariant vector sub-space W n,1 in a vector space V n,1 . Thus we get For a general vector space V n,m with m polarizations, we can define the gauge invariant vector sub-space as the intersection of kernels of all possible linear maps G i 's as, This means that a vector in W n,m would vanish for any linear map G i . This is exactly the sub-space where lives all gauge invariant coefficients c gauge-inv of (2.8) and the expansion basis b gauge-inv of (2.9). Let us compute the dimension of W n,m , for example when m = 2. Generally, for any two linear spaces U 1 , U 2 , we have the following relation for the dimension, Apply this relation to the vector space of kernels, i.e., The first two terms in RHS can be computed by (3.7), while in order to compute the third term, we need to use the following proposition 7 , PROPOSITION 1 : any two kernels of linear maps G i 's satisfy the splitting formula, 14) and its generalization, PROPOSITION 1 EXTENDED : the kernels of linear maps G i 's satisfy the generalized splitting formula, Together with (3.8), we can rewrite (3.13) as Recursively using (3.12), we are able to generalize above result to arbitrary m. For simplicity let us denote U i := Ker G i , and when m = 3 we get In the second line, the first three terms have already been computed, while in order to compute the fourth term we need to use the following proposition 8 PROPOSITION 2 : three kernels of linear maps G i 's satisfy the distribution formula, and its generalization PPROPOSITION 2 EXTENDED : the kernels of linear maps G i 's satisfy the generalized distribution formula, Together with (3.12), we can rewrite (3.17) as is not true. For example, in a two-dimension space U , let us choose U1, U2, U3 to be line y = 0, x = 0 and x = y respectively. Then U1 + U2 is the whole XY-plane, and (U1 + U2) ∩ U3 is the line x = y. While in the RHS, U1 ∩ U3 and U2 ∩ U3 are just the origin (0, 0). So the RHS is a point.
In equation (3.20), in order to compute the dimension dim W n,3 := dim(Ker G 1 ∩ Ker G 2 ∩ Ker G 3 ), we need the result of dim(Ker G 1 + Ker G 2 + Ker G 3 ), which by proposition 1 extended (3.15) it equals to dim Ker G 123 . Using (3.8), we get Notice that the numerical factors 1, 3, 3, 1 are nothing but 3 i for i = 0, 1, 2, 3. Let us proceed further to arbitrary m. With proposition 1 extended and proposition 2 extended, equations (3.13) and (3.20) are exactly the same as the principle of inclusion-exclusion. By the well-known principle of inclusion-exclusion, we get where the second summation is over all subsets with s indices. It is also well-known that starting from the principle of inclusion-exclusion we can arrive at By proposition 1 extended, we can write where the dimension of vector space V n,m can be computed via 9 The counting of (3.27) can be carried out as follows. Firstly we select i pairs of , and there are m 2i choices, while each left can be contracted with (n − 2) momenta after ( · kn) by momentum conservation. For 2i 's, the number of different contractions is (2i)! 2 i (i!) .
Hence the dimension of arbitrary gauge invariant vector space W n,m can be computed by formula (3.26) and (3.27 In paper [25] the same result has been provided up to n = 7. 10 . Comparing with that result, our calculation shows more efficiency than that of solving linear equations of gauge invariance directly. Furthermore, several examples of dim W n,m and dim W n+m,m with arbitrary n but definite value of m are listed below as dim W n+m,m n − 2 (n − 1) 2 + 1 n 3 + 3n (n + 1) 4 + 6(n + 1) 2 + 3

Gauge invariant vectors
The dimension of gauge invariant vector space characterizes the minimal number of vectors to expand an arbitrary vector, while the explicit form of vector is not constrained. For the working experiences of EYM amplitude expansion with one, two and three gravitons [21], we get the insight that the coefficients appearing therein could be recast in a manifestly gauge invariant form as linear combinations of multiplications of fundamental f -terms.
Here the fundamental f -terms stand for two types of Lorentz contractions of field strength f µν i = k µ i ν i − µ i k ν i and external momenta, with at most two f i 's, Fundamental f -terms: This observation can be generalized beyond m = 3, and it can be stated as follows. For any vectors in gauge invariant vector space W n,m with m < n 11 , Every vector in W n,m can be rewritten in a manifest gauge invariant form, which is a linear combination of the multiplication of fundamental f -terms with the rank of f i being m.
We shall prove this statement by induction. The cases with m = 1, 2, 3 have already been shown to be true in [21]. Following the idea of induction, we assume that this statement is true for all s < m, and prove that it must be true for m. 10 In paper [25], there are two types of spaces. The other one is the space with at least one contraction between polarization vectors in polynomials, i.e., polynomials without monomial ( ·k) m , which is exactly the vector space that Yang-Mills amplitudes live in. Its dimension is (n − 3)!. 11 We should emphasize the condition m < n, which is different from previous discussion where m could equal to n. Proof of the statement in this subsection can not be trivially generalized to the m = n case, so if results in this subsection could be applied to the case m = n is still a question for us.
A polynomial h n,m ∈ W n,m with m polarizations 1 , 2 , . . . , m can be generally written as where momentum conservation has been applied to eliminate 1 · k n , so that all ( 1 · i ), ( 1 · k i ) appearing in h n,m are linearly independent. Polynomials T 1i ∈ V n,m−2 and i · T 1i , T 1i ∈ V n,m−1 . Since h n,m ∈ W n,m , by definition we have G a h n,m = 0 , ∀(1 ≤ a ≤ m) .
where we have considered the fact that h n,m does not contain ( 1 · k n ). With above result we can rewrite h n,m as We also need to consider the gauge invariance of h n,m with respect to polarization vector 1 , Then we get After substituting above results back to h n,m , we get So h n,m is already manifestly gauge invariant for polarization vector 1 . In fact, we can also choose to eliminate other coefficients in (3.34) and introduce different poles in denominator of h n,m . We can also generate another set of equations by considering the operator relations [T i1n , G a ] = 0 with i = m + 1, · · · , n − 2 and a = 2, · · · , m. Applying them to h n,m produces which means T 1i 1 is gauge invariant for 2 , 3 , · · · , m . By assumption of induction, T 1i 1 can be written as a linear combination of multiplication of fundamental f -terms. Because h n,m and T 1i 1 are gauge invariant for a with a = 2, · · · , m, and (k n−1 f 1 f i 1 T 1i 1 )'s are linearly independent, T 1i 1 is also gauge invariant for all its own polarization vectors. Again by assumption of induction, any (Af i 1 T 1i 1 ) can also be written in a manifest gauge invariant form with only f appears. Thus as a linear function of (k n−1 · f 1 · f i 1 · T 1i 1 ) and T 1i 1 , the polynomial h n,m can also be written in a manifest gauge invariant form, and we have proven the first part of our statement.
To complete our proof, we need to apply above procedure to ( i 1 · T 1i 1 ) in (3.29) and rewrite it as where in the last summation i 2 can equal to 1. Let us again apply operator equations [T ai 1 n , G a ] = T ai 1 , with a = 2, · · · , m and a = i 1 , which generates a set of equations, (3.39) Then we apply [T ji 1 n , G a ] = 0 with j = m + 1, m + 2, · · · , n − 1, 1 and a = 2, · · · , i 1 − 1, i 1 + 1, · · · , m to ( i 1 · T 1i 1 ), which leads to G a T 1i 1 j = 0. It says that T 1i 1 i 2 is gauge invariant for its own polarization vectors, and it can be written as linear combination of multiplication of fundamental f -terms. For the same reason as before, we conclude that T 1i 1 i 2 is also gauge invariant for its own polarization vectors. Continuously applying the same procedure to T until to the last polarization vector, we would arrive at To further reduce the expression (k · f · · · f · k) to the fundamental f -terms, we should get help from the following identities, where A, B, C could be any strings. More explicitly, applying above identity to expression with three f 's, we get So any f -term with any number of f i 's can be reduced to fundamental f -terms, while at the same time T (1i 1 ···i s−1 )is has been reorganized as a linear combination of multiplication of fundamental f -terms. This ends the proof of statement by induction method. Before ending this subsection, let us take a look on another gauge invariant f -term that mentioned in [21], i.e., the trace Tr(f a 1 f a 2 · · · f a k ) = f µν a 1 f a 2 ,νρ · · · f σ a k ,µ . It can be expanded as where identity (3.42) has been used in the derivation. Combining the first and third term as well as the second and fourth term, we can get A simple example is Tr(f a 1 f a 2 ) = 2(k a 2 f a 1 f a 2 A)/(Ak a 2 ). So this type of gauge invariant f -terms, which is originally viewed as a new type different from (kf · · · f k), are also composed by fundamental f -term.

Gauge invariant basis
Any gauge invariant vector in W n,m could be an element to form a gauge invariant basis b gauge-inv in the EYM amplitude expansion (2.9). However, in order to turn a subset of W n,m to a complete basis, we should choose a set of vectors satisfying the following two properties, 1. all vectors in the set are linearly independent, 2. the number of vectors in the set equals to the dimension of gauge invariant vector space.
Note that the fundamental f -terms are not completely independent from each other. For instance, using (3.42) it is easy to see that So one can always reduce any fundamental f -terms to the following form, From the definition of f µν i , it's easy to get In the case of A EYM n,m , the momentum list is {k 1 , . . . , k n , k h 1 , . . . , k hm } while the polarization vector list is { h 1 , . . . , hm }, so by default the above subscripts a, b ∈ {h 1 , . . . , h m } and i ∈ {1, . . . , n, h 1 , . . . , h m }. After using momentum conservation to eliminate k n , we can restrict the fundamental f -terms to be Using above fundamental f -terms, we can construct a set of vectors as with the convention The linear independence of these vectors (3.50) is obvious. In order to demonstrate that they form an independent basis of W n+m,m , we should show the total number of these vectors equals to dim W n+m,m according to property 2. We can count the total number of independent vectors with respect to specific s as According to (3.26) and (3.27), the dimension of W n+m,m is Noticing the relation we immediately get dim W n+m,m = #(vectors) defined in (3.50). Hence the set of vectors defined in (3.50) satisfies the required two conditions and could be chosen as an expansion basis for A EYM n,m in (2.9). In practice we would prefer a basis with minimal dimension, then we define the fundamental f -terms as where K a := a i=2 k i . The vectors in the expansion basis can be constructed from above fundamental f -terms as with the convention They contribute to a complete set of expansion basis, and a general EYM amplitude can be expanded into this basis as where H/h i is the set of gravitons excluding h i , and the three sets . . , h γr } with 2p + q + r = m are a splitting of all gravitons. The reduced summation runs over all possible splittings a ∪ b ∪ c = H and the prime means that terms with index circle should be excluded 12 . We can see that all the information of polarization vectors hκ is encoded in B as expected.

Determining expansion coefficients via differential operators
We have defined the gauge invariant expansion basis, and the next step is to determine the expansion coefficients. As earlier mentioned, the EYM amplitude can be expanded schematically in the form, or more explicitly see (3.59). The expansion coefficients are linear combinations of Yang-Mills amplitudes A YM n+m . To use (4.1) efficiently, a crucial point is to find a way to distinguish vectors in the gauge invariant basis from each other. Inferred from the explicit form of vectors in (3.1), we notice that the signature of vectors is the structure ( · ) p ( · k) q , where each ( · ) and ( · k) is linearly independent. This motives us to consider two kinds of differential operators as .
Applying these operators on RHS of (4.1) will vanish all terms but those containing corresponding ( · ) and ( ·k). While applying these operators on LHS of (4.1), the physical meaning will be different. Applying T ab on single-trace EYM amplitudes produces multi-trace EYM amplitudes which would sightly complicate the amplitude expansion, however applying T ah i b on single-trace EYM amplitude produces another single-trace EYM amplitude but with one less graviton. It would transform the graviton h i to a gluon h i and insert the gluon in between positions of gluons a, b respecting to the order-ordering. So each time applying an insertion operator on (4.1), the number of gravitons is reduced by one, while a multiplication of m insertion operators would transform the LHS of (4.1) to Yang-Mills amplitudes completely, as expected 13 .
In fact, we can go a step further and define a differential operator as multiplication of m properly chosen insertion operators. When applying the differential operator on (3.59), there would be one and only one new vector in RHS of (3.59) survive, As a consequence, we get one linear equation with only one unknown variable, and the corresponding expansion coefficient can be computed directly as a function of A YM n+m that generated by differential operator applying on RHS of (3.59) 14 . The problem of EYM amplitude expansion is then translated to the construction of properly defined differential operators, which would be the main purpose of this section. We find it very helpful to use quivers to represent the gauge invariant basis and differential operators for our purpose.

The gauge invariant basis and its quiver representation
The definition of insertion operator (4.2) indicates that a differential operator would only affect the Lorentz contraction ( h · k), so all other types of Lorentz contractions (k · k) and ( · ) can be treated as unrelated factors. In order to characterize the structure of ( h · k) in a gauge invariant vector, we can assign a quiver, i.e., directed graph, to it 15 . We call a directed graph representing all ( · k)'s of a vector as ( k)-quiver 13 Alternatively, we could also apply less insertion operators to generate a set of linear equations, and recursively use the expansion of single-trace EYM amplitude with less number of gravitons into Yang-Mills amplitudes.
14 Similar idea of selecting only one unknown variable at each step has already been considerably applied in the OPP reduction method [30] for one-loop amplitude. 15 The idea of using arrows to represent Lorentz contractions has already been applied in literatures [31,32], where all types of Lorentz contractions are considered. However, we are only interested in Lorentz contraction of the type · k in this paper.
of the vector. In a quiver, we use a directed solid line to represent ( h i · k h j ) with an arrow pointing to a graviton momentum k h i , and a directed dashed line to represent ( h i · k j ) with an arrow pointing to a gluon momentum k j as 16 As for the fundamental f -term (3.48), which can be expanded as then we can assign three ( k)-directed graphs to it as 17 Since each graph denotes a multiplication of ( · k) terms, hence applying the following derivatives , we will get non-vanishing results. Similarly, for (k 1 · f h i · k) we can assign two ( k)-directed graphs as 16 From now on, we will identify an directed line with its corresponding ( · k) term, and sometimes when we refer to a specific directed line connecting two nodes from a to b, we will use the label (ab). 17 Notice that there are four terms in the expansion of k1 · f h i · f h j · k1, while the h i · h j term is the most crucial signature to distinguish it from other fundamental f -terms. However in this paper we only consider insertion operators so that · is out of our sight. Figure 1. The quiver representation of fundamental f -terms.
Note that the factor ( h i ·k 1 ) exists in both (k 1 ·f h i ·f h j ·k 1 ) and (k 1 ·f h i ·k), so the action of derivative ∂ h i ·k 1 on them both are non-zero. Consequently, we prefer to eliminate the dashed lines representing h i ·k 1 in the graphs of ( k)-quiver to obtain a simple presentation. Furthermore, to represent a fundamental f -term by one graph and distinguish .55) and (3.56) are represented by quivers in Fig.1. To distinguish these quivers from ( k)-quiver, we will call them basis quivers or just quivers. We should emphasize that from a basis quiver it's easy to recover all corresponding ( k)-quivers by replacing any one solid or dashed arrow ( h i · k) in the graph by a dashed arrow ( h i · k 1 ), i.e., from Fig.1 to (4.5),(4.7). However, given a ( k)-quiver, it is hard to tell which basis quiver it comes from, especially when there are many ( · k 1 ) lines. The fact that there is no one-to-one correspondence between basis quivers and ( k)-quivers causes some technical difficulties in the construction of differential operators. Fortunately, for a gauge invariant vector, its basis quiver and ( k)-quivers do possess a common property: they all contains m and only m lines (counting both dashed line and solid line), since each line carries one h i .
Note that the basis quiver for F h i h j is a colored loop, where colors are to remind us that it is an overlapping of three ( k)-quivers after eliminating dashed lines. We call such a colored loop as pseudo-loop. In general, there can also appear real loops. For example, the F h 2 However as explained in [21], the terms with indices or part of indices forming a closed circle will not present in the expansion of EYM amplitude, although such terms do appear in the gauge invariant basis. So we will exclude basis quivers with real loops in practical computation. Next let us come to consideration of the quiver representation of gauge invariant basis. As shown in (3.57), a vector in gauge invariant basis is a multiplication of fundamental f -terms, with indices following the convention (3.58). Since each h i appears only once in a vector, we can conclude that each point labelled by h i in the basis quiver of a gauge invariant vector has at most one out-going line · · · · · · · · · · · · single pesudo-loop one tree with one pesudo-loop two trees with one pesudo-loop single dashed line tree with one dashed line but possibly several in-coming lines. Consequently, all pseudo-loops are topological disconnected from each other. The point labelled by K a is connected by only in-coming lines but not out-going lines, hence all such points are also topological disconnected from each other. Furthermore, pseudo-loops can not be connected with points labelled by K a either. So a quiver graph could have many disconnected components, whose number is at least p and at most p + r. A line for F h j h i can be connected to one and only one disconnect component.
With above analysis, let us discuss the possible structures appearing in a quiver representation for gauge invariant basis (4.8). Firstly, since each F aγ i hγ i is represented by a dashed line with arrow pointing to K aγ i , a dashed line can never be followed by a pseudo-loop or solid line. Secondly, each F . . , h γr } it will be followed by a dashed line, while if h β i ∈ {h α 1 , . . . , h α 2p } it will be followed by a pseudo-loop. However if h β i ∈ {h β 1 , . . . , h βr }/{h β i }, for instance h β i = h β j it will be followed by another solid line, and the latter could further be followed by pseudo-loop, dashed line or solid line. The connection of solid lines should at the end stops at a dashed line or pseudo loop, otherwise it would form a real circle which should be excluded.
To summarize, the quiver representation of a vector in gauge invariant basis could contain the following sub-structures, These two examples illustrate our previous discussions very well. There are three disconnected components for the first one, and two for the second one. In the second graph, two dashed lines is connected to one node representing the fundamental f -terms F 4 h 3 , F 4 h 6 . All directed solid lines stop at pseudo-loops or dashed lines.
In fact, we can give a more precise description of the structures of basis quivers by using the concept of rooted tree [33]. The quiver of a vector in gauge invariant basis consists of some disconnected components and each component contains only one pseudo-loop or node K a i . If we focus on a disconnected component with node K a i , this is exactly a rooted tree with the root being node K a i . More precisely, it is a directed rooted tree with an orientation towards the root, i.e., the direction of all lines in the tree directs to the root from leaves, as illustrated in the previous two examples. For the disconnected component with a pseudoloop, we should split pseudo-loop into two colored lines and produce two graphs. For each graph, we take the node with only in-coming lines as the root, thus we obtain two rooted trees from a disconnected component with a pseudo-loop. The picture of rooted trees will help us to construct the differential operators and understand many properties of our algorithm later.

Constructing differential operators
In quiver representation, a vector will be non-vanishing under a derivative ∂ h i ·k only if its ( k)-quiver representation contains a solid or dashed line corresponding to h i · k, where k can either be a graviton momentum or a gluon momentum. A vector in gauge invariant basis is a multiplication of some ( · k)'s, hence by constructing a differential operator as a proper multiplication of some derivatives ∂ h i ·k 's, we expect all but one vector being vanishing, so it can select a particular non-vanishing vector in gauge invariant basis. The differential operator can be constructed by three types of insertion operators (4.2), which we will introduce next. The first type of insertion operator takes the form, where k a is the momentum of a gluon. A vector would vanish under T ah i (a+1) unless its ( k)-quiver representation contains a dashed line corresponding to h i · k a or h i · k a+1 . Applying this insertion operator to the three types of fundamental f -terms we immediately get, and The above results tell us that if the basis quiver of a vector in gauge invariant basis contains a dashed line representing F a hγ i , then there exists a differential operator containing an insertion operator T ah i (a+1) that will select this vector and other vectors also having the same dashed line. The relation (4.12) can be graphically represented as, The second type of insertion operators takes the form where the Lorentz contraction of a polarization vector with a graviton momentum has been included. Since by definition the momentum k n does not appear in fundamental f -terms, when applying T h j h i n to them only the derivative ∂ h i ·k h j works. Explicitly, we get represented in quivers as At first sight, both F h i h j and F h j h i are non-vanishing under T h j h i n , so we may conclude this insertion operator is not sufficient to distinguish the two terms. However, we should note that the insertion operator is actually a selecting operator used to distinguish a vector with the contraction ( h i · k j ) from others, rather than distinguishing F h i h j and F h j h i . According to this point of view, it is immediate to accept that, since the quivers of F h i h j and F h j h i both contain the solid line from h i to h j , the insertion operator F h i h j would have non-vanishing results acting on them.
In order to construct a differential operator that can distinguish F h i h j from F h j h i , we need to consider another type of insertion operator. Since in the corresponding ( k)-quiver of F h i h j there are always two lines linked together, a solid line (h i h j ) and a dashed line (k 1 k h i ) or (k 1 k h j ), we can multiply T h j h i n by an additional insertion operator containing derivative ∂ h j ·k 1 . For this purpose, we can construct the operator T 1h j 2 T h j h i n , and applying it to F h i h j we have It is non-vanishing only when i = i, j = j, hence it distinguishes the pseudo-loop of F h i h j from all other pseudo-loops. However T 1h j 2 causes some additional troubles, for it will give non-vanishing results when actting other fundamental f -terms as because the two f -terms both contain a contraction ( · k 1 ). Combining these results, we will have 18 It means that although T 1h j 2 T h j h i n is able to distinguish one pseudo-loop from the others, it would mix contributions from vectors without pseudo-loop. However, it is not a problem at all, if we try to solve the coefficients of basis in a proper algorithmic way. We can firstly compute the coefficients of F With above discussions, let us present a general picture of associating a differential operator for a vector in the gauge invariant basis. The first step is to construct the corresponding ( k)-quivers from the vector's basis quiver. Because of the trouble caused by ( h i · k 1 ) in (4.16), we should avoid these graphs with dashed line( h i · k 1 ). So we only consider these graphs got by the following method, the dashed line in basis quiver is mapped to ( h i · K a ) while the solid line will be mapped to ( h i · k h j ), and each pseudo-loop has two possible ways of mapping, and can be mapped either to a solid line ( h i · k h j ) connecting a dashed line( h j · k 1 ) or to a solid line ( h j · k h i ) connecting a dashed line ( h i · k 1 ). We can randomly choose one to represent a pseudo-loop. After above mapping, we generate a corresponding ( k)-quiver. Now come to the crucial points: the ( k)-quiver is a collection of rooted trees. Pseudo-loops have already been mapped to lines, and they become different branches attaching to the root k 1 through dashed lines. As a comparison, the other rooted trees which comes from components without pseudo-loops have roots K a instead of k 1 . The collection of rooted trees for a vector in gauge invariant basis can be algebraically represented as the embedded structure where at each level we write as {root : leaf 1; ...; leaf m}. For example, the second quiver in (4.9) can be represented as With the view of rooted trees, we can construct corresponding differential operators as • (a) Assign operator T ah i (a+1) to each dashed line (h i K a ), which uniquely picks up the corresponding dashed line in the basis quiver.
• (b) Assign operator T h j h i n to each solid line (h i h j ), which uniquely picks up the corresponding solid line in the basis quiver.
• (c) Assign operator (k 1 · k h i )T 1h i 2 to each dashed line (h i k 1 ).
Above rules can be represented pictorially as (4.19) The corresponding differential operator for a vector in gauge invariant basis is defined by multiplying all assigned operators together. Since there is one-to-one mapping between ( k)-quivers and differential operator, we call the ( k)-quivers constructed according to above rule as D-quivers. We want to emphasize that, a D-quiver in fact represents two things at the same time: (1) it defines an unique differential operator, (2) it is a special choice of ( k)-quiver, which can be associated to a given basis quiver.
Above discussion can be summarized as the following expression stating the mapping from a given vector to a differential operator, with the subscripts following convention (3.58). There are several technical points we want to explain. Firstly, the mapping rule is defined such that Secondly, although insertion operators are commutative, when applying on EYM amplitudes we need to follow proper ordering to make the physical picture clear. We shall apply insertion operators of the type T ahγ a , T 1hα2 first, then the types T hαh α n and T h β h β n . More explicitly, the ordering of applying insertion operators is from the roots to the leaves in the D-quiver opposite to the direction of arrows. In fact, we can make the result more concrete when acting D i on A EYM n,m . As mentioned, each D i can be represented by a D-quiver as the collection of rooted trees. For example, the D-quiver for a differential operator is Applying it to A EYM n,12 leads to . This example contains all crucial points we want to clarify, so let us give more explanations, especially about the similarity between shuffle structure in (4.23) and the rooted tree structure in (4.22).
• Firstly, let us consider the tree with root k 1 . It is connected to two branches {h 1 : {h 2 , h 4 }; h 3 } and {h 5 , h 6 }. Applying T 1h 1 2 and T 1h 5 2 will produce the structure where the subscript R is denoted for a "restricted shuffle", meaning that when making shuffle permutation for three sets, the first element of the third set should be placed after the first element of other two sets. Applying T h 5 h 6 n from the first branch will give us {h 5 , h 6 } as while applying insertion operators from the second branch will give (4.26) • Second let us consider the rooted tree with root K 4 , which also contains two branches. Applying T 4h 8 5 and T 4h 9 5 on the sub-structure {2, 3, ..., n − 1} R in (4.26) results in Having defined the differential operator D i for a vector in gauge invariant basis as in (4.20), we can apply it to the equation (4.1) and get a linear equation for expansion coefficient of a particular B i . However, for a vector with pseudo-loops, in general we will meet D i [B j ] = 0 for some j = i. In this case, we get a set of linear equations which can be solved directly. For a EYM amplitude with a large number of gravitons and gluons, the size of linear equations would become too large to be solved. Thus it is better to find a way to avoid solving linear equations. Before doing so, let us take a careful look on the linear equations, i.e., equations D i [B j ] = 0 with different B j 's under the same D i .
By inspecting D-quivers and corresponding operators, we find that there are two types of mixing-up for basis quivers. This conclusion comes from a key point that, while operators T ah i (a+1) or T h j h i n is able to select a particular dashed line or solid line uniquely in the basis quiver, the operator (k 1 · k h i )T 1h i 2 fails to do so. As a consequence, the contributions of different basis quivers will mix together when they are mapped to the same D-quiver. The reason why this situation can happen is that each pseudo-loop has two possible ways of mapping, each of which generates a D-quiver, so it's possible that two pseudo-loops share the same D-quiver. For example, let us consider the following four basis quivers B i which in total map to five D-quivers, Hence if we map the basis quiver B 1 to D-quiver D 2 and construct the corresponding differential operator D 2 , then applying D 2 on B 2 will produce non-zero result, which mix up the coefficients of B 1 , B 2 . Above example can be described in a more general way. Considering a quiver basis with pseudo-loop F hα 2i−1 hα 2i and the mapping (k 1 · k hα 2i )T hα 2i hα 2i−1 n T 1hα 2i 2 , if in the quiver basis there is also a factor F hα 2i h β , then we can always find a new basis quiver obtained by replacing from the original basis quiver 19 , which has non-zero contribution under above differential operator. We can do the replacement independently for each pseudo-loop. So if there are κ i solid lines connecting to the node h α 2i , the total number of new quiver basis which have non-zero contribution under D i will be ( p i=1 (κ i + 1) − 1). These quiver basis have an important property that D i [B j ] = 1 for all j, and it will be important in the later construction of linear combination of D i 's.
Then let us consider the second type of mixing-up originating from identity (4.16). Although these basis will not produce the same D-quiver 20 , they could produce the same ( k)-quiver by replacing a dashed line ( · K a ) or a solid line ( · k h j ) with ( · k 1 ). For example, applying D 1 on the following two basis quivers all yields non-zero results, .
(4.29) 19 Such operation of replacing can be simply realized by exchanging two subscripts hα 2i−1 and h β . 20 Please recall that the collection of D-quivers is a subset of all ( k)-quivers.
Note that B 1 can itself be a rooted tree or a branch of a rooted tree, while B 2 can only be a branch of a rooted tree. Thus in this type, a pseudo-lop branch has been mixed up with a tree branch. Explicitly, for quiver basis with pseudo-loop F hα 2i−1 hα 2i with mapping rule (k 1 · k hα 2i )T hα 2i hα 2i−1 n T 1hα 2i 2 , we can always find new quiver basis by replacing with arbitrary p = α 2i 21 . Since each pseudo-loop can do the replacement independently, there are in total (2 p − 1)(n − 2 + m − 1) new quiver basis. For the generated quiver basis, applying D i on them would produce a factor (−k hα 2i · (k 1 + K a )) or (−k hα 2i · k hp ) respectively according to (4.17). This is consistent with the mass dimension counting. Furthermore, the generated quiver basis have corresponding differential operators (4.29) which vanishes the original quiver basis with pseudo-loop. Thus the mixing between pseudo-loops and rooted trees will not be a problem if we solve the linear equation in a proper order.
We have discussed the mixing of basis under given differential operators, and let us continue to discuss how to disentangle them. The trouble comes from the first type of mixing-up, where a differential operator is not able to select only one vector. Our solution is to construct a linear combination of differential operators such that it can be applied to vanish all but one vector. Let us start from a simple example by considering the four basis quivers (4.28) and their corresponding D-quivers. It is easy to see that  It means D i uniquely selects a vector in gauge invariant basis. Generalizing this example, we can construct the linear combination of differential operators as follows.
• For a given basis B i , it could be mapped to many D-quivers. We randomly choose one D-quiver. For example, • For the branch connecting to root k 1 , there are two nodes originating from the corresponding pseudoloop. One node h a is connected to k 1 by a dashed line and we denote the other node by h b . We can separate this branch into two parts by disconnecting node h a from node h b while keeping the connecting line to h b , and denote these two parts by H a , H b . For example, • In graph H a , k 1 is connected to h a , and the arrows in solid and dashed lines all run from leaves to root k 1 . We can alternatively draw all possible rooted trees by connecting k 1 to another node other than h a by necessarily changing the direction of arrows for certain solid lines so that the arrows still run from leaves to root k 1 . Each new rooted tree defines a multiplication of operators denoted by D Ha,j with j = 1, ..., k where k is the total number of nodes excluding h a in the rooted tree H a . Then we define the linear combination of rooted trees as where s(j) is the number of solid lines connecting node h j to the node h a . For example, • Multiplying D Ha with the operators corresponding to H b produces the required operator that selecting a particular branch connecting to root k 1 . For example we get the linear combination • A basis quiver would have many branches, and for each branch we similarly construct a linear combination of operators. Multiplying all these linear combinations with those from rooted trees with root K a 's, we get the required differential operator selecting a particular vector B i in gauge invariant basis relating to the first type of mixing-up.
We shall emphasize that, after expanding above defined differential operator we will get many D-quivers, while these D-quivers would possibly produce non-zero results when applying to the vectors relating to the second type of mixing-up. This again infers that we should solve coefficients of vectors with fewer pseudo-loops first. We also remark that, although we have constructed the linear combination of operators to select a particular basis of the first type of mixing-up, when the size of linear equations is small it is quite favorable to solve them directly using the differential operators defined in (4.20). The reason is that, while it is much simpler for computing coefficients of the first type using linear combination of differential operators, it complicates contributions from the second type since all the D-quivers after expanding the linear combination would possibly produce non-zero results for the basis of the second type.

Algorithm for the evaluation of expansion coefficients
After clarifying the structure of differential operators, the next step is to apply them to the computing of expansion coefficients for the generic expansion formula (3.59). For vectors of gauge invariant basis defined in (3.57), the algorithm is implemented order by order, starting from p = 0 to the largest value p. For a given p, we start from the largest r to r = 0. The value of p denotes the number of pseudo-loops in a vector, hence when p = 0 the basis quiver contains only solid and dashed lines without any pseudo-loop. Such vector can be mapped to an unique D-quiver representing the following differential operator, Recalling identities (4.11), (4.13) and (4.14), a vector B j is non-vanishing only when its D j -quiver is the same as that given by (4.33). Thus the differential operator (4.33) uniquely selects one vector in gauge invariant basis while all others vanish, and the expansion coefficient can be solved by an univariate linear equation. Furthermore, the differential operator is normalized to one, hence the expansion coefficient can be directly computed by applying differential operator (4.33) on the EYM amplitude, leading to  Note that T aγ i hγ i (a γ i +1 ) inserts h γ i between a γ i and a γ i + 1 relative to the color-ordering, while T h β i h β i n inserts h β i between n and another graviton h β i . Hence in the resulting Yang-Mills amplitudes, the legs h i 's can never appear in the positions before 2 or after n, and all Yang-Mills amplitudes are in the BCJ basis with legs 1, 2, n fixed. An example of evaluating (4.35) has been discussed in (4.23).
Heading to p = 1 case, the differential operator for vector with one pseudo-loop is defined as 22 , with indices following convention (3.58), and the total number of differential operators is In differential operators (4.36), the insertion operator T 1hα 2 2 will contribute a derivative ∂ hα 2 ·k 1 relating to momentum k 1 . In its quiver, there is only one branch with root k 1 , and as we have analyzed, applying these differential operators on vectors will produce non-zero results only if the D-quiver of vector contains only one or no branch with root k 1 . So all vectors with two or more pseudo-loops will vanish under (4.36). Furthermore, when applying (4.36) on vectors without pseudo-loop, there could be non-zero contribution. However it is not an issue since all coefficients of such vectors have been solved a priori by differential operators (4.33) and they enter into the linear equations as known parameters.
For vectors with one pseudo-loop, there are in general more than one vectors being non-vanishing under a specific differential operator (4.36), as shown in (4.30). So we need to apply a complete set of differential operators to generate an algebraic system of linear equations, and solving expansion coefficients from this algebraic system. Alternatively, we can also apply the differential operator constructed by rule (4.32), i.e., a special linear combination of differential operators in (4.36). Then an expansion coefficient can be determined by an univariate linear equation again. Nevertheless, we can compute the coefficient of vector with one pseudo-loop as, (4.37) where the summation in curly bracket represents a linear combination of differential operators constructed following the rule (4.32). Note that the insertion operator T 1hα 2 2 inserts h α 2 in between 1 and 2, so the resulting Yang-Mills amplitudes are no longer in the BCJ basis with legs 1, 2, n fixed. BCJ relations are required in this step to write all Yang-Mills amplitudes into BCJ basis. While Yang-Mills amplitudes from contributions of vectors with no pseudo-loops are still in BCJ basis. Now let us proceed to the vectors with p pseudo-loops. According to the same argument with one pseudo-loop, by applying corresponding differential operator, all vectors with (p + 1) or more pseudo-loops 22 As mentioned, using the simple rule (4.20) we might need to solve algebraic systems of linear equations. While using a more complicated combination of differential operators as (4.32), the algebraic system is decoupled to univariate linear equations. will vanish. While for different vectors with p pseudo-loops, a linear combination of differential operators constructed by rule (4.32) is able to uniquely select a vector from all other vectors with p pseudo-loops. However, these differential operators still produce non-zero results when applying on vectors with (p − 1) or fewer pseudo-loops. In order to solve these linear equations, all coefficients of vectors with (p − 1) or fewer pseudo-loops should be solved a priori and enter these linear equations as known parameters. This inspires us to solve linear equations order by order from vectors with p = 0 to p = m 2 pseudo-loops. The differential operators relating to vectors with p pseudo-loops in gauge invariant basis are given as, (4.38) with indices following convention (3.58), and the total number of differential operators is The expansion coefficients of vectors with p pseudo-loops reads, Again, the insertion operator p i=1 T 1hα 2i 2 inserts h α 2 i 's in between legs 1 and 2, and we need to rewrite the resulting Yang-Mills amplitudes into BCJ basis by BCJ relations.
The algorithm for evaluation of expansion coefficients can be summarized as follows, --Start of Algorithm --STEP 0: Apply differential operators (4.33) on EYM expansion formula (3.59) to generate (m + n − 3) m linear equations, and solve expansion coefficients from these equations 23 . The result is directly given by (4.35).  23 In fact, solving equations is not necessary in this step. The expansion coefficients have been uniquely determined by (4.35), and the remaining thing to do is to explicitly work out the differential operators on A EYM n,m according to (4.35).
coefficients from these equations, and rewrite Yang-Mills amplitude into BCJ basis by BCJ relation. . . . STEP m 2 : Repeat the previous step but with p = m 2 differential operators.
--End of Algorithm -- The total number of repeated steps in the algorithm depends on the number of gravitons but not the gluons, while the total number of equations is much more sensitive to m than to n. Table 1 shows the number of linear equations to be solved in the algorithm for some EYM amplitudes. Comparing the total number of equations for a fixed m, for example A EYM 5,4 and A EYM 15,4 , we see the latter is about 44 times larger than the former when number of gluon increases ten. While comparing the total number of equations for a fixed n, for example A EYM 5,4 and A EYM 5,8 , we see the latter is about 85902 times larger than the former when number of graviton increases four. Hence the size of algebraic system is significantly controlled by m. One also notice that the number of equations decreases rapidly as moving to the next step in the algorithm. A large amount of equations are solved in Step-0, where expansion coefficients are explicitly defined by acting differential operators on EYM amplitudes. So in some sense it is trivial. For step p = 0 in the algorithm, the number of equations decreases significantly compared to the previous step, however nontrivial contributions from previous steps and BCJ relations would make results involving. Nevertheless, in each step the linear equation system is decoupled, and an expansion coefficient is trivially solved via an univariate linear equation. Step-0 1296 14641 65536 262144 4826809 34012224 100000000 2562890625 25600000000 Step -1  216  726  1536  61440  428415  1574640  28000000  318937500  1792000000 Step -2  3  3  3  2880  7605  14580  2100000  10631250  33600000 Step -3  0  0  0  15  15  15  42000  94500  168000 Step

Demonstration of EYM amplitude expansion in gauge invariant vector space
In order to demonstrate the EYM amplitude expansion in gauge invariant basis and the algorithm for determining expansion coefficients, in this section we present the expansion of EYM amplitudes with up to four gravitons. Expansion of EYM amplitudes with one, two and three gravitons to Yang-Mills amplitudes in BCJ basis has been discussed in paper [21], however here it receives a more systematic analysis in the language of gauge invariant vector space. While expansion of EYM amplitude with four gravitons to Yang-Mills amplitudes in BCJ basis is a new result.

The expansion of EYM amplitude with one and two gravitons
Let us start with A EYM n,1 (1, . . . , n; h 1 ). This amplitude lives in the gauge invariant vector space W n+1,1 , and the dimension of this space is (n − 2) according to (3.26). Hence A EYM n,1 can be expanded in a complete set of gauge invariant basis with (n − 2) gauge invariant vectors, as The expansion coefficient according to (4.35) is calculated as, . . . , a 1 , h 1 , a 1 + 1, . . . , n) , where the graviton h 1 is transformed to a gluon and inserted between a 1 , a 1 + 1 by T a 1 h 1 (a 1 +1) . Hence In comparison with the result in [21], we can reformulate above result as, where the shuffle permutation ¡ is defined in (B.2) and Y p as well as X p are defined in (B.4).
Let us continue to A EYM n,2 (1, . . . , n; h 1 , h 2 ). The dimension of gauge invariant vector space W n+2,2 is dim W n+2,2 = (n − 1) 2 + 1. The vectors in gauge invariant basis and their quiver representations are shown below as, is a real loop and should be excluded from the basis, while there is only one vector with pseudo-loop. Following the algorithm, Step-0 is to compute the coefficients of expansion basis with no pseudo-loops, i.e., by formula (4.35). Applying differential operators T a 1 h 1 (a 1 +1) T a 2 h 2 (a 2 +1) and T hσ 2 hσ 1 n T aσ 2 hσ 2 (aσ 2 +1) on A EYM n,2 respectively leads to A YM n+2 (1, . . . , a σ 1 , h σ 1 , a σ 1 + 1, . . . , a σ 2 , h σ 2 , a σ 2 + 1, . . . , n) , a σ 1 < a σ 2 σ∈S 2 A YM n+2 (1, . . . , a 1 , h σ 1 , h σ 2 , a 1 + 1, . . . , n) , a 1 = a 2 , (5.5) where σ = {σ 1 , σ 2 } is a permutation of {1, 2}, and the summation is over all elements of S 2 . In Step-1, we substitute above solutions back to the expansion formula and get, and there is only one unknown variable C[F h 1 h 2 ]. If applying differential operator (k 1 · k h 2 )T h 2 h 1 n T 1h 2 2 on both sides of above formula, in the RHS the non-vanishing contribution comes from vectors F h 2 h 1 F a h 2 and F h 1 h 2 , and according to (4.17), (4.15) we get In the LHS we get, Then we arrive at Yang-Mills amplitudes in the second term is already in the BCJ basis with legs 1, 2, n fixed while those in the first term is not. So we need to rewrite the first term in BCJ basis as, with K a 1 ···am = 1≤i<j≤m k a i · k a j . Combining above results together, we finally obtain Summing over all expansion basis with corresponding coefficients (5.5), (5.6) and (5.11), we get the expected EYM amplitude expansion. In fact, all contributions of vectors with no pseudo-loops computed in Step-0 can be rearranged in a compact expression as, and as we shall see, this is a general property for EYM amplitudes with arbitrary gravitons. After rearrangement of terms, we can rewrite the expansion of EYM amplitude with two gravitons in a rather compact form as,
Then we proceed to Step-1, and compute the expansion coefficients for vectors with one pseudo-loop. After substituting solutions in Step-0 back to the expansion formula, we get where the first summation runs over all possible splitting {α 1 , α 2 , γ 1 } of {1, 2, 3}, while the second summation not only runs over all splitting {α 1 , α 2 , β 1 } but also all possible values of β 1 . Terms in the first summation correspond to the first three quivers with one pseudo-loop in Fig.3, while terms in the second summation correspond to the remaining six quivers. As mentioned, for a fixed value of p, we should start from terms with larger r, i.e., terms in the first summation. As argued in the previous section, when applying a defined differential operator, only the corresponding vector survives and all others vanish. This means there is no mixing contributions between different pseudo-loop of the first type. For example, applying differential operator T h 3 h 2 n T 1h 3 2 T a 1 h 1 (a 1 +1) on formula (5.14), the only surviving vector with one pseudo-loop is F h 2 h 3 F a 1 h 1 . However vectors with no pseudo-loops would contribute, and from our previous general argument we can determine the non-vanishing vectors to be F h 3 where the relation (4.17) has been used. Working it out explicitly, we get .., a 1 , h 1 , a 1 + 1, ..., n − 1}, n) Terms in the first and second lines are similar to the one given in (5.9), hence we can borrow the result (5.11) to here and immediately work out the summation as, The other two terms with r = 1 can be obtained by permutation of above result.
Now we move to the vectors with p = 1, r = 0. As discussed, a defined differential operator (4.36) would possibly mix contributions of many vectors with one pseudo-loop, and in general we should solve an algebraic system of linear equations to compute all of them. However, in the current simple example we can intentionally choose a differential operator to avoid the mixing of vectors. For instance, in order to compute the coefficient of vector F h 2 h 3 F h 2 h 1 we should choose the differential operator T h 3 h 2 n T 1h 3 2 T h 2 h 1 n . If instead we choose the other differential operator T h 2 h 3 n T 1h 2 would be non-vanishing and their contributions will mix together. Hence we apply T h 3 h 2 n T 1h 3 2 T h 2 h 1 n on formula (5.14), and compute the coefficient as, Yang-Mills amplitudes in the second term are already in BCJ basis with legs 1, 2, n fixed, while those in the first term should be rewritten to BCJ basis by applying BCJ relations. Similar computations can be inferred from (5.9) and (5.18), and consequently all coefficients of vectors with one pseudo-loop can be computed. Summing up all above results we get the complete expansion of A EYM n,3 , which is consistent with results given in [21].

The expansion of EYM amplitude with four gravitons
EYM amplitude A EYM n,4 (1, . . . , n; h 1 , h 2 , h 3 , h 4 ) lives in gauge invariant vector space W n+4,4 , and it can be expanded as linear combination of dim W n+4,4 = (n+1) 4 +6(n+1) 2 +3 vectors. All vectors in gauge invariant basis and their quivers are shown in Fig.4. Among them, there are in total 6 × (n − 2) 2 + 44(n − 2) + 87 vectors with real loops which should be excluded. For the remaining vectors, we can compute their expansion coefficients following the algorithm. Again in Step-0, we compute the coefficients of vectors with no pseudo-loops by formula (4.35). We shall not write down the explicit coefficient for each basis but present the summation of them in a compact expression as 24 , Then let us continue with Step-1, to compute expansion coefficients of vectors with one pseudo-loop. As shown in Fig.4, there are in total seven distinct topologies, and the last one should be excluded. For the other six topologies, according to rules (4.20) we assign each of them with a differential operator respectively, and represent differential operators in quiver representation as , , , , , , 24 Note that the result of Step-0 can be similarly generalized to arbitrary points. where without ambiguity we have ignored the dashed line (h a 1) corresponding to (k 1 k ha )T 1ha2 , which is always linked to the ending point of the cyan line. The first two quivers of differential operators are consistent with the rules (4.19), and they are sufficient to distinguish the corresponding vectors uniquely. For the third and fourth quivers of differential operators, noticing the choice of direction of cyan line we know that they are also able to determine the expansion coefficients without mixing contributions from other vectors with one pseudo-loop. However, the last two types of vectors do mix together under the defined differential operators. It can be seen that, with the sixth quiver of differential operators it is able to distinguish the sixth type of vectors. However with the fifth quiver of differential operators, contributions from the fifth type of vectors would be mixed up with those from the sixth type of vectors. Although we can disentangle all vectors by constructing linear combination of differential operators as in formula (4.32), in the current simple example we have alternative way of solving equations. By firstly solving the coefficients of vectors of the sixth topology and then solving the vectors of the fifth topology but with the former solutions as known inputs, we are able to compute all coefficients order by order. Furthermore, we want to emphasize that, the differential operators also pick up contributions from vectors with no pseudo-loops, and we should compute all coefficients of vectors with no pseudo-loops before computing of vectors with one pseudo-loop.
Let us analyze these six topologies one by one. For the first topology, the corresponding differential operator also picks up following contributions in Step-0, , .
For instance, using differential operator ( Applying differential operator on A EYM n,4 produces Yang-Mills amplitudes and using BCJ relations they can be rewritten into BCJ basis. For the second topology, the corresponding differential operator picks up following contributions in , .
For instance, using differential operator ( Applying differential operator on A EYM n,4 produces Yang-Mills amplitudes For the third topology, the corresponding differential operator picks up following contributions in Step-0,

, .
For instance, using differential operator ( Applying differential operator on A EYM n,4 produces Yang-Mills amplitudes For the fourth topology, the corresponding differential operator picks up following contributions in Step-0, . For instance, using differential operator ( Applying differential operator on A EYM n,4 produces Yang-Mills amplitudes According to our discussion, we will consider the sixth topology before the fifth. The corresponding differential operator picks up following contributions in Step-0, . For instance, using differential operator (k 1 · k h 4 )T h 4 h 3 n T 1h 4 2 T h 2 h 1 n T h 3 h 2 n we can compute the coefficient of Applying differential operator on A EYM n,4 produces Yang-Mills amplitudes Then come to the last piece. Besides the contribution from the sixth topopoly, the differential operator corresponding to the fifth topology also picks up following contributions in Step-0, . Let's consider an example, the differential operator ( Applying differential operator on A EYM n,4 produces Yang-Mills amplitudes Above computations provide all expansion coefficients for gauge invariant basis with one pseudo-loop based on the solutions in Step-0 and the BCJ relations. Let us continue to Step-2, where there are only three different vectors According to the rule, we define differential operators for them respectively as It can be checked directly that each differential operator picks up only one vector with two pseudo-loops, while it also picks up following contributions in Step-0 and Step-1, , , For instance, Using differential operator is the contribution from expansion in Step-0, and is the contribution from expansion in Step-1. While applying differential operator on A EYM n,4 produces Yang-Mills amplitudes, Then using BCJ relations for A YM (1, α 1 , 2, . . . , n), A YM (1, α 1 , α 2 , 2, . . . , n), A YM (1, α 1 , α 2 , α 3 , 2, . . . , n) and A YM (1, α 1 , α 2 , α 3 , α 4 , 2, . . . , n) we can rewrite all Yang-Mills amplitude into BCJ basis with legs 1, 2, n fixed. Collecting all above results, we get the required EYM amplitude expansion. Because the final result is complicated we would not present the explicit expression for A EYM n,4 (1, · · · , n; h 1 , h 2 , h 3 , h 4 ). However we have numerically checked the algorithm up to A EYM 6,4 and find agreement with CHY formalism.

Conclusion
There are already quite a lot well-formulated results for expansion of EYM amplitudes to Yang-Mills amplitudes in KK basis, however a compact expression or even a recursive formula for expansion to Yang-Mills amplitudes in BCJ basis is still in pursuit. The latter expansion is generally much more complicated as conventionally expected. In the KK basis the expansion coefficients of Yang-Mills amplitudes are only polynomials of polarizations and momenta, and they are constrained to explicit compact expressions by gauge invariance. In the BCJ basis, the expansion coefficients of Yang-Mills amplitudes are instead rational functions, whose explicit form is much more difficult to determine. This is the reason that we consider using differential operators to determine expansion coefficients in paper [21]. This paper is motivated by the problem of expanding EYM amplitudes to Yang-Mills amplitudes in BCJ basis by differential operators. We have implemented an algorithm to systematically perform the expansion and compute the expansion coefficients. However the EYM amplitude is not directly expanded to BCJ basis but instead to a basis in gauge invariant vector space, as schematically shown in formula (4.1). After determining the expansion coefficients, we transform Yang-Mills amplitudes to BCJ basis by BCJ relations. Expanding EYM amplitude in a manifest gauge invariant form for both expansion basis and their coefficients is a very interesting point of view, and differential operators can be naturally introduced into the problem. It contributes to our major results.
The first major part of this paper is devoted to the construction of gauge invariant basis and their corresponding differential operators. A systematic algorithm is built upon the properties of applying differential operators on different basis. To construct a complete set of manifestly gauge invariant polynomials as the expansion basis, we start from the most general vector space V n,m with m ≤ n, where all possible polynomials of Lorentz contractions among polarizations and momenta live in this space, obeying some additional conditions. Then we define some linear mapping G i , which is a realization of gauge invariant condition for a polarization. By taking the interaction of kernels of all possible G i 's, we construct the gauge invariant sub-space W n,m from V n,m , which is the vector space containing all gauge invariant polynomials. This is also the space where the expansion basis of EYM amplitude lives. We present the formula for computing the dimension of W n,m , which indicates the number of gauge invariant vectors a EYM amplitude would be expanded to. We also find that the gauge invariant vectors can be realized by linear combinations of multiplications of fundamental f -terms. Above results at the end help us to construct a linearly independent and complete basis combinatorially for EYM amplitude expansion.
After clarifying the structure of gauge invariant expansion basis, we further construct differential operators from multiplication of insertion operators. The differential operators are constructed such that when applying a differential operator on an expression only one particular vector in gauge invariant basis is non-vanishing while all others vanishing. In order to do so, we start with analyzing the structures of gauge invariant basis and find the quiver representation for them. With the help of quiver representation, we summarize all possible components appearing in gauge invariant vectors, and provide mapping rules for writing a differential operator directly from a gauge invariant vector, as multiplication of three basic types of insertion operators. Based on above results, an algorithm for expansion of EYM amplitudes is implemented, with the idea of solving algebraic systems of linear equations order by order. To demonstrate the algorithm, we present the expansions of EYM amplitudes with up to four gravitons in the language of gauge invariant basis, which are all consistent with CHY formalism numerically.
Although the algorithm for expanding tree-level single-trace EYM amplitude to Yang-Mills amplitudes in BCJ basis has been laid down thoroughly in this paper, it still inspires further works to do in future. Firstly, the expansion coefficients of BCJ basis demands an explicit and possibly compact formulation. It is a rather difficult problem, but we have found some clues in results (5.9) and (5.18) already, and hope it could help to figure out the general picture. Secondly, in this paper we only deal with single-trace EYM amplitudes, while discussions can be generalized to multi-trace EYM amplitudes by using trace operator T i j . We think this generalization should be straightforward.
Thirdly, in this paper we are focusing on EYM amplitudes, so the parameters of vector space V n,m is constrained to m < n. However, the case m = n is also very interesting in physics since Yang-Mills amplitudes live in this space. Another interesting example is the deformed Yang-Mills theory with F 3 term [34,35]. Although the dimension of W n,m still holds for m = n, the explicit form of vectors in gauge invariant basis should be reconsidered since we are not able to trivially exclude momentum k n in all expression by momentum conservation. Furthermore, for Yang-Mills amplitude an additional constraint should be applied to the vector space, i.e., there should be at least one ( · ) contraction, and let us denote the vector space by W n,m . The new vector space W n,m can help us to understand the implication of gauge invariance in Yang-Mills amplitudes more deeply, along the line of former studies in papers [25][26][27]. It is also a curious problem about how to write Yang-Mills amplitudes in a manifestly gauge invariant form. Maybe it can also help us to understand more about the Pfaffian in the integrand of CHY formula, and provide a new point of view for BCJ relations.
In order to do so, it is suffice to show The proof of Ker G 1 + Ker G 2 ⊆ Ker G 1 G 2 is trivial. For each v ∈ Ker G 1 + Ker G 2 , it can always be written as Thus the action of G 1 G 2 on v is where we have used the commutative of G i , i.e., G 1 G 2 = G 2 G 1 . Hence v ∈ Ker G 1 G 2 , and consequently Ker G 1 + Ker G 2 ⊆ Ker G 1 G 2 . The proof of Ker G 1 + Ker G 2 ⊇ Ker G 1 G 2 is not so easy and we will prove it by induction. Let us start from the vector space V n,2 , i.e., containing only two polarizations 1 , 2 . A polynomial h n,2 in V n,2 can be written as where momentum conservation has been used to eliminate the appearance of k n . For h n,2 ∈ Ker G 1 G 2 , by imposing G 1 G 2 h n,2 = 0 we get From above equation we can solve α 1 and substitute it back to h n,2 . After reorganization of terms, we get Since the appearance of f i , it is easy to see that G i v i = 0. Hence v 1 ∈ Ker G 1 and v 2 ∈ Ker G 2 . This shows that if h n,2 ∈ Ker G 1 G 2 , there is also h n,2 ∈ Ker G 1 + Ker G 2 . Now let us assume that for all vector spaces V n,s , s < m, if a polynomial h n,s ∈ Ker G 1 G 2 , then it can always be separated into two parts, one part belonging to Ker G 1 and the other belonging to Ker G 2 . For a polynomial h n,m in the vector space V n,m , it can be expanded to where T mi ∈ V n,m−2 and i · T mi , T mi ∈ V n,m−1 . For h n,m ∈ Ker G 1 G 2 , by definition we have 0 = h n,m | 1 →k 1 2 →k 2 = ( m · k 1 )T (2) m1 + ( m · k 2 )T The result (A.10) tells us that all T mi , T mi , i = 3, . . . , m−1 and T mi , i = m+1, . . . , n−1 belong to Ker G 1 G 2 , and by the induction they belong to Ker G 1 + Ker G 2 . For the remaining terms in (A.7), i.e., h n,m = ( m · 1 )T m1 + ( m · 2 )T m2 + ( m · k 1 )( 1 · T m1 ) + ( m · k 2 )( 2 · T m2 ) . (A.11) After adding 0 = ( m · 1 )(k 1 · T m1 ) − ( m · 1 )(k 1 · T m1 ) + ( m · 2 )(k 2 · T m2 ) − ( m · 2 )(k 2 · T m2 ) at the RHS of above equation, we can reorganize h n,m to be h n,m = ( m · 2 )(T m2 + k 2 · T m2 ) + ( m · f 1 · T m1 ) + ( m · 1 )(T m1 + k 1 · T m1 ) + ( m · f 2 · T m2 ) . (A.12) Using the result (A.9) we get Thus h n,m belongs to Ker G 1 + Ker G 2 . So finally we have proven that Ker G 1 + Ker G 2 ⊇ Ker G 1 G 2 is valid in any vector space V n,m , and the proposition 1 is proven.

B Explicit BCJ coefficients
In this appendix, we provide some explanation for notations in (2.2). For convenience we also collect some explicit BCJ coefficients which are used in the computation. In formula (2.2), we have The K is defined as Definition of W needs further explanations. Given two ordered sets Ξ = {ξ 1 , ξ 2 , ..., ξ n } and β = {β 1 , ..., β r } where set β is a subset of Ξ, for a given element p ∈ Ξ with its position K in Ξ, i.e., ξ K = p, we define Furthermore, since p has split set β into two subsets β L p and β R p , i.e., the collections of elements on the LHS and RHS of p respectively, we can define Next we provide some examples. We consider the BCJ basis with legs 1, 2 being fixed in the first two positions and leg n in the last position in the color-ordering. For an arbitrary amplitude with one or two gluons inserted between legs 1, 2, we have the BCJ relations (k p · k 1 + k q · (Y q + k p ))(k p · (Y p + k q )) K 1pq K 1p A YM n+1 (1, 2, {3, . . . , n − 1} ¡ {q, p}, n)