Expansion of EYM amplitudes in gauge invariant vector space

Motivated by the problem of expanding the 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 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 a gauge invariant vector space by imposing constraints on the 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 this gauge invariant basis with expansion coefficients being linear combinations of the Yang-Mills amplitude, manifesting the gauge invariance of both the expansion basis and coefficients. With the help of quivers, we compute the expansion coefficients via differential operators and demonstrate the general expansion algorithm using several examples.


Mills amplitudes to
. Regarding the gravity amplitude, the Kawai-Lewellen-Tye (KLT) relations [19], which originally state that a closed string amplitude is a combination of the products of two open string amplitudes, degenerate to similar relations between gravity and Yang-Mills amplitudes in the field theory limit. In addition, the BCJ double copy conjecture reveals another new way of constructing the gravity amplitude from Yang-Mills amplitudes based on the exciting idea of color-kinematic duality [17,20,21].
In addition to these relations, the amplitudes of Einstein-Yang-Mills (EYM) theories where gravitons are allowed to interact with gauge bosons have also been investigated from many aspects [13,14,22,23]. In particular, in [14], a generalized KLT relation is proposed from the study of Cachazo-He-Yuan (CHY) formalism [10][11][12][13][14], schematically formulated for the tree-level single-trace EYM amplitude 1) as with being the amplitudes of Yang-Mills-scalar theory and the momentum kernel, defined in [24][25][26]. Parallel to the study of monodromy relations in string theory, in [27], the authors present a new relation formulating the EYM amplitude with n gluons and one graviton as a linear combination of -point Yang-Mills amplitudes in a compact expression. This result has been generalized to situations with more than one graviton [28,29] and double color traces [28] in the framework of CHY formalism. Furthermore, in [30], by studying the constraints of gauge invariance, a compact recursive formula is presented for the expansion of EYM amplitudes with m gravitons in terms of the KK basis of color-ordered Yang-Mills amplitudes; the result has also been proven in the CHY formalism [31] and generalized to multi-trace amplitudes [32]. For the purpose of the current paper, we recall the expansion of EYM amplitudes to color-ordered Yang-Mills amplitudes in the KK basis, as in [30,32], where is a set of m gravitons, and stands for the shuffle permutations between two ordered sets , i.e., permutations of keeping the respective orderings of and . In this expansion, legs and n are always fixed in the first and last positions in the color-ordering. Hence, using the recursive formula, at the end, the EYM amplitude would be expanded to the basis of Yang-Mills amplitudes with legs and n being fixed. The coefficient of each Yang-Mills amplitude is a linear combination of , which are polynomial functions of polarization vectors and momenta whose precise definition can be found in [30].

S[σ|σ]
A R S n−3 While the expansion of the EYM amplitude in the KK basis of Yang-Mills amplitudes has been solved completely, as the KK basis is not the minimal basis of colorordered 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? At first glance, it seems that this question has already been solved by the generalized KLT relation (1.1). However, in (1.1), the momentum kernel and are difficult to compute; we also need to sum over all permutations. Hence, the generalized KLT relation does not work well with respect to practical computation. One could also start with expression (1.2) and reformulate the KK basis to the BCJ basis using the BCJ relations. However, computation of several examples is sufficient 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 [33], a new method is proposed by introducing differential operators into this problem. The differential operator was originally applied to the study of the relationships among the amplitudes of different theories [34], and later, a series of studies showed how to apply differential operators to the expansion of the EYM amplitude to the KK basis [33,35,36]. Then, differential operators were naturally applied to the expansion of the EYM amplitude into the BCJ basis, being limited to some simple cases where the EYM amplitudes contain one, two, or three gravitons. However, a systematic method for a generic EYM amplitude with n gluons and m gravitons is still needed.
In this study, we attempt to fulfill this request by providing a systematic method for computing the expansion coefficients of the EYM amplitude with m gravitons in the BCJ basis. In addition to the use of differential operators, we also require the principle of gauge invariance. Because the Yang-Mills amplitudes of the BCJ basis are linearly independent, if we can write an EYM amplitude as a linear combination of Yang-Mills amplitudes of the BCJ basis, the gauge invariance of polarization tensors of gravitons would be transformed partially into the gauge invariance of expansion coefficients, which contain one half of the polarization vectors of the polarization tensors. Hence, the gauge invariance places strong constraints on the form of the expansion coefficients. In fact, the gauge invariance principle has already played an important role in the study of scattering amplitude. It is expected that the gauge invariance could completely determine the amplitudes of certain field theories [8,37,38], and further exploration can be found from various perspectives [30,34,[39][40][41][42]. In particular, as demonstrated in [30], it is the constraints of gauge invariance that make a compact formula available for expansion of EYM amplitude in the KK basis. However, the potential applications of gauge invariance have still not been fully exploited. In this paper, we would like to propose a different understanding of gauge invariance. Similar to what we have performed Chinese Physics C Vol. 44, No. 12 (2020) 123104 1) Hereafter we will always abbreviate tree-level single-trace EYM amplitude as EYM amplitude for simplicity.

123104-2
for the symmetries in the amplitudes of super-Yang-Mills theory, because the principle of gauge invariance is a strong constraint for gauge theory, we prefer to make it manifest at the level of scattering amplitudes.
With the new understanding of gauge invariance, in this study, we will show how to expand the general EYM amplitude into the BCJ basis of Yang-Mills amplitudes systematically. This paper is organized as follows. In §2, we review some background. In §3, we introduce the gauge invariant vector space living in a general vector space consisting of 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 representation, which is the description of the mathematical structures of these vectors and operators. With the help of quivers, we implement a systematic algorithm to compute expansion coefficients. In §5, we illustrate our method using several explicit examples, e.g., the EYM amplitudes with up to four gravitons for the purpose of clarifying some subtleties. In §6, we conclude our discussion and point out some problems that remain to be solved in the future. Detailed proofs of some propositions as well as some explicit BCJ coefficients in the BCJ relations are presented in the appendices.

Expansion of EYM amplitudes to Yang-Mills amplitudes in BCJ basis
In this section, we review some background knowledge, which is useful in the subsequent discussion of expanding the EYM amplitude to the BCJ basis of Yang-Mills amplitudes. First, as reviewed in [33], an arbitrary color-ordered Yang-Mills amplitude can be expanded to the BCJ basis with three particles being fixed in certain positions related to the color-ordering, as follows: A n (1, β 1 , · · · , β r , 2, α 1 , · · · , α n−r−3 , n) = The expansion coefficients, i.e., the BCJ coefficients, were first conjectured in [17] and later proven in [18], using the expression Second, we review the differential operators that are originally introduced in [34]. An important differential operator is the insertion operator, defined as Physically, it represents changing a graviton k into a gluon and inserting it between i and in the color ordering of gluons. If two gluons are not adjacent, for instance , we will have 4) and its physical meaning is also clear 1) . Another important operator is the gauge invariance differential operator, defined as . (2.5) It has a physical meaning of imposing gauge invariance, i.e., changing . For an arbitrary polynomial of polarization vectors and momenta, if it vanishes under oper- ator , we can conclude that it is gauge invariant for polarization vector . Gauge invariance operators are commutative, i.e., , so the result of a multiplication of sequential operators does not depend on the ordering, and we can denote a sequential gauge invariance operator as The insertion operator and gauge invariance operator satisfy the following commutative relation, with , and it is valid after application to any functions of polarization vectors and momenta 2) .
Finally, let us present a general discussion on the expansion of the EYM amplitude to the BCJ basis. For particles with spin, the corresponding Lorentz representations are constructed by polarizations, e.g., the polarization vector for gluons and polarization tensor for gravitons. When expanding the EYM amplitude to the BCJ basis, the polarization tensor of gravitons is factorized into two parts, . The part is inherent because of the polarization vector of gluons in the Yang-Mills basis, while the other part, , is absorbed into the expansion coefficients. More explicitly, the expansion coefficients are rational functions of momenta and polarization vectors . A crucial difference between expanding to the KK basis Chinese Physics C Vol. 44, No. 12 (2020) 123104 are not in the same trace, it has no clear physical meaning. 2) For detailed description of these differential operators and their relations please refer to the paper mentioned before.

123104-3
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., This observation inspires us to consider another form of expansion: In equation (2.8), independent Yang-Mills amplitudes are taken to be the expansion basis, and each coefficient as a function of momenta and polarization vectors should satisfy the conditions of gauge invariance for all with . In equation (2.9), represents the expansion basis, and the expansion coefficients become a linear combination of , with coefficients being rational functions of momenta. The latter form has already appeared in [33]; to distinguish the two different bases, we call gauge invariant building blocks 1) .

Building up expansion basis in gauge invariant vector space
As mentioned earlier, in the expansion of EYM amplitudes, the gauge invariant coefficients as well as expansion basis are crucial. They are polynomial functions of polarization vectors that vanish under conditions of gauge invariance. In this section, we would like to start from the most general vector space and localize a gauge invariant subspace of it. The expansion basis we are looking for is living in this exact subspace.

Gauge invariant vector space and its dimension
Let us start from the most general polynomial , constructed by Lorentz contractions of n momenta and m polarizations with . By Lorentz invariance and the multi-linearity of , this polynomial must be formed schematically as where for each monomial, the degree of is m, and each appears once and only once, while the coefficients are rational functions of Mandelstam variables of momenta. If we take all monomials as a generating set, 2) then we can This operator establishes a linear mapping between different vector spaces as where in the resulting vector space, the polarization does appear and is replaced by , denoted . This linear map is surjective 3) by noticing the reduction of , i.e., We can successively apply different gauge invariant operators , with , and establish a mapping chain of vector spaces. Because all are commutative, the result does not depend on the ordering of successive application, and we can denote the mapping chain as The superscripts label the removed polarization vectors in the vector space. Note that different orderings of applying produce different mapping chains, which eventually lead to the same vector space, so (3.4) in fact represents a collection of mapping chains.
The kernel of the linear map is defined by Physically, this means that the vectors of kernel are gauge invariant for the i-th particle. Using the fact that the linear map is surjective (4.3), by the fundamental theorem of linear mapping [43], we have Then, the dimension of the kernel can be computed using the difference in the dimensions of vector space as When applying more than one , this relation can be generalized to Chinese Physics C Vol. 44, No. 12 (2020) 123104 1) Although we already know the formulation (2.9) is more suitable for applying differential operators, in the previous work we are not able to push the discussion further since the discussion of building blocks are too difficult at that time.
2) These monomials are not linearly independent. There are relations between them generating by momentum conservation and transverse condition . Furthermore, we consider only the parity even case, i.e., without total antisymmetric tensor .
3) The property of surjectivity is the cornerstone in our discussion. For the vector space of polynomials without term surjectivity of the map no longer holds. 123104-4 For example, let us consider the simplest case , Vector space is the field of rational functions of Mandelstam variables, so the basis is simply , and . For a vector space with only one polarization, the kernel consists of all vectors vanishing under the gauge invariant operator. This is the gauge invariant vector sub-space in a vector space . Thus, we have For a general vector space with m polarizations, we can define the gauge invariant vector sub-space as the intersection of kernels of all possible linear maps as This means that a vector in would vanish under any linear map . This is exactly the sub-space where all gauge invariant coefficients of (2.8) and the expansion basis of (2.9) live.
Let us attempt to compute the dimension of and start with the case . Generally, for any two linear spaces , we have the following relation for the dimension 1) , . (3.12) Applying this relation to the vector spaces of kernels, i.e., , we get dim W n,2 := dim (Ker . (3.13) The first two terms in the RHS can be computed using (3.7), and to compute the third term, we need to use the following proposition 2) , G i PROPOSITION 1: any two kernels of linear maps satisfy the splitting formula, 14) and its generalization, Recursively using (3.12), we need to generalize the above result to arbitrary m. For simplicity, let us denote , and when , we have In the second line, the first three terms have already been computed, while to compute the fourth term, we need to use the following proposition.

3)
G i PROPOSITION 2: three kernels of linear maps satisfy the distribution formula, and its generalization, G i PROPOSITION 2 EXTENDED: the kernels of linear maps satisfy the generalized distribution formula, Together with (3.12), we can rewrite (3.17) as In equation (3.20), to compute the dimension , we need the result of , which by proposition 1 extended (3.15) is equal to . Using (3.8), we get Then, Chinese Physics C Vol. 44, No. 12 (2020) 123104 Suppose are subspaces of V, then the sum of is defined as the set of all possible sums of elements of , explicitly . We should note that the definition of sum is different of direct sum, a sum is a direct sum if and only if , and for direct sum . 2) Proof of proposition 1 and proposition 2 can be found in Appendix A.
3) In general is not true. For example, in a two-dimension space U, let us choose to be line , and respectively. Then is the whole XY-plane, and is the line . While in the RHS, and are just the origin . So the RHS is a point.
Notice that the numerical factors are nothing but for . 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 obtain (3.23) 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 can be computed via In [37], the same result has been provided up to 2) .
Compared with that result, our calculation shows more efficiency than that shown by solving linear equations of gauge invariance directly. Furthermore, several examples of and with arbitrary n but definite values of m are listed below.

Gauge invariant vectors
The dimension of gauge invariant vector space characterizes the minimal number of vectors required to expand an arbitrary vector, while the explicit form of the vector is not constrained. From the working experiences of EYM amplitude expansion with one, two, and three gravitons [33], we get the insight that the coefficients appearing therein could be recast in a manifestly gauge invariant form as linear combinations of multiplication of fundamental f-terms. Here, the fundamental f-terms stand for two types of Lorentz contractions of field strength and external momenta, with at most two , This observation can be generalized beyond and can be stated as follows. For any vector in gauge invari- dim W n+m,m n − 2 (n − 1) 2 + 1 n 3 + 3n (n + 1) 4 + 6(n + 1) 2 + 3 Chinese Physics C Vol. 44, No. 12 (2020) 123104 1) The counting of (3.27) can be carried out as follows. Firstly we select i pairs of , and there are choices, while each left can be contracted with momenta after by momentum conservation. For 's, the number of different contractions is .
2) In this paper, there are two types of spaces being considered. The another one is the space with at least one contraction between polarization vectors in polynomials, i.e., polynomials without monomial , which is exactly the vector space that Yang-Mills amplitudes live in. Its dimension is . We shall prove this statement by induction. The cases with have already been shown to be true in [33]. Following the idea of induction, we assume that this statement is true for all and prove that it must be true for m.
where momentum conservation has been applied to eliminate , so that all appearing in are linearly independent. Polynomials and . Because , by definition, we have From the operator equation (2.7), we explicitly have with . Applying this to generates a set of equations as follows: where we have considered the fact that does not contain . With the above result, we can rewrite as We also need to consider the gauge invariance of with respect to the polarization vector , Then, we get After substituting the above results back into , we get Therefore, is already manifestly gauge invariant for polarization vector . In fact, we can also choose to eliminate other coefficients in (3.34) and introduce different poles in the denominator of .
We can also generate another set of equations by considering the operator relations with and . Applying them to produces which means that is gauge invariant for . By assumption of induction, can be written as a linear combination of the multiplication of fundamental fterms. Because and are gauge invariant for with , and are linearly independent, is also gauge invariant for all of its own polarization vectors. Again, by assumption of induction, any can also be written in a manifest gauge invariant form where only f appears. Thus, as a linear function of and , the polynomial 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 the above procedure to in (3.29) and rewrite it as where in the last summation, can be equal to . Let us again apply the operator equations , with and , which generates a set of equations, Therefore, becomes for its own polarization vectors and can be written as a linear combination of the multiplication of fundamental fterms. For the same reason as before, we conclude that is also gauge invariant for its own polarization vectors. Continuously applying the same procedure to until the last polarization vector, we arrive at To further reduce the expression to the fundamental f-terms, we can get help from the following identities, could be any string. More explicitly, applying the above identity to the expression with three f s, we get Therefore, any f-term with any number of can be reduced to fundamental f-terms, while at the same time, has been reorganized as a linear combination of multiplication of fundamental f-terms. This ends the proof of statement by induction.
Before ending this subsection, let us take a look at another gauge invariant f-term that is mentioned in [33], i.e., the trace . It can be expanded as where identity (3.42) has been used in the derivation. Combining the first and third terms, as well as the second and fourth terms, we obtain A simple example is . Therefore, this type of gauge invariant f-term, which is originally viewed as a new type different from , is also composed of the fundamental f-term.

Gauge invariant basis
Any gauge invariant vector in could be an element to form a gauge invariant basis in the EYM amplitude expansion (2.9). However, to turn a subset of 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 is equal 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 Therefore, one can always reduce any fundamental fterms to the following form, From the definition of , it is easy to obtain In the case of , the momentum list is , while the polarization vector list is , so by default, the above subscripts and . After using momentum conservation to eliminate , we can restrict the fundamental f-terms to be Using the 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. To demonstrate that they form a real basis of , we should show that the total number of these vectors is equal to 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 is Noting the relation we immediately have , as 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 in (2.9). In practice, we would prefer a basis with minimal dimensions; then, we can define the fundamental f-terms as The vectors in the expansion basis can be constructed from the above fundamental f-terms as , p, q, r ∈ N and 2p+q+r = m , with the convention .

(3.58)
They contribute to a complete set of the expansion basis, and a general EYM amplitude can be expanded into this basis as

59)
where is the set of gravitons excluding , and the three sets , , and with are a splitting of all gravitons.
represents a particular vector in the expansion basis , represents the coefficient of the corresponding vector, and the reduced summation runs over all possible splittings , where the prime denotes that terms with the index circle should be excluded. The more discussion of index circle can be found in [33]. We can see that all the information regarding polarization vectors is encoded in , 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 mentioned earlier, the EYM amplitude can be expanded schematically in the form, or more explicitly, as in (3.59). The expansion coefficients are linear combinations of Yang-Mills amplitudes . To use (4.1) efficiently, a crucial point is to find a way to distinguish vectors in the gauge invariant basis from each other. Inspired by the explicit form of vectors in (3.1), we note that the signature of vectors is the structure , where and could be linearly independent. This motivates us to consider two kinds of differential operators: Applying these operators to the RHS of (4.1), all terms will vanish except those containing corresponding and . While applying these operators to the LHS of (4.1), the physical meaning will be different. Applying to single-trace EYM amplitudes produces multi-trace EYM amplitudes, which would complicate the amplitude expansion; however, applying to a single-trace EYM amplitude produces another singletrace EYM amplitude but with one less graviton. will transform the graviton to a gluon and insert the gluon in the positions between gluons respecting the color-ordering. Therefore, each time an insertion operator is applied to (4.1), the number of gravitons is reduced by one; then, a multiplication of m insertion operators would transform the LHS of (4.1) to Yang-Mills amplitudes completely, as expected 1) . In fact, we can take one step further and define a differential operator as a multiplication of m properly chosen insertion operators. When applying the differential operator to (3.59), in addition to some vectors with already known coefficients, there would be one and only one vector with unknown coefficient remaining in the RHS of (3.59), and all other vectors vanish: As a consequence, we have a linear equation with only one unknown variable, and the corresponding expansion coefficient can be computed directly as a function of , generated by a differential operator applied to the RHS of (3.59) 2) . The problem of EYM amplitude expansion is then translated to construction of properly defined differential operators, which would be the major purpose of this section. Surprisingly, we find it very helpful to use quivers to represent the gauge invariant basis and differential operators for our purpose.

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 , so all other types of Lorentz contractions and can be treated as unrelated factors. To characterize the structure of in a gauge invariant vector, we can assign a quiver, i.e., directed graph, to it 3) . Therefore, in this subsection, we first define the quiver representation of atomic factors such as ; second, we give the quiver representations of fundamental f-terms; finally, we consider the quiver representation of gauge invariant vectors and discuss the properties of their quiver representation.
We call a directed graph representing all of a vector the -quiver of the vector. In a quiver, we use a directed solid line to represent with an arrow pointing to a graviton momentum and a directed dashed line to represent with an arrow pointing to a gluon momentum , as [45], [46]

4)
As for the fundamental f-term (3.48), which can be expanded as

4)
Chinese Physics C Vol. 44, No. 12 (2020) 123104 1) Alternatively, we could also apply less insertion operators to generate a set of linear equations of single-trace EYM amplitudes, and recursively use the expansion of single-trace EYM amplitude with less number of gravitons into Yang-Mills amplitudes.
2) The idea of selecting only one unknown variable at each step is similar with that of the OPP reduction method for one-loop amplitude.
ϵ · k 3) The idea of using arrows to represent Lorentz contractions has already been applied in previous literatures, where all types of Lorentz contractions are considered. However, we are only interested in Lorentz contraction of the type in this paper.
(ϵ · k) (ab) 4) From now on, we will identify an directed line with its corresponding term, and sometimes when we refer to a specific directed line connecting two nodes from a to b, we will use the label and b is called the head, a called the tail.

123104-10
(ϵk) (ϵk) its -quiver representation consists of three -directed graphs, as 1) Because each graph denotes a multiplication of terms, when applying the derivatives to , we will get non-vanishing results. Similarly, for , their -quivers are where we have distinguished two cases, being the momentum of a graviton and being a momentum of a gluon.
Note that the factor exists in both and , so the action of the derivative on them both is non-zero. Consequently, we prefer to eliminate the dashed lines representing in the graphs of -quivers to obtain a simple presentation. Furthermore, to represent one fundamental f-term by only one graph and distinguish from , we combine the two solid arrows in (4.5) into a loop. Finally, the fundamental f-terms , , defined in (3.54), (3.55), and (3.56) are represented by quivers in Fig. 1. To distinguish these quivers from -quivers, we will call them basis quivers or just quivers.
We should emphasize that from a basis quiver, it is easy to recover all corresponding -quivers by replacing any one solid or dashed arrow in the graph by a dashed arrow , i.e., from Fig. 1 to (4.5), (4.7). (ϵk) However, given an -quiver, it is difficult to determine which basis quiver it comes from, especially when there are many lines. The fact that there is no oneto-one correspondence between basis quivers and quivers causes some technical difficulties in the construction of differential operators. Fortunately, for a gauge invariant vector, its basis quiver and -quivers do possess a common property: they all contain m and only m lines (counting both dashed lines and solid lines), as each line carries one .
Note that the basis quiver for is a colored loop, where the colors are to remind us that it is an overlapping of three -quivers after eliminating dashed lines. We call such a colored loop a pseudo-loop. In general, there are also real loops. For example, the containing a monomial and containing a monomial can be represented as However, as explained in [33], the terms with indices or part of indices forming a closed circle will not be present in the expansion of the EYM amplitude, although such terms do appear in the gauge invariant basis. Therefore, we will exclude basis quivers with real loops in practical computation.
Next, let us consider the quiver representation of a vector in the gauge invariant basis. As shown in (3.57), such a vector is a multiplication of fundamental f-terms, expressed as Because each appears only once in a vector, only one Chinese Physics C Vol. 44, No. 12 (2020) 123104 Notice that there are four terms in the expansion of , while the term is the most crucial signature to distinguish it from other fundamental fterms. However in this paper we only consider insertion operators so that is out of our sight.

123104-11
exists; therefore, we can conclude that each point labelled by in the basis quiver of a gauge invariant vector has at most one out-going line but possibly several incoming lines. Consequently, all pseudo-loops are topologically disconnected from each other. The point labelled by is connected by only in-coming lines but not outgoing lines; hence, all such points are also topologically disconnected from each other. Furthermore, pseudo-loops cannot be connected with points labelled by either. Therefore, a quiver graph could have many disconnected components, whose number is at least p and at most , as several dashed lines can be connected to the same node , while a solid line for can be connected to one and only one disconnected component.
With the above analysis, let us discuss the possible structures appearing in a quiver representation for a vector in a gauge invariant basis (4.8). First, because each is represented by a dashed directed line with an arrow pointing to , its head can never be connected with a pseudo-loop or a solid line. Second, each is represented by a solid line with an arrow pointing to , , its head is linked with a dashed line, while if , its head is linked with a pseudo-loop, and if , for instance , its head is linked with another solid line, and the head of the latter is further linked with a pseudo-loop, a dashed line, or a solid line. A succession of solid lines should finally stop at a dashed line or a pseudo loop; otherwise, it would form a real loop that should be excluded.
To summarize, the quiver representation of a vector in a gauge invariant basis could contain the following substructures: 1. only a single dashed line, 2. a dashed line linked with a tree consisting of solid lines, 3. only a single pseudo-loop, 4. a pseudo-loop connected to a tree consisting of solid lines on one side, 5. a pseudo-loop connected to two trees consisting of solid lines on both sides, as shown in Fig. 2. Two examples of quiver representations for two vectors in the gauge invariant basis of are shown as follows: The two examples illustrate our previous discussion very well. There are three disconnected components for the first one and two for the second one. In the second graph, two dashed lines are connected to one node, representing the fundamental f-terms , . 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 using the concept of a rooted tree [47]. The quiver of a vector in a gauge invariant basis consists of some disconnected components, and each component contains only one pseudo-loop or a node . If we focus on a disconnected component with node , it is exactly a rooted tree, with the root being the node . More precisely, it is a directed rooted tree with an orientation toward the root, i.e., all lines in the tree are directed to the root from the leaves, as illustrated in the previous two examples. For the disconnected component with a pseudo-loop, we could split the pseudo-loop into two colored lines, resulting in two sub-graphs. For each sub-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. This picture of rooted trees will help us to construct the differential operators and understand many properties of our algorithm later.

Constructing differential operators
Because a vector in the gauge invariant basis is a polynomial of , it will be non-vanishing under the action of a derivative only if its -quiver representation contains a solid or dashed line corresponding to . Hence, by constructing a differential operator as a proper combination of some derivatives , we expect under its action that ideally only one vector is non-vanishing, so it can select a particular non-vanishing vector in a gauge invariant basis. Although in fact we cannot do this, we succeed in dividing the computation of coefficients of a gauge invariant basis into many steps, and in each step, by applying an appropriate differential operator, only one new vector is non-vanishing, except some vectors whose coefficients are already known. The goal in this subsection is to construct such differential operators.
The expected differential operators can be constructed using three types of insertion operators (4.1). The first type of insertion operator takes the form where is the momentum of a gluon. A vector is nonzero under if its -quiver contains a dashed line corresponding to or . Applying this insertion operator to the fundamental f-terms, we have 11) and The above results tell us that if the basis quiver of a vector in a gauge invariant basis contains a dashed line representing , then a differential operator containing the insertion operator will select out this vector and other vectors containing the same dashed line. The relation (4.12) can be graphically represented as,  14) represented in quivers as Because both and are non-vanishing under , we may conclude that this insertion operator is not sufficient to distinguish these two terms. However, we shall note that the insertion operator is actually a differential operator that works through the smaller pieces, i.e., Lorentz contractions , rather than the fundamental f-terms and . According to this view, it is easy to accept that and are non-vanishing under the action of , as their quivers both contain a solid line from to .
To construct a differential operator that can distinguish from , we need to consider a third type of The key difference in these two terms is that has two polarization vectors, while has only one. In other words, in the -quiver of , there are always two lines linked together, a solid line and a dashed line or , so we can multiply by an additional insertion operator containing the derivative ; under such operators, always vanishes. Then, choosing the operator and applying it to , we have It is easy to see that the operator satisfies our require- ment for distinguishing and , and it also distinguishes the pseudo-loop of from all other pseudoloops. However causes some additional difficulty, as there will be some multiplication of fundamental fterms that do ont vanish, such as (4.16) This means that although is able to distinguish one pseudo-loop from the others, it would mix contributions from vectors without pseudo-loops. However, this is not a problem at all, if we attempt to solve the coefficients of the basis in multiple steps. We can first compute the coefficients of and using the differential operators and respectively, under which has no contribution at all. Then, we can apply to compute the coefficient of , and we treat the coefficients of , as known input.
After the above discussion, we can roughly give a general picture of constructing a differential operator to select a particular vector in the gauge invariant basis through the quiver representation. The major idea is to construct a new special -quiver from a vector's basis quiver, which can be used to construct the expected differential operators. A reasonable method for determining these new -quivers is as follows: a dashed line in the basis quiver of a vector suggests that there is also a dashed line in the new -quiver, but representing and a solid line in the basis quiver also suggests that there is a solid line in the new -quiver representing , while for a pseudo-loop in the basis quiver, we can choose to construct either a solid line connected to a dashed line or a solid line connected to a dashed line in the new -quiver. We are free to take any one of the two choices when meeting a pseudo-loop. Finally, we obtain a new -quiver, which is used to construct differential operators. As we discussed at the end of the previous subsection, the -quiver is a collection of rooted trees. The disconnected component of a pseudo-loop in the basis quiver of a vector has been split into two branches, where each branch is a rooted tree with the root being and is suitable according to our choice, and the components without pseudo-loops directly give us rooted trees. Furthermore, a collection of rooted trees can be algebraically represented as an embedded structure, where at each level, we write . 2)

(ϵk)
Second, having obtained the desired -quivers, we can construct the corresponding differential operators using the following rules: 1. assign an operator to each dashed line in the new -quiver, which uniquely picks up the corresponding dashed line in a vector's basis quiver;

assign an operator
to each solid line in the new -quiver, which uniquely picks up the corresponding solid line in a vector's basis quiver; 3. assign an operator to each dashed line in the new -quiver. These rules can be represented graphically as (4.18) (ϵk) (ϵk) Therefore, the corresponding differential operator for a vector in a gauge invariant basis is defined by multiplying all assigned operators in the new -quiver together; then, we call the -quivers constructed according to the above rules D-quivers. We want to emphasize the following: (1) there is a one-to-one map between D-quivers (ϵk) and differential operators, so one quiver defines a unique differential operator; (2) a D-quiver is a special quiver, which can be associated with a given basis quiver.
Finally, the above discussion can be summarized as the following map, stating that from a given vector to a corresponding differential operator, where . There are several technical points we wish to explain. First, the mapping rule is defined such that Second, although insertion operators are commutative, when acting on EYM amplitudes, we need to choose a proper order to make the physical meaning clear. We shall apply insertion operators of the type first; then, we apply the types and . More explicitly, the ordering of applying insertion operators is from the roots to the leaves in the D-quiver opposite to the direction of the arrows. In fact, we can make the result more concrete when acting on . As mentioned, each can be represented by a D-quiver as a collection of rooted trees. For example, the D-quiver for a differential operator is Then, the rooted trees can be written as Applying this to leads to multiplied by . This example contains all crucial points we wish to clarify, so let us give more explanation, especially about the similarity between the shuffle structure in (4.22) and the rooted tree structure in (4.21).
• First, let us consider the tree with root . It is connected to two branches, and . Applying and will produce the structure where the subscript R denotes a "restricted shuffle," meaning that when making a shuffle permutation for three sets, the first element of the third set should be placed after the first element of the other two sets. Applying from the first branch will give us , as

Applications of differential operators
Having defined the corresponding differential operator for a vector in a gauge invariant basis as in (4.19), we can apply it to equation (4.1) and obtain a linear equation for the expansion coefficient of a particular , as well as other coefficients. However, for a vector with pseudo-loops, in general, we will meet for some . In this case, we have a set of linear equations. For an EYM amplitude with a large number of gravitons and gluons, the number of linear equations will become too large to be easily solved. Thus, it is better to find a way to avoid solving a large number of linear equations.
To find a such method, we need to analyze the behaviors of different under the action of , i.e., equations with different under the same . By inspecting D-quivers and corresponding operators, we find that there are two types of problems that cause difficulties in solving linear equations.
The first problem comes from a key observation that, while operator or is able to select a particular dashed line or solid line uniquely in the basis quiver, the operator fails to do so. As a consequence, the contributions of different basis quivers will mix together when they can produce the same D-quivers. The reason for this is that each pseudo-loop of the vectors' basis quiver has two possible ways of generating Dquivers, so it is possible that two basis quivers with pseudo-loops generate the same D-quiver. For example, let us consider the following four basis quivers , which generate five D-quivers in total.
Hence, if we choose as the corresponding differential operator of the basis quiver , then after applying to these five vectors, is also non-zero, in addition to , which means that the coefficients of are mixed together in the linear equation given by .
The above phenomenon is a general one. Assuming the basis quiver of a vector in the gauge invariant basis has a pseudo-loop connected with a solid line , and the corresponding differential operator of the pseudo-loop is , then we can almost always find a new vector in the basis having a factor 1) , which is non-zero under the same differential operator. We can do this operation independently for each pseudo-loop in a vector. If there are solid lines connected to the node , the total number of vectors that is non-zero under the corresponding differential operator of the pseudo-loop will be . The results of these vectors under the action of the differential operator are for being a vector of the set; this fact will be important in the subsequent construction of a linear combination of .
Now, let us consider the second problem originating from identity (4.16). Although the basis quivers of some vectors will not produce the same D-quiver 2) , they could give the same -quiver by replacing a dashed line or a solid line by . For example, applying to the following two basis quivers yields nonzero results, Note that can be a rooted tree by itself or a rooted tree obtained by splitting a pseudo-loop, while can only be a rooted tree obtained by splitting a pseudo-loop. Thus, in this case, a branch of a disconnected component with a pseudo-loop is mixed with a disconnected component without a pseudo-loop. Explicitly, for a vector with a pseudo-loop and the corresponding operator for the pseudo-loop , we can always find some new vectors by replacing with , or with arbitrary 3) . Because the replacement for each pseudo-loop is independent, there are a total of new vectors, and applying to these new vectors would produce or , respectively, according to (4.16). This is consistent with the counting of mass dimension. However, these new vectors have their corresponding differential operators (4.28) under which the original vector with a pseudo-loop vanishes. Thus, the second problem is easy to deal with if we solve the linear equations of unknown coefficients in the proper order.
We have discussed two types of problems in detail, and the second type is easily solved, so let us continue to discuss how to deal with the first one. The first type of problem originates from the fact that under the action of a differential operator, several vectors with pseudo-loops in the gauge invariant basis do not vanish at the same time, so their coefficients are mixed together in the linear equations. Our solution is to construct a linear combination of differential operators such that under its action, only one vector is non-vanishing. Let us start from the simple example (4.27), where it is easy to get If we define some new differential operators as , then This means selects a unique vector from the entangled vectors, and the linear equations of the coefficients of these vectors are easily solved. Generalizing this example, we can construct the linear combination of differential operators as follows.

B i
• For a given vector , we can obtain many Dquivers in general, but we choose only one D-quiver freely. For example, • For the D-quiver whose root is , there are two nodes coming from the original pseudo-loop. If the node connected to by a dashed line is denoted , then another node is denoted by . We can separate this D-  . Then, we define a new differential operator by where is the number of steps for moving from the node to the node . For example, • Multiplying by the differential operator corresponding to gives us the expected operator, which will select only one particular vector from the set of vectors entangled with the original vector. For example, we get the linear combination D D B i • A basis quiver of a vector would have many disconnected components, and for each disconnected component with a pseudo-loop, we can apply the same procedure and similarly construct a corresponding operator as a linear combination of some operators D. Multiplying all these operators with the operators obtained from disconnected components without pseudo-loops, we obtain the final differential operator, which will select a particular vector in a gauge invariant basis without the first type of problem.
We should emphasize that, after obtaining these differential operators using the above method, if we apply them to the expansion, there are still some issues resulting from the second type of problem. This suggests that we should solve coefficients of vectors with fewer pseudo-loops first. We also remark that, although we have provided the method to solve the problem of mixing of some vectors in solving the linear equations of coefficients, when the size of linear equations is small, it is quite favorable to solve them directly using the origin-al differential operators defined in (4.19). The reason is that, while it is much simpler for computing coefficients of the mixed vectors by using differential operators constructed by the above method, it may be complicated for the cases we meet in the second type of problems since some vectors with less pseudo-loops are non-vanishing under the action of these operators for the second type of problems.

Algorithm for evaluation of expansion coefficients
After clarifying the structure of differential operators, the next step is to apply them for computing the expansion coefficients for the generic expansion formula (3.59). For vectors of the gauge invariant basis defined in (3.57), the algorithm is implemented order by order, starting from to the largest value p. For a given p, we start from the largest r and proceed to . The value of p denotes the number of pseudo-loops in a vector; hence, when , the basis quiver contains only solid and dashed lines, without any pseudo-loops. Such a vector can be mapped to a unique D-quiver representing the following differential operator, Recalling identities (4.11), (4.13), and (4.14), a vector is non-vanishing only when its -quiver is the same as that given by (4.32). Thus, the differential operator (4.32) uniquely selects one vector in the gauge invariant basis while all others vanish, and the expansion coefficient can be solved using an univariate linear equation. Furthermore, the differential operator is normalized to one, 33) and hence, the expansion coefficient can be directly computed by applying differential operator (4.32) on the EYM amplitude, leading to Note that inserts between and relative to the color-ordering, while inserts between n and another graviton . Hence, in the resulting Yang-Mills amplitudes, the legs can never appear in the positions before or after n, and all Yang-Mills amplitudes are in the BCJ basis with legs fixed. An example of evaluating (4.34) has been discussed in (4.22).
For the case, the differential operator for a vector with one pseudo-loop is defined as 1) , 35) with indices following convention (3.58), and the total number of differential operators is In differential operators (4.35), the insertion operator will contribute a derivative relating to momentum . In its quiver, there is only one branch with root , 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 . Therefore, all vectors with two or more pseudo-loops will vanish under (4.35). Furthermore, when applying (4.35) on vectors without a pseudoloop, there could be a non-zero contribution. However, it is not an issue because all coefficients of such vectors have been solved a priori by differential operators (4.32), and they enter into the linear equations as known parameters.
For vectors with one pseudo-loop, there are in general more than one vector that is non-vanishing under a specific differential operator (4.35), as shown in (4.29). Therefore, we need to apply a complete set of differential operators to generate an algebraic system of linear equations and solve the expansion coefficients from this algebraic system. Alternatively, we can also apply the differential operator constructed by rule (4.31), i.e., a special linear combination of differential operators in (4.35). Then, an expansion coefficient can be determined by a univariate linear equation again. Nevertheless, we can compute the coefficient of a vector with one pseudo-loop as where the summation in curly brackets represents a linear combination of differential operators constructed following rule (4.31). Note that the insertion operator inserts between and , so the resulting Yang-Mills amplitudes are no longer in the BCJ basis with legs fixed. The BCJ relations are required in this step to write all Yang-Mills amplitudes into the BCJ basis, while Yang-Mills amplitudes from the contributions of vectors with no pseudo-loops are still in the BCJ basis.
Now, let us proceed to the vectors with p pseudoloops. Using the same argument used with one pseudoloop, by applying the corresponding differential operator, all vectors with or more pseudo-loops will vanish. For different vectors with p pseudo-loops, a linear com-Chinese Physics C Vol. 44, No. 12 (2020) 123104 1) As mentioned, using the simple rule (4.19) we might need to solve algebraic systems of linear equations. While using a more complicated combination of differential operators as (4.31), the algebraic system is decoupled to univariate linear equations.

123104-18
(p − 1) bination of differential operators constructed using rule (4.31) 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 applied to vectors with or fewer pseudo-loops. To solve these linear equations, all coefficients of vectors with or fewer pseudo-loops should be solved a priori, and these linear equations will be entered as known parameters. This inspires us to solve linear equations order by order from vectors with to pseudo-loops. The differential operators relating to vectors with p pseudo-loops in the gauge invariant basis are given as 37) 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 inserts between legs and , and we need to rewrite the resulting Yang-Mills amplitudes into the BCJ basis using the BCJ relations.
The algorithm for evaluation of expansion coefficients can be summarized as follows, --Start of Algorithm --(m + n − 3) m STEP 0: Apply differential operators (4.32) to EYM expansion formula (3.59) to generate linear equations and solve the expansion coefficients from these equations 1) . The result is directly given by (4.34). STEP 1: Substitute solutions of Step-0 back in equation (3.59); then, apply differential operators (4.35) to the resulting formula to generate linear equations. Solve the expansion coefficients from these equations and rewrite the Yang-Mills amplitude into the BCJ basis using the BCJ relations.

STEP p:
Substitute the solutions of all previous steps back into equation (3.59); then, apply differential operators (4.37) on the resulting formula to generate linear equations. Solve the expansion coefficients from these equations and rewrite the Yang-Mills amplitude into the BCJ basis using the BCJ relations.

STEP
: Repeat the previous step but with differential operators. 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 lists 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, and , we see that the latter is approximately 44 times larger than the former when the number of gluons increases by ten. While comparing the total number of equations for a fixed n, for example, and , we see that the latter is approximately 85902 times larger than the former when the number of gravitons increases by four. Hence, the size of the algebraic system is significantly controlled by m. One can also note that the number of equations decreases rapidly while moving to the next step in the algorithm. A large amount of equations are solved in Step-0, where the expansion coefficients are explicitly defined by applying the differential operators to the EYM amplitudes. Therefore, in some sense, it is trivial. For step in the algorithm, the number of equations decreases significantly compared with the previous step; however, non-trivial contributions from previous steps and BCJ relations Chinese Physics C Vol. 44, No. 12 (2020) 123104 A EYM n,m 1) In fact, solving equations is not necessary in this step. The expansion coefficients have been uniquely determined by (4.34), and the remaining thing to do is to explicitly work out the differential operators on according to (4.34).

123104-19
would make the results complicated. Nevertheless, in each step, the linear equation system is decoupled, and an expansion coefficient is trivially solved via a univariate linear equation.

Demonstration of EYM amplitude expansion in gauge invariant vector space
To demonstrate the EYM amplitude expansion in the 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 the BCJ basis has been discussed in [33], however, here, it receives a more systematic analysis in the language of gauge invariant vector space, while expansion of the EYM amplitude with four gravitons to Yang-Mills amplitudes in the BCJ basis is a new result.

Expansion of EYM amplitude with one and two gravitons
Let us start with . This amplitude lives in the gauge invariant vector space , and the dimension of this space is , according to (3.26). Hence, can be expanded in a complete set of the gauge invariant basis with gauge invariant vectors, as The expansion coefficient according to (4.34) is calculated as where the graviton is transformed to a gluon and inserted between by . Hence, In comparison with the result in [33], we can reformulate the above result as A EYM n,1 (1, . . . , n; h Let us proceed to . The dimension of the gauge invariant vector space is . The vectors in the 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 a pseudo-loop. Following the algorithm, Step-0 is to compute the coefficients of the expansion basis with no pseudo-loops, i.e., , , and , using equation (4.34). Applying differential operators and to respectively leads to 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 -4  0  0  0  0  0  0  105  105  105   Total  1515  15370  67075  326479  5262844  35601459  130142105  2892553980  27425768105 Chinese Physics C Vol. 44, No. 12 (2020) 123104 123104-20 . . . , a 1 , h σ 1 , h σ 2 , a 1 + 1, . . . , n), a 1 = a 2 , (5.5) is a permutation of , and the summation is over all elements of . In Step-1, we substitute the above solutions back to the expansion formula and get and there is only one unknown variable, . Applying the differential operator to both sides of the above equation, in the RHS, the non-vanishing contribution comes from vectors and , and according to (4.16) and (4.17), we get In the LHS, we get Then, we arrive at The Yang-Mills amplitudes in the second term are already in the BCJ basis, with legs fixed, while those in the first term are not. Therefore, we need to rewrite the first term in the BCJ basis, as . Combining the above results, we finally obtain

(5.11)
Summing over the total 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 follows: 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 the EYM amplitude with two gravitons in a rather compact form:

Expansion of EYM amplitude with three gravitons
The EYM amplitude lives in the gauge invariant vector space . Because , it is considered to be expanded into terms. Among these gauge invariant vectors, there are terms containing real loops that should be excluded. We therefore need to compute expansion coefficients. The expansion basis and their quiver representations are shown in Fig. 3. Following the algorithm, in Step-0, we consider the gauge invariant vectors with no pseudo-loops using equation (4.34).
Applying the differential operators to the expansion formula of , we immediately obtain where is a permutation of , and the summation is over all elements of . Applying to , we get where is a splitting of . Applying to , we get with . As mentioned, after summing over all the results produced in Step-0, we get a compact expression, Recalling the compact expression (5.12) for the EYM amplitude with two gravitons, we confirm that is it always possible for the total contribution of Step-0 to be written in a compact form. Then, we proceed to Step-1 and compute the expansion coefficients for vectors with one pseudo-loop. After substituting the solutions in Step-0 back into the expansion equation, we get where the first summation runs over all possible splittings of , while the second summation runs over not only all splittings but also all possible values of . 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 are no mixing contributions between different pseudo-loops of the first type. For example, applying the differential operator to equation (5.14), the only surviving vector with one pseudo-loop is . However, vectors with no pseudo-loops would contribute, and from our previous general argument, we can determine the non-vanishing vectors to be and . Hence, we have , (5.15) where the relation (4.16) has been used. Working this out explicitly, we have (5. 16) Terms in the first and second lines are similar to those given in (5.9); hence, we can use the result (5.11) here and immediately work out the summation as The other two terms with can be obtained by permutation of the above result.
We now proceed to the vectors with . As discussed, a defined differential operator (4.35) 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, to compute the coefficient of vector we should choose the differential operator . If instead we choose the other differential operator , both vectors and would be non-vanishing, and their contributions will mix together. Hence, we apply to equation (5.14) and compute the coefficient as follows: ficients of vectors with one pseudo-loop can be computed. Summing up the above results, we get the complete expansion of , which is consistent with the results given in [33].

Expansion of EYM amplitude with four gravitons
The EYM amplitude lives in the gauge invariant vector space , and it can be expanded as a linear combination of vectors. All vectors in the gauge invariant basis and their quivers are shown in Fig. 4. Among them, there are in total vectors with real loops, which should be excluded. We can compute the expansion coefficients of the remaining vectors using the algorithm. Again, in Step-0, we compute the coefficients of vectors with no pseudo-loops using equation (4.34). We shall not write down the explicit coefficient for each basis; rather, we present their summation in a compact expression as 1) , Then, let us proceed with Step-1, to compute the 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.19), we assign each of them to a respective differential operator, and we represent differential operators in quiver representation as where without ambiguity, we have ignored the dashed line corresponding to , which is always linked to the ending point of the cyan line. The first two quivers of differential operators are consistent with rules (4.18), and they are sufficient to distinguish the corresponding vectors uniquely. For the third and fourth quivers of differential operators, noting the choice of direction of the 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. The sixth quiver of differential operators 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 combinations of differential operators as in equation (4.31), in the current simple example, we have an alternative way of solving equations. By first 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 wish 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 those 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 the following contributions in Step-0, For instance, using the differential operator , we can compute the coefficient of as ) .

A EYM
n,4 Applying the differential operator on produces Yang-Mills amplitudes and using the BCJ relations, they can be rewritten into the BCJ basis. For the second topology, the corresponding differential operator picks up the following contributions in Step-0, For instance, using the differential operator , we can compute the coefficient of as ) .

A EYM
n,4 Applying the differential operator on produces Yang-Mills amplitudes  For the third topology, the corresponding differential operator picks up the following contributions in Step-0, For instance, using the differential operator , we can compute the coefficient of as ) .

(5.20)
A EYM n,4 Applying the differential operator on produces Yang-Mills amplitudes is the contribution from expansion in Step-0, and is the contribution from expansion in Step-1. Applying the differential operator on produces Yang-Mills amplitudes

Conclusion
There are already quite a few well-formulated results for expansion of EYM amplitudes to Yang-Mills amplitudes in the KK basis; however, a compact expression or even a recursive formula for expansion to Yang-Mills amplitudes in the BCJ basis is still required. The latter expansion is generally much more complicated than conventionally expected. In the KK basis, the expansion coefficients of Yang-Mills amplitudes are only polynomi-als 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 why we considered using differential operators to determine expansion coefficients in [33]. This paper is motivated by the problem of expanding EYM amplitudes to Yang-Mills amplitudes in the BCJ basis using 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 the BCJ basis but rather to a basis in gauge invariant vector space, as schematically shown in equation (4.1). After determining the expansion coefficients, we transform Yang-Mills amplitudes to the BCJ basis using the BCJ relations. Expanding EYM amplitude in a manifest gauge invariant form for both the expansion basis and its coefficients is a very interesting point of view, and differential operators can be naturally introduced into the problem. This contributes to our major results.
The first major part of this paper is devoted to the construction of a gauge invariant basis and the corresponding differential operators. A systematic algorithm is built upon the properties of applying differential operators on different bases. To construct a complete set of manifestly gauge invariant polynomials as the expansion basis, we start from the most general vector space with , where all possible polynomials of Lorentz contractions among polarizations and momenta live in this space, obeying some additional conditions. Then, we define a linear mapping , which is a realization of gauge invariant condition for polarization. By taking the interaction of kernels of all possible , we construct the gauge invariant sub-space from , 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 , which indicates the number of gauge invariant vectors to which an EYM amplitude would be expanded. We also find that the gauge invariant vectors can be realized by linear combinations of multiplications of fundamental f-terms. These results finally help us construct a linearly independent and complete basis combinatorially for EYM amplitude expansion.
After clarifying the structure of the 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 the gauge invariant basis is non-vanishing, while all others vanish. To do so, we start by analyzing the structures of the gauge invariant basis and finding the quiver representation for them. With the help of the quiver representation, we summarize all possible components appearing in the 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 the 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 the gauge invariant basis, which are all consistent with CHY formalism numerically.
Although the algorithm for expanding the tree-level single-trace EYM amplitude to the Yang-Mills amplitudes in the BCJ basis has been thoroughly provided in this paper, it still inspires further work for the future. First, the expansion coefficients of the BCJ basis demand 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 they can help to determine the general picture. Second, in this paper, we only deal with single-trace EYM amplitudes, while discussions can be generalized to multi-trace EYM amplitudes using the trace operator . We think this generalization should be straightforward. Third, in this paper, we focus on EYM amplitudes, so the parameters of vector space are constrained to . However, the case is also very interesting in physics because Yang-Mills amplitudes live in this space. Another interesting example is the deformed Yang-Mills theory with the term [48,49]. Although the dimension of still holds for , the explicit form of vectors in the gauge invariant basis should be reconsidered because we are not able to trivially exclude momentum in all expressions by momentum conservation. Furthermore, for the Yang-Mills amplitude, an additional constraint should be applied to the vector space, i.e., there should be at least one contraction, and we can denote the vector space by . The new vector space can help us to understand the implication of gauge invariance in Yang-Mills amplitudes more deeply, along the lines of previous studies in [8,37,38]. It is also a curious problem to write Yang-Mills amplitudes in a manifestly gauge invariant form. Perhaps, it can also help us to understand more about the Pfaffian in the integrand of the CHY formula and provide a new point of view for the BCJ relations.
We are grateful to Kang Zhou, ZhongJie Huang, and Yiwen Lin for discussions about this work. Xiao-Di Li would like to thank Yi-Jian Du for his enlightening discussions and kind hospitality at Wuhan University.
From the above equation, we can solve and substitute it back into . After reorganization of terms, we get Because of the appearance of , it is easy to see that . Hence, and . This shows that if , there is also .   Result (A.10) tells us that all and belong to , and by induction, they belong to . For the remaining terms in (A.7), 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 to the RHS of the above equation, we can reorganize to be ) . (A.12) Using result (A.9), we get Thus, belongs to . Therefore, finally, we have proven that is valid in any vector space , and proposition 1 is proven.
Proof of proposition 2: We wish to prove the following distribution formula of linear maps , To do so, it is sufficient to show To show (A.15), we note that any vector can always be written as Thus, we can check that Referring to proposition 1 (A.1), the above result shows that . Hence, (A.15) is derived.
where has been eliminated using momentum conservation. Now, we impose the condition that . Imposing , we obtain the equation T mi ∈ V n,m−2 ϵ i · T ′ mi , T ′′ mi ∈ V n,m−1 Because , , by assumption, they satisfy (A.16). Now, we consider the remaining terms in (A.7); after reorganization of the terms, we get   . Thus, we have successfully separated into two parts satisfying (A.16) in general vector space , and proposition 2 is proven.
For completeness, let us present the proof of (3.23) and (3.24) as follows: