Expansion of Einstein-Yang-Mills Amplitude

In this paper, we provide a thorough study on the expansion of single trace Einstein-Yang-Mills amplitudes into linear combination of color-ordered Yang-Mills amplitudes, from various different perspectives. Using the gauge invariance principle, we propose a recursive construction, where EYM amplitude with any number of gravitons could be expanded into EYM amplitudes with less number of gravitons. Through this construction, we can write down the complete expansion of EYM amplitude in the basis of color-ordered Yang-Mills amplitudes. As a byproduct, we are able to write down the polynomial form of BCJ numerator, i.e., numerators satisfying the color-kinematic duality, for Yang-Mills amplitude. After the discussion of gauge invariance, we move to the BCFW on-shell recursion relation and discuss how the expansion can be understood from the on-shell picture. Finally, we show how to interpret the expansion from the aspect of KLT relation and the way of evaluating the expansion coefficients efficiently.


Introduction
A fairly non-trivial relation between Einstein-Yang-Mills (EYM) amplitude and pure Yang-Mills amplitudes was proposed in [1] recently, where the amplitude of n gluons coupled to a single graviton is expanded as linear sum of (n + 1)-point pure gluon amplitudes in an elegant formulation, which is different from the earlier proposed relations that express n-gluon m-graviton amplitudes by (n + 2m)-gluon amplitudes [2][3][4][5][6].
It is now widely known that, the non-trivial relations among amplitudes are important both in the practical evaluation and the analytical study, while the U (1)-relation, Kleiss-Kuijf(KK)relation [7] and especially the Bern-Carrasco-Johansson(BCJ) relations [8] among amplitudes of the same field theory have received considerable investigations in the past few years, and inspired the color-kinematics duality for gravity and Yang-Mills amplitudes [9,10]. As an analogous scenario, where amplitudes of two originally seemingly unrelated theories take part in novel identity, recall that the famous Kawai-Lewellen-Tye (KLT) relation [11] was proposed quite a long time ago, which formulates a closed string amplitude as products of two open string amplitudes, and in the field theory limit it expands a pure gravity amplitude as bi-linear sum of Yang-Mills amplitudes. The newly proposed linear EYM amplitude relation was also inspired by the study of string theory, where monodromy relations for mixed closed-open string amplitudes, previously been applied to the study of BCJ relations [3,12,13], has been considered.
Because of its compact and simple nature, a substantial research interests has been drawn to the study of EYM amplitude relations and to its generalizations [14][15][16][17][18][19] 1 . In particular most of the discussions are based on the Cachazo-He-Yuan (CHY) formulation [21][22][23][24][25], by genuinely reformulating the CHY-integrand in an appropriate form. Notably, explicit expressions for EYM amplitude relations with arbitrary number of gluons coupled to up to three gravitons were provided in [14]. The technique for reformulating the CHYintegrands in these papers developed into a systematic explanation in [26], and it is revealed therein that the cross-ratio identity and other off-shell identities of integrands [27,28] are crucial tools for deforming CHYintegrands into alternative forms corresponding to different field theories. These powerful tools benefit from the integration rule method [29][30][31][32] developed for the purpose of evaluating CHY-integrand without referring to the scattering equations. The idea of integration rule and cross-ratio identity method was to decompose arbitrary CHY-integrand using cross-ratio identities into those corresponding to cubic-scalar Feynman diagrams dressed with kinematic factors. By carefully organizing terms one can identify the resulting CHY-integrands as amplitudes of certain field theories, hence the amplitude relations, as was done in [14,15]. In fact, there is more about EYM amplitude relations from the perspective of CHYframework. Starting from CHY-integrand of a theory, it is always possible to reformulate it to another form by cross-ratio and other off-shell relations, for instance the Yang-Mills-scalar (YMs) amplitude can be expanded as linear sum of bi-adjoint cubic-scalar amplitudes. We shall discuss this later in this paper.
As it is very often, on-shell technique can prove to be a powerful tool for the purpose of understanding non-trivial amplitude relations within field theory framework. One such example is the on-shell proof of BCJ relations [33,34]. The central idea is to deduce physical identities only from general principles such as locality, unitarity and gauge invariance. This is particularly true with the advent of Britto-Cachazo-Feng-Witten(BCFW) on-shell recursion relation [35,36], which utilizes the first two. In most cases, the BCFW recursion relation computes the amplitude in a way such that only contributions from finite local single poles are summed over, which requires a vanishing behavior in the boundary of BCFW complex parameter plane. This is exactly the case for BCJ relations of Yang-Mills amplitudes. However, for generic situations, the amplitude as a rational function of BCFW parameter z is not vanishing in z → ∞ and the boundary contributions can not be avoid. This is a problem one would meet when applying BCFW recursions to the EYM amplitude relations, and such subtlety complicates the on-shell proof of EYM amplitude relations. The evaluation of boundary contributions is generically a difficult problem, but many methods have been proposed to deal with it. Noteworthily systematic algorithm has also been proposed recently [37][38][39][40][41][42][43] so that at least in principle it is indeed possible to systematically study the EYM amplitude relations using BCFW recursion relations. On the other hand it is also known that very often gauge invariance can become a very handy tool in constraining the specific analytic form of the scattering amplitude. Recent progresses have pushed the gauge invariance principle forward and indicate that, the gauge invariance along with cubic graph expansion are enough to determine the amplitudes [44][45][46][47][48][49]. In a less but still quite challenging situation, we claim that the gauge invariance should uniquely determine the EYM amplitude relations, and from which we can explicitly write down the expansion for EYM amplitude with arbitrary number of gravitons.
As the number of gravitons increases and that of gluons decreases, in the extremal limit we would come to the amplitude with pure gravitons. This is the important problem of expanding gravity amplitude as pure Yang-Mills amplitudes. Furthermore, with the philosophy of decomposing CHY-integrands, the same argument applies to the Yang-Mills amplitudes which would be expanded as pure bi-adjoint cubicscalar amplitudes. This is exactly the cubic-graph expansion of Yang-Mills amplitude which makes the color-kinematics duality manifest [26]. The EYM amplitude relation combined with CHY-integrand, more specifically the Pfaffian expansion, would produce the non-trivial expansion for Yang-Mills amplitude as cubic-scalar graphs, as well as expansion for gravity amplitude as pure Yang-Mills amplitudes and eventually the cubic-scalar graphs. This provides a way of computing the BCJ numerators, which is usually considered to be very difficult [46,[50][51][52][53][54][55][56][57][58][59][60][61][62]. When KLT relation is in action, the EYM amplitude relation can be connected to the BCJ numerator problem. We will learn more about this in later sections.
In this paper we examine the EYM amplitude relations from the perspectives of CHY-formulation, BCFW on-shell recursion, KLT relation, and through the contruction of BCJ numerators. This paper is organized as follows. In §2, we present the general theoretical playground of non-trivial amplitude relations from the CHY-formulation, and explain the expansion of amplitudes as the expansion of Pfaffian of CHYintegrand. In §3, we facilitate the principle of gauge invariance to determine the EYM amplitude relations for gluons coupled to arbitrary number of gravitons. In §4, we generalize the EYM amplitude relations to pure Yang-Mills amplitudes and apply the non-trivial relation to the computation of BCJ numerators. In §5, we provide the on-shell proof of some EYM amplitude relations by BCFW recursion relations. In §6 we study in the language of KLT relations. Conclusion is presented in §7 and some useful backgrounds are summarized in the Appendix.

Amplitude relations from the perspective of CHY-formulation
The non-trivial relation revealed recently between EYM amplitudes and pure Yang-Mills amplitudes [1,14,15] has an intuitive interpretation in the CHY-framework. In fact, the CHY-formulation tells more beyond the EYM amplitudes. In the CHY-formula, it is the so called CHY-integrand I CHY that describes specific field theories. The CHY-integrand is an uniform weight-4 rational function of n complex variables z i for n-point scattering system, i.e., with the 1/z 4 i scaling behavior in the z i → ∞ limit. For almost all known theories, the weight-4 CHY-integrand can be factorized as two weight-2 ingredients, formally written as (2.1) Let us then define two new weight-4 CHY integrands as follows where PT(α) is the Parke-Taylor factor Supposing the two CHY-integrands I CHY L , I CHY R also describe certain physical meaningful field theories and produce the corresponding color-ordered amplitudes A L (α), A R (β) after CHY-evaluation, then by CHYconstruction [21][22][23] we could arrive at the following generalized KLT relation, where A is the amplitude of specific field theory determined by the theories of A L , A R , while S n denotes permutations on n elements and S[σ| σ] is some kinematic kernel. The summation is over S n−3 permutations of sets {2, . . . , n − 2}, depending on our choice of legs k 1 , k n−1 , k n being fixed. The expression (2.4) denotes a general expansion for the original amplitude A defined by CHYintegrand (2.1). If for a specific σ ordering, we sum over all S n−3 permutations of σ and define the result as then the original amplitude can be expressed as where C( σ) serves as the expansion coefficients. Similarly, if for a specific σ ordering we sum over all S n−3 permutations of σ and define the summation as then the original amplitude can be expanded as The expressions (2.6) and (2.8) have provided two different expansions of the original theory. There are several general remarks regarding the expansion in above, • Firstly, the expansion is into a chosen (n−3)! BCJ basis, and the corresponding expansion coefficients C( σ) and C(σ) would also be unique. However, as we will discuss soon, sometimes it is better to expand the original amplitude into the (n − 2)! KK basis. Because of the BCJ relations among color-ordered partial amplitudes, the expansion coefficients in the KK basis will not be unique and depend on the generalized gauge choice in the BCJ sense.
• Secondly, with the amplitude expansion formula in hand, the next is to compute the expansion coefficients. For this purpose, there are several approaches. The first approach is to use the definitions (2.5) and (2.7) directly. However, in general it is very hard to evaluate the summation for generic n-point situation, and only in certain special case a direct evaluation is possible, which we shall explain later. The second approach seeds back to the expression (2.1), and the major idea is to expand the weight-2 ingredients I L or I R into the PT(α) factor of n elements. In fact, this is the approach followed in [14,15]. The expansion can be systematically achieved by successively applying cross-ratio identities to the CHY-integrands, where in each step a gauge choice should be taken in the cross-ratio identity. In general, such expansion leads to a result with (n − 1)! cyclic basis. Then one can use the KK relation to rewrite it into the (n − 2)! KK basis. As already mentioned, the gauge dependence remains in the expansion coefficients at each step, and it would disappear only after using the BCJ relations to rewrite all into (n − 3)! BCJ basis.
Besides the above two direct evaluation methods for expansion coefficients, there are also some indirect ways. For example, one can propose some ansatz for the expansion coefficients, then prove and generalize it by on-shell recursion relations. One can also use some general considerations, for instance the gauge invariance or the soft behavior, to determine the coefficients [47][48][49].
In this paper, we will investigate the expansion from these different views.
• Thirdly, although in most theories, the CHY-integrand is given by products of two weight-2 ingredients as (2.1), for some theories the CHY-integrand is defined by the product of four weight-1 ingredients. So there are various combinations of them to form weight-2 parts. In other words, there are possibilities to have more than two expansions given in (2.6) and (2.8). It would be interesting to survey the consequence of different combinations for these theories.
After above general discussions, now we focus on our major topic in this paper, i.e., the single trace part of EYM theory, whose CHY-integrand is defined as for scattering system of r gluons and s gravitons with r + s = n, and α = {α 1 , . . . , α r }. We can define two new CHY-integrands as where σ, σ ∈ S n . It is easy to tell that the I CHY L is the CHY-integrand of Yang-Mills-scalar(YMs) theory and I CHY R is the CHY-integrand of Yang-Mills theory. Correspondingly, the amplitude A L is the color-ordered YMs amplitude A YMs r,s with r scalars and s gluons, which has two trace structures associated with the two PT-factors, while the amplitude A R is color-ordered Yang-Mills amplitude A YM n . One thing to emphasize is that the scalar carries two groups (one gauge group and one flavor group) and has bi-adjoint scalar-cubic interactions.
An immediate consequence from (2.4) reads The expansion (2.11) is into the BCJ basis with (n − 3)! independent Yang-Mills amplitudes. However, as it will be clear soon, an expansion into (n − 2)! KK basis is more favorable, and we would present it here as with the expansion coefficients (2.14) The expansion coefficient in (2.13) is the desired quantity we want to compute in this paper. As we have discussed in previous paragraph, these coefficients are determined by only one weight-2 ingredient in the CHY-integrand in (2.1). This means that while keeping the same weight-2 ingredient, we have the freedom to change the other weight-2 ingredient. As an implication of such modification, we could work out the expansion for different field theories but with the same expansion coefficients. This freedom could simplify our investigation of expansion coefficients. For example, in the context of EYM amplitude as Yang-Mills amplitudes, we can change the Pf Ψ n in (2.9) as PT n (ρ). The resulting CHY-integrand describes a Yang-Mills-scalar amplitude with r scalars and s gluons, and the weight-2 ingredients are now I L = PT r (α)Pf Ψ s and I R = PT n (ρ). With the same philosophy as in (2.10), we can define two new CHY-integrands as Hence we arrive at the following expansion with the same expansion coefficients as in (2.12). This non-trivial relation expresses any single trace color-ordered amplitude of Yang-Mills-scalar theory as linear combination of color-ordered amplitude of bi-adjoint scalar φ 3 theory. After studying the expansion of single trace part of EYM theory to YM theory, we will briefly discuss the expansion of gravity theory to YM theory. The CHY-integrand of gravity theory is If expanding the reduced Pfaffian by cross-ratio identities, we will get by (2.18). As already pointed out in papers [23,54,57,[63][64][65][66], the coefficients n(1, σ, n) in the expansion (2.19)(hence also the one in the expansion (2.20)) is nothing but the DDM basis for the BCJ numerator of YM amplitude. While in the expansion (2.8), i.e., A = σ∈S n−3 C(σ)A L (n − 1, n, σ, 1), suppose we can rewrite the (n − 3)! BCJ basis into (n − 2)! KK basis, then identifying the resulting formula with the one given by (2.20), and equaling the expansion coefficients of the same KK basis, we will get the BCJ numerator n(1, σ, n) as linear combination of C(σ). Thus here we provided a new way of computing the BCJ numerators via the computation of amplitude expansion (2.8). Although in (2.20) we have taken gravity amplitude as example, the same consideration can be applied to large number of theories, and the BCJ numerators of those theories can also be identified as the expansion coefficients after rewriting (2.8) into KK basis. With the above general framework for studying the non-trivial relations among different theories, we will start our exploration from next section. As we will see, these relations encode many surprising structures and connect many important topics, such as the construction of BCJ numerators mentioned above, the boundary contribution of amplitudes under BCFW on-shell recursion relation, the DDM chain and BCJ relations, etc,.

The gauge invariance determines the amplitude relations
A physical amplitude should be gauge invariant, i.e., vanishes on the condition i → k i . If considering the amplitudes with gravitons and expressing the graviton polarization states as products of two Yang-Mills polarization states ± i ± i , then it should also vanish by setting one of the polarization vector as k i . The gauge invariance is an important property of amplitude, and of course it is also valid in the amplitude relations. As we have mentioned in previous section, there are different approaches to study expansion coefficients. In this section, we will demonstrate how to use the gauge invariance to fix coefficients, which is the same spirit spelled out in [47][48][49].

With single graviton
To motivate our discussion, let us start with the single trace EYM amplitude with one graviton with the known expansion where the summation ¡ is over all shuffles σ ¡ σ, i.e., all permutation sets of σ ∪ σ while preserving the ordering of each σ, σ. The color-ordering of external legs in A YM n+1 has cyclic invariance. However if we conventionally fix the leg 1 in the first position, then every external leg could have a definite position in the color-ordering. In this convention, we can define Y p (also X p in the following paragraph) as the sum of momenta of all the gluons at the left hand side (LHS) of leg p, given the definite color-ordering of color-ordered YM amplitudes.
A clarification of the definition Y p is needed here for the future usage. In the (n + 1)-point pure Yang-Mills amplitude, the gluon legs has two different correspondents in the EYM amplitude, i.e., The gluon legs k i , i = 1, . . . , n are also gluons in EYM amplitude while gluon leg p is originally graviton leg in EYM amplitude. Y p is specifically defined as the sum of momenta in the gluon subset of EYM amplitude at the LHS of p, while we also define another quantity X p as the sum of all momenta at the LHS of p no matter it is in the gluon subset or graviton subset of EYM amplitude. X p , Y p would be different when considering EYM amplitude with more than one gravitons, but in the current discussion they are the same.
Let us now return to the relation (3.1). In the LHS, imposing any gauge condition i → k i for gluon legs would vanish the EYM amplitude, while any Yang-Mills amplitudes in the right hand side(RHS) also vanish under the same gauge condition. For the graviton polarization ± p ± p , setting either p → p would vanish the EYM amplitude. In the RHS, the graviton polarization is distributed in two places: one is in the Yang-Mills amplitude and the other, in the expansion coefficients. The vanish of RHS for the former case is trivial, while for the latter case, it vanishes due to the fundamental BCJ relations This consequence is rather interesting. For the non-trivial relation (3.1) to be true and the gauge invariance be not violated, we eventually end up with BCJ relations. On the other hand, if we assume A EYM n,1 can be expanded as linear combination of Yang-Mills amplitudes in the KK basis for convenience, and the expansion coefficients should be certain sum of p · k i to compensate the extra graviton polarization, then α (p · x p )A YM n+1 (1, α 2 , . . . , α n−1 , α p , n) should be in the BCJ relation form! The lesson learned from the EYM amplitude with one graviton suggests that, while expressing EYM amplitudes as linear combination of Yang-Mills amplitudes, the gauge invariance strongly constraints the possible form of expansion coefficients. This property motivates us to find the expansion of the single trace multi-graviton EYM amplitude with more than one graviton by gauge invariance, i.e., we want coefficients to make the expression to zero under each gauge condition In order to construct the non-trivial relations for EYM amplitudes with generic number of gravitons, we start with the following two ansatz, • Ansatz 1: The single trace EYM amplitude A EYM n,m with m gravitons can always be expanded into EYM amplitudes A EYM n+m−m ,m with m < m gravitons.

(3.5)
These two ansatz come from lower-point known results. The first ansatz is in fact a general statement saying that a recursive construction for EYM amplitude expansion exists. While the second ansatz is presented in an explicit expression which has an obvious BCJ-like relation form. The validation of ansatz 2 in fact can be verified by BCFW recursions. In the expression (3.5), we emphasize again that X h i is defined to be the sum of all momenta in the LHS of leg h i , no matter those legs representing gluons or gravitons in the EYM amplitude.
Bearing in mind that any EYM amplitude expansion relations should follow the above mentioned two ansatz, we are now ready to extend the introductory one graviton example to EYM amplitudes with arbitrary number of gravitons. However, before presenting the general algorithm, let us familiarize ourselves by studying EYM amplitudes with two and three gravitons.

Expressing n-gluon two-graviton EYM amplitudes as Yang-Mills amplitudes
The algorithm for producing general EYM amplitude relations is to expand A EYM n,m as A EYM n+1,m−1 successively until A EYM n+m,0 ≡ A YM m+n . Note that the gravitons are colorless, and it has no color-ordering in EYM amplitudes. But in our construction of EYM amplitude relation by gauge invariance principle, it is necessary to specify a graviton in each step of expansion A EYM n,m → A EYM n+1,m−1 , which in the A EYM n,m amplitude it denotes a graviton but in the A EYM n+1,m−1 amplitude it denotes a gluon, such that we can apply gauge invariance principle with that graviton. Furthermore, it requires us to select one arbitrary graviton to start with. Now let us outline the arguments that lead to the correct expansion of EYM amplitude with two gravitons A EYM n,2 (1, 2, . . . , n; p, q). In the first step, let us specify the graviton h p , and following the Ansatz 1 let us propose the following terms that would contribute to the expansion of A EYM n,2 , (3.6) In fact, these proposed terms (3.6) are reminiscent of the expression (3.1) for expanding the single trace EYM amplitude with one graviton. This is not yet the complete expansion expression for A EYM n,2 , but we will explain soon after how to recover the remaining terms. Let us investigate the gauge invariance of gravitons h p and h q for the proposed terms (3.6). The gauge invariance for h q is obvious since A EYM n+1,1 (· · · ; q) vanishes under q → q gauge condition. However, T 1 is not gauge invariant under p → p due to the expansion coefficients p · Y p , and there are some missing terms in order to produce the complete expansion for A EYM n,2 . Let us proceed to expand the A EYM n+1,1 in T 1 with the known result (3.1), which gives  Note that in the permutation shuffle {2, . . . , n − 1} ¡{p}¡{q}, the position of leg q would be either at the LHS of p or RHS of p. But the leg p denotes a graviton in A EYM n,2 . So the expansion coefficient is q · X q but not q · Y q (remind that X q is the sum of all momenta in the LHS of q, while Y q is the sum of all momenta in the LHS of q excluded the leg p, which means that if p is at the RHS of q, X q = Y q , but if p is at the LHS of q, X q = Y q + p).
From the Ansatz 2 (3.5), we know that in the A EYM n,2 expansion, the correct terms with coefficients ( q · •)( q · •) must be Comparing T 1 (3.7) with the correct result (3.8), it is easy to see that for those terms with p in the LHS of q, Y p = X p so that the corresponding terms in T 1 and (3.8) are the same. While for those terms with p in the RHS of q, we have Y p + q = X p . So in order to reproduce the correct result (3.8), we should add another contribution is exactly equivalent to the correct result (3.8). Remind that Ansatz 2 gives correct answer for contributing terms without ( · ) coefficients for EYM amplitude expansion, and up to now, we have reformulated the correct result as (3.10) which guarantees an easy generalization.
Of course, (3.10) is still not yet the complete expansion for A EYM n,2 , since those terms with coefficient ( p · q ) are still missing. Let us propose that the complete expansion takes the form 11) and the remaining task is to determine T 3 . It can be determined either by gauge condition p → p or by q → q, however the latter is much more convenient since T 1 is already manifestly gauge invariant for leg q. Setting q → q, we have which has a solution Hence we get the non-trivial relation for EYM amplitude with two gravitons as (3.14) If we define the tensor the above EYM amplitude expansion can be reformulated in a more compact form as (3.16) From expression (3.16) we can infer some important features. Firstly, for terms in the first line, leg p denotes a gluon and leg q denotes a graviton, while for terms in the second line, leg q denotes gluon instead of graviton. This difference leads to the difference of expansion coefficient, such that the Y p factor in the first line has been replaced by the factor F q · X q . Or we can say a factor F q is inserted. As we would see shortly after, this is a general pattern for EYM amplitudes involving more gravitons.
Secondly, in the expression (3.16), the gauge invariance for leg q is manifest, while gauge invariance for leg p is not manifest and requires further checking. Although it can be checked directly, we will follow another approach. Note that the whole result should be symmetric under switching p ↔ q. For the terms with kinematic factors ( · k)( · k), this symmetry is manifest since it is given by (3.8). For the terms with kinematic factors ( p · q ), the result is not manifestly symmetric. In order to shown the symmetry, we need to use the generalized BCJ relation. Let us divide the set {2, ..., n − 1} into two ordered subsets α = {a 1 , ..., α m } and β = {β 1 , ..., β t } such that m + t = n − 2, then the general BCJ relation is given by [12,34] where the first summation is over all shuffles, and X β i is the summation of all momenta of legs at the LHS of leg β i . Using (3.17) with β = {q, p} it is easy to see that Hence the symmetry of legs q, p for the terms with factor ( p · q ) is apparent. Since the gauge invariance for leg q is satisfied, by the symmetry, the gauge invariance for leg p is also satisfied.
The above discussion allows a systematical generalization to EYM amplitude with any number of gravitons. Before doing so, let us introduce a new quantity which would be useful in later discussions. Assuming the gravitons have been split into two subsets α, β, where α is the ordered length-m 1 set in the gluon side and β is a length-m 2 set in the graviton side whose ordering is not relevant, we define which is in fact the first line of (3.14). While if α = {q, p}, β = ∅, m 1 = 2, m 2 = 0, and we have (3.22) which reproduces the second line of (3.14).

Expressing n-gluon three-graviton EYM amplitudes as Yang-Mills amplitudes
Now let us explore the details by the EYM amplitude with three gravitons A EYM n,3 (1, . . . , n; p, q, r). Our purpose is to construct a recursive algorithm for EYM amplitude expansion which is manifestly gauge invariant in each step for gravitons (except the initial one), and the terms without ( h i · h j ) factor matches the Ansatz 2 (3.5). In the current case it is Again, the starting point is specifying an arbitrary graviton for expanding A EYM n,3 → A EYM n+1,2 and without lose of generality we choice p. The proposed contributing terms are To match the correct result (3.23), we need to add terms where leg q or r is at the LHS of p. This means that we need to add terms A EYM n+2,1 where leg q or leg r now denotes gluon and its position is at the LHS of leg p. For A EYM n+2,1 terms where leg p, q are gluons but leg r is graviton, the added term should be . . , n − 1} ¡ {q, p}, n; r) .

(3.25)
These terms introduce the missing terms for T [{p}|{q, r}] in order to match the result (3.23), however the gauge invariance for q is still broken. In order to keep gauge invariance for legs q, r at every step, we should further modify (3.25) by adding terms with ( p · q ) coefficients, and the resulting terms should not alter the matching with result (3.23). From experiences gained in the previous subsection, we can propose the following modification which are manifestly gauge invariant for q, r.
Similarly, for A EYM n+2,1 terms where legs p, r are gluons but leg q is graviton, the proposed gauge invariant term should be Emphasize again that the above proposals are based on the gauge invariant principle, the Ansatz 1 and the Ansatz 2.
Notice that in (3.26) and (3.27), we have ( q · Y q ) or ( r · Y r ) instead of ( q,r · X q,r ). So in order to arrive at a complete matching with result (3.23), we should further add A EYM n+3,0 terms where all p, q, r are gluon legs. For G[{q, p}|{r}], the Ansatz 1 guides us to propose additional terms as However, these terms are not gauge invariant for leg r, and according to gauge invariance principle we need to modify (3.28) as which reproduces the correct result (3.23) yet is gauge invariant manifestly. The ( p · F q · r ) coefficient in the second line of (3.29) is so we can see clearly that, the second line of (3.29) only introduces terms with ( h i · h j ) factor which will not contribute to the (3.23) terms.
Similarly, for G[{r, p}|{q}], we need to add the following gauge invariant terms    . The correctness of terms without ( h i · h j ) is guaranteed by construction, while the terms with ( h i · h j ) factor are determined by gauge invariance in each step. The gauge invariance for q, r is then manifest at each term, except for the leg p. It is also easy to see that, in each step when leg h i in the amplitude denotes a gluon while in the previous step it denotes a graviton, the corresponding gauge invariant term is just to insert a F h i into the kinematic factor in an appropriate position. It corresponds to replacing k µ (3.20), we can define a new quantity Before presenting the algorithm for general EYM amplitude relations, we give a remark on the gauge invariance of p. It is not apparent, but one can show the full S 3 symmetry among three gravitons after using various BCJ relations. Hence the gauge invariance of leg p should be indicated by the symmetry.

A constructive algorithm for producing general EYM amplitude relations
The basic idea of constructive algorithm for producing general EYM amplitude relations A EYM n,m is to write down the contributing terms A EYM n+1,m−1 , A EYM n+2,m−2 , . . . , A EYM n+m,0 successively, and the explicit expression corresponding to A EYM n+m ,m−m relies on A EYM n+m −1,m−m +1 recursively. Briefly speaking, provided we have written down the contribution of A EYM n−m 2 ,m 2 , where the gravitons in this EYM amplitude are labeled as β = {β 1 , . . . , β m 2 }. Then by specifying a graviton, say h β i , Step m: For each G in the previous step, specify one graviton and record the contribution   [47][48][49], for n-point Yang-Mills amplitudes, manifest gauge invariance for (n − 1) points will be enough to guarantee the correctness of the result, so the gauge invariance of the n-th point. We believe the same conclusion can be made for EYM theory by similar argument. If we buy this argument, result (3.39) must be the right expression. As a demonstration, let us briefly present the non-trivial relations for EYM amplitude with four gravitons A EYM n,4 (1, . . . , n; h 1 , h 2 , h 3 , h 4 ). The contributions in each step are abbreviated as follows, .

Expanding to pure Yang-Mills amplitudes: ordered splitting formula
The recursive construction given in (3.39) is easy to implement and one can eventually get an expansion with pure Yang-Mills amplitudes. In this subsection, we will present the related discussion.
To familiarize ourselves with this problem, let us start with some examples. The first example is the one with two gravitons. After substituting (3.1) into the first term of (3.16), we get A EYM n,2 (1, 2, . . . , n; p, q) = where it is crucial to use X q instead of Y q in the first term of the expansion, since to the leg q, leg p is actually a gluon. Although the expression (  we need to specify a graviton leg which would be the gluon leg in A EYM n+2,1 , and our choice is leg q. Secondly, we have defined a new notation Z h i . To define Z h i , we shall introduce a new concept, i.e., the ordered splitting of m elements. To define the ordered splitting, we must first define an ordering of m elements, for example, h 1 ≺ h 2 ≺ · · · ≺ h m (we will call it ordered gauge ). Once the ordered gauge is fixed, the ordered splitting is then defined by the following ordered set of subsets {α 1 , . . . , α t } satisfying following conditions, • Denoting R α i as the last element of the ordered subset α i (or named the pivot), then R α 1 ≺ R α 2 ≺ · · · ≺ R αt according to the ordered gauge(it defines the ordering of subset α i in the set {a 1 , . . . , a t }), • In each subset, all other elements must be larger than R α i according to the ordered gauge. However, there is no ordering requirement for all other elements.
To better understand the definition of ordered splitting, we take the set {p, q, r} with ordered gauge p ≺ q ≺ r as an example to write down all ordered splitting, Now let us define the notation Z h i . It is easy to notice that, the six lines in (3.43) are one-to-one mapped to the six ordered splitting of {p, q, r} with ordered gauge p ≺ q ≺ r. The Z h i is the sum of momenta of legs satisfying the following two conditions: (1) legs at the LHS of the leg h i in the colorordered Yang-Mills amplitudes, (2) legs at the LHS of the label chain defined by the ordered splitting. The label chain for a given ordered splitting is the ordered set {1, 2, . . . , n − 1, α 1 , . . . , α t , n}. For instance, for the ordered splitting {{p}, {q}, {r}} in the first line of (3.43), the label chain is {1, 2, . . . , n − 1, p, q, r, n}, and for {{p}, {r, q}} in the second line of (3.43), the label chain is {1, 2, . . . , n − 1, p, r, q, n}.
With the understanding of Z h i , it is easy to see that all Y h i appearing in (3.43) is equal to Z h i , so we can rewrite (3.43) as A EYM n,3 (1, . . . , n; p, q, r) We have numerically checked the above relation, by comparing with A EYM n,3 directly evaluated with the CHY definition and found agreements in the lower-point examples up to A EYM 3,3 . In addition, when expanding the amplitude A EYM n,3 (1, . . . , n; p, q, r) into terms of pure Yang-Mills ones by (3.44), (3.42) and (3.1) and considering BCJ relations, we obtain the same result proposed in [14].
With the above result (3.44), it is ready to outline the rule for generalizing (3.43) for arbitrary number of gravitons, • Decide an ordered gauge a priori, and write down all possible ordered splitting.
• For each ordered splitting {α 1 , α 2 , . . . , α t }, write down a factor ( e 1 · F e 2 · · · F e |α i |−1 · F e |α i | · Z e |α i | ) for each subset α i = {e |α i | , . . . , e 2 , e 1 }, and product the factors for all α i s in the ordered splitting. This is the desired coefficients for the color-ordered amplitudes with color-ordering defined by the corresponding ordered splitting.
• Sum over the results from all possible ordered splitting, and we get the desired expansion of EYM amplitudes into pure Yang-Mills amplitudes.
A final remark says that, because of the choice of ordered gauge a priori, in the expansion the gauge invariance of gravitons is not as obvious as the one given by the recursive expansion in the previous subsection.

Expanding to pure Yang-Mills amplitudes: KK basis formula
The formula (3.44) and its generalizations provide an expansion of EYM amplitudes into Yang-Mills amplitudes in the framework of ordered splitting. However, to get the explicit expression for BCJ numerators, we need an expansion based on KK basis. From the framework of ordered splitting to the KK basis is somehow straightforward, and the only problem is to identify all the corresponding ordered splitting that can produce the desired KK basis. More explicitly, we need to reconstruct the ordered splitting from a given color-ordering in KK basis. It is easy to propose such algorithm, and we would like to demonstrate with a six-graviton example.
Assuming the ordered gauge is • Start with the lowest element in the ordered gauge, here h 1 , and write down all possible subsets α i ⊂ {h 1 , . . . , h 6 } which contains h 1 as its last element respecting the given color-ordering in KK basis. Since h 5 , h 4 are at the LHS of h 1 , we can write down four subsets, • Repeat until we get the complete ordered splitting.
For our current example, we can write the recursive construction starting from α h 1 ,4 as follows, This gives two ordered splitting. After generating all possible ordered splitting regarding the given color-ordering in KK basis, we can readily write down its coefficient in KK basis as the sum of those from the ordered splitting. However, we remark that, the definition of Z h i relies on the ordered splitting, so the explicit expression of Z h i in different ordered splitting is different. Thus we need to be careful when summing over results of all possible ordered splitting.
Finally, let us present an complete example with three gravitons to demonstrate the algorithm, which is shown as follows, • Choose the ordered gauge as p ≺ q ≺ r.
As we have emphasized, Z q has different meaning in the first and third terms. It will contain k r in the third term, but not in the first term. • The above example clearly shows how to reproduce the EYM amplitude expansion from the ordered splitting formula to the KK basis. Although for EYM amplitudes with many gravitons the procedure would be quite involved, but with the help of computer, it would not be a serious problem. And the expansion in KK basis is ideal for the study of BCJ numerators, which we will mention in the next section.

The BCJ numerator of Yang-Mills theory
In this section, we will apply the story of EYM theory to pure Yang-Mills theory, to provide an algorithmic construction of BCJ numerators for Yang-Mills amplitudes. Let us start from some necessary backgrounds. In paper [67], an expansion of Pfaffian is introduced as, where the sum is over permutations S n , and I, J, . . . , K are closed cycles from splitting of the permutation. The factor Ψ I is defined as for I containing more than one z i s, and for I containing only one z i , which is explicitly the diagonal term of C matrix. Although the definition (4.1) is proposed for the full 2n × 2n Pfaffian matrix, it is also valid for the sub-matrix of Pfaffian. In the case of Yang-Mills theory, we are interested in the reduced Pfaffian, where the λ and ν-th rows and columns in Ψ have been removed. The prime in the indicates that the sum is taken over all p ∈ S n permutation such that ν is changed into λ, which we call constraint permutation. The constraint permutation p has been decomposed to closed cycles I, J, . . . , K. In this paper, we take 1, n as the gauge choice, so the constraint permutation is the closed cycle I = (1i 2 i 3 · · · i |I| n), and With above gauge choice, we can decompose the constraint permutation sum in (4.4) as three summations. The first summation is the summation of splitting {2, 3, . . . , n − 1} into a gluon subset A and a graviton subset B, while both subsets could be empty. The second summation is the permutation over gluon subset A. The third summation is the permutation over graviton subset B. This gives 3 where in the second line we have used (4.1). 3 Naively, there will be a sign factor (−) 1 2 n B (n B +1) if using the result (4.1). However, comparing with our explicit example and the argument using the BCFW on-shell recursion relation, it seems this sign factor should not be there.
If we multiply (4.6) by another copy of reduced Pfaffian, then the expression in the curly bracket in the third line of (4.6) can be identified as CHY-integrand of single trace EYM amplitude with gluon legs {1, A, n} and graviton legs B. Thus the above relation between reduced Pfaffian and Pfaffians for fewer points, exactly establishes the same relation between Einstein gravity and EYM/Yang-Mills amplitudes. As already pointed out in §2, the relation between gravity and YM amplitudes shares the same kinematic coefficients for relations between Yang-Mills amplitude and cubic-scalar amplitudes. Since we have already worked out the EYM amplitude relations in the previous section, we can use the known results for EYM amplitude to directly write down the basis of BCJ numerators for Yang-Mills theory, without going through again the expansion of Yang-Mills amplitudes into cubic-scalar amplitudes. Now we present the algorithm for finding the BCJ numerator of one KK basis with the color-ordered A YM (1, h 2 , h 3 • With the knowledge of ρ(A), we generate an ordered subset O(B) = {h 2 , . . . , h n−1 }/ρ(A). Then we choose the standard ordered gauge 2 ≺ 3 ≺ · · · ≺ n − 1, and reconstruct the ordered splitting as the story in (3.46). Each ordered splitting contributes a factor C i .
• Collecting results for all possible ordered splitting, we get the coefficient, As a demonstration of above algorithm, let us compute the KK basis of BCJ numerators for four-point Yang-Mills amplitude. For the numerator n (1|23|4) , we can generate the following terms, Summing over all the results, we get Similarly, for n (1|32|4) we can reconstruct the following terms, Summing over all results, we get In paper [26], the BCJ numerators for four-point Yang-Mills are computed by directly manipulating the Pf Ψ using cross-ratio identities, and the result is given by where t is the gauge choice in cross-ratio identities. If taking t = 1/2, we get the relabeling symmetric expression between 2 ↔ 3. While if taking t = 0, we get which has a perfect agreement with results (4.8) and (4.9) up to a factor −2.
Although only four-point example is explicitly displayed, this algorithm can directly be applied to obtain BCJ numerators for tree level Yang-Mills amplitudes with more than four points.

Inspecting the amplitude relations through BCFW recursions
Using gauge invariance, we have determined the non-trivial relations between any EYM amplitudes and Yang-Mills amplitudes, with the formula as shown in (3.39) for generel EYM amplitudes. However, it would be interesting to inspect it again by BCFW recursions. As mentioned before, the relation between color-ordered EYM amplitudes and color-ordered Yang-Mills amplitudes is the same as that between colorordered Yang-Mills-scalar amplitudes and pure scalar amplitudes. So for simplicity, we would use BCFW recursion to inspect the following two non-trivial relations, (1, {2, . . . , n − 1} ¡ {q, p}, n) .

(5.2)
For more complicated relations where A YMs n,m with m > 2, of course we can also investigate them by BCFW recursions. While m ≥ 2, we need to consider both contribution from finite poles and boundary contribution, and although for m > 2 cases the analysis would be more involved, the techniques for computing boundary contribution are no more than the m = 2 case. Hence their analysis is similar to the relations of A YMs n,2 .
In fact, verify the relation (5.1) by BCFW recursion is trivial. We only need to apply BCFW on both sides of (5.1), while a k 1 , k n shifting is sufficient to detect all the contributions for A YMs n,1 and A φ 3 n+1 . By using the relation (5.1) for A YMs n ,1 with n < n, we can prove it inductively. While the starting point of induction, i.e., the relation for A YMs 2,1 , can be verified explicitly. For the two-scalar one-gluon amplitude A YMs 2,1 (1, 2; p), we explicitly have is the three-point scalar amplitude with cubic vertex defined as f a 1 a p a 2 . We have absorbed the factor i √ 2 into p . This operation just changes the normalization factor and will not affect the following discussions.
We will not repeat the trivial proof of (5.1) here, but jump to the more typical one (5.2). While applying BCFW deformation, for (5.1) only residues of finite poles will contribute, but for (5.2) the boundary contribution is un-avoidable. So we need to compare both sides of (5.2) with the finite pole contributions as well as boundary contribution, and the later will be computed by analyzing Feynman diagrams [68][69][70]. Now let us study the relation for A YMs n,2 (1, . . . , n; p, q) with two gluons by BCFW recursion. We shift the momenta of two scalars k 1 and k n , i.e., where q satisfies k 1 · q = k n · q = q 2 = 0. We shall emphasize that the on-shell condition of scalar k n is not really used in the following proof, thus this relation is also valid for amplitudes with k n off-shell. As the standard BCFW recursion arguments, under this deformation, amplitude A as a rational function of z can be written as A =

Finite Poles
Res A(z) z + Boundary terms. (5.5) To prove the relation (5.2), we should confirm, (1) the sums over finite poles for the LHS and RHS of (5.2) match with each other, (2) the boundary terms for the LHS and RHS of (5.2) match with each other.

Contributions of finite poles
Let us first treat the contributions from finite poles. According to BCFW recursion, the contribution in the LHS of (5.2) can be written as where A(α|β) is short for A L (α, P α ) i P 2 α A R (−P α , β), and P α is the sum of all momenta in set α, which is the propagator in between the left and right sub-amplitudes. The contribution from the finite poles in the RHS of (5.2) can be given by the sum of F 1 , F 2 defined as  and where remind again Y p is defined as the sum of the momenta at the LHS of p in its corresponding subamplitude which has scalar origin in A YMs n,2 . While p is the the right sub-amplitude A R , the intermediate propagator P α appears as the first leg of A R which is at the LHS of p, and we write it explicitly out of Y p to avoid ambiguity.
Let us now compare the result in (5.6) and F 1 , F 2 . From inductive assumption, we know A YMs n,1 for any n and the lower-point amplitude A YMs n ,2 with n < n satisfy (5.1), (5.2). So we have, 1. The sum of A L in (5.7) is equal to the A L in the first term of (5.6) by relations of YMs amplitude with one gluon.
2. In the (5.13), , while for the latter factor, Noting that this term plus (5.10), and that p · P 1...i + p · q = p · P 1...i,q , we get which is equal to the second term of (5.6) by relations of YMs amplitude with one gluon for A R .

The
which equals to the fourth term of (5.6) by relations of lower-point YMs amplitude with two gluons.
Hence, all contributions of finite poles in the LHS and RHS of (5.2) under (k 1 , k n )-shifting match with each other. Figure 1. Boundary of the left hand side of relation for amplitudes with more than two scalars and two gluons.

The boundary contributions
Now let us discuss the boundary contributions. We would consider the situations with n > 2 and n = 2 separately, since there are more subtleties in the latter situation.
The case with n > 2: in this case, the boundary term of the BCFW recursion under (k 1 , k n )-deformation comes from the diagram as shown in Fig.(1.a). In the LHS of (5.2), it is given by f a n a 1 e i s 1n A YMs n−1,2 (2, . . . , n − 1, P e 1,n ; p, q) , (5.14) where A YMs n−1,2 (P e 1,n , 2, . . . , n − 1; p, q) denotes sum of all possible diagrams with one off-shell scalar line P e 1,n , (n − 2) on-shell scalars k 2 , . . . , k n−1 as well as two on-shell gluons p, q. In the RHS of (5.2), it is given by we re-express the boundary term in the RHS of (5.2) as ¡ ( p · Y p )f a n a 1 e i s 1n A YMs n,1 (2, {3, . . . , n − 1} ¡ {p}, P e 1,n ; q) Remind that our proof of (5.2) does not rely on the on-shell condition of the right-most scalar k n , hence (5.2) is also valid for amplitudes with off-shell k n . Assuming the validation of (5.2) for YMs amplitude with n < n gluons, we simply get the sum (5.17) as f a n a 1 e i s 1n A Y M s n−1,2 (2, . . . , n − 1, P e 1,n ; p, q), which is identical to the boundary contribution in the LHS of (5.2).
The case with n = 2: this case is much more subtle. The boundary contributions in the LHS of (5.2) come from the diagrams as shown in Fig.(1.b), Fig.(1.c), while the boundary contributions in the RHS of (5.2) come from the diagrams as shown in Fig.(1.d), Fig.(1.e), Fig.(1.f).
According to the Feynman rules for Yang-Mills-scalar amplitudes, we can compute the three terms for the RHS of (5.2) as On the other hand, we can compute the two terms for LHS of (5.2) as If we re-write the second line in the result of Fig.(1.b) by Jacobi identity then the matching of boundary contribution in both sides of (5.2) can be easily checked.
With above discussions, we have confirmed the non-trivial relations between YMs amplitude and pure scalar amplitudes (hence the EYM amplitude and Yang-Mills amplitudes) by BCFW recursion relations. The proof of relations for YMs amplitude with more than two gluons requires more labors, but the strategy is similar, which includes comparing the contributions from finite poles and boundary contributions. We will not discuss it further.

Inspecting the amplitude relations through KLT relation
In the following discussions we will demonstrate that, at least in the first few simplest scenarios, the newly discovered multi-graviton relations [1,14,15] can be readily understood from the perspective of KLT relations. It was demonstrated in [72] that the KLT relation provides a much more perturbation-friendly construction of the EYM amplitudes, which would be otherwise difficult to calculate in viewing of the infinite vertices that constitute the linearized gravity Feynman rules. In this setting, EYM amplitude factorizes into a copy of pure gluon amplitude and a copy that gluon interacts with scalars, through which the color dependence is introduced. To have simpler expression, we will use the (n − 2)! symmetric KLT relation first introduced in [73][74][75], where the numerator in the expression defined using gluon scalar currents carries both kinematic and color factors. The formula defined in (6.2) has provided a way of evaluating the numerator n(1, α, n). However, it is obvious that, directly calculating all currents and then making the sum is not an efficient method. There are two alternative methods to compute the coefficients n(1, α, n).
The first is to carry out the summation step by step as was done in [61,71]. The idea is to divide the full S n−2 permutation sums appearing in (6.2) into (n − 2) blocks of S n−3 permutation sums, such that in each block we can pull out a format of BCJ sums. Then one can use the Fundamental BCJ relation for currents to simplify the expression and arrive at a similar sum as the one given in (6.2) but with only S n−3 permutation sums. Iterating this procedure several times, we can finally compute the coefficients. Establishing the Fundamental BCJ relation for currents is a crucial point for this method, and we will show how to do this in the later sections. The second method is, however, less straightforward. When expanding the amplitude into KK basis with the formulation given in the second line of (6.1), it is shown in [23,54,57,66] that, the coefficients n(1, α, n) are nothing but the numerators of Del Duca-Dixon-Maltoni (DDM) basis provided we write the whole A YMs (1, β, n) amplitude into BCJ form (i.e., numerators satisfying the Jacobi relations). Using this aspect, the problem is translated to computing the BCJ numerators of DDM basis by any conventional methods.
The purpose of this section is to show that, the newly discovered EYM amplitude relations can also be fitted in the framework of KLT relations. The methods that developed in the computation of BCJ numerators in various theories [61,71] are also well-suited in the analysis of EYM amplitude expansion, with only a few modification. This connects the problem of EYM amplitude expansion with many other theories. In the following discussions, we will use both methods developed years ago for computing the BCJ numerators to address the problem of constructing the expansion coefficients n(1, α, n).

The case with single gluon
For the purpose of being self-contained we list the color-ordered Feynman rules for gluon-scalar interaction presented in [72] in the Appendix A. Consider first the scalar Yang-Mills amplitudes when there is only one gluon. Note that a (color-stripped) gluon propagator does not transmit the color/flavor of scalars attached to its two ends, so that for single trace part of the partial amplitude, gluon lines cannot be internal or the color factors carried by the scalars at its two ends factorize. A consequence is that all single gluon amplitudes are consisting of cubic graphs. For example at four points when, say leg 3, is the gluon line there are only three cubic graphs in the KK sector, up to anti-symmetry of the three-vertices,

.
It is very important to notice that the color-kinematics duality is ensured by the vanishing of the sum of their numerators (even at the off-shell level) This observation (i.e., only cubic vertex is allowed and the gauge invariance), when generalized to higher points, indicates that the Feynman diagrams provide the desired BCJ form. In particular, an n-point DDM half-ladder numerator n(1, 2, 3, · · · i, p g , i+1, · · · , n) is therefore given by the corresponding Feynman diagrams as k 2 12 k 2 123 · · · = f 1,2, * f * ,3, * . . . f * ,i,a × δ ab × f b,i+1, * · · · f * ,n−1,n So what are the allowed DDM numerators? The Yang-Mills-scalar theory has a gauge group and a flavor group. The flavor ordering fixes the ordering of scalars, thus the only allowed freedom is the location of gluon leg along the DDM-chain. In other words, the desired expansion coefficients in (6.1) are nothing but the one given in (6.4) with all possible insertions of gluon legs, and we find agreement with the new single graviton relation (up to an overall factor).

The four-point gluon-scalar amplitude involving two gluons
Next we consider Yang-Mills-scalar amplitudes involving two gluons. For this case, since Feynman diagrams will involve the four-point vertex, the BCJ form will not be manifest for Feynman diagrams and the computation will be more complicated. Thus in this subsection, we will follow the method of summing over the color-ordered KK basis.
At four points we have the following two KK basis amplitudes, where n s , n t , n u and n 4 denote the factors and Now that with quartic graph present, the original Jacobi identity inevitably needs to be modified if colorkinematics duality is to remain holding. To better keep track of how this is done we write the BCJ sum of the two KK basis amplitudes in terms of the factors just introduced so that every term appears in the sum has a clear graphical interpretation. Also for future reference we analytically continue one of the scalar legs, say leg 4, and write s 21 A YMs 2,2 (1, 2 g , 3 g , 4) + (s 21 + s 23 )A YMs 2,2 (1, 3 g , 2 g , 4) (6.12) = s 12 n s s 12 − n t s 23 + n 4 + (s 21 + s 23 ) − n u s 13 + n t s 23 + n 4 = (n s + n t + n u + (s 12 − s 13 )n 4 ) + k 2 4 − n u s 13 + n 4 .
The fact that BCJ sum vanishes in the on-shell limit suggests that the Jacobi identity is modified as n s + n t + n u + (s 12 − s 13 )n 4 = 0 up to terms proportional to k 2 4 . A careful inspection shows that these terms actually cancel completely. Plugging equations (6.8) to (6.11) into the left hand side of this modified Jacobi sum, we see that (neglecting an overall factor (−i)δ a 1 a 4 /2), = 0 . (6.14) We obtain the numerator by feeding the off-shell continued BCJ sum just computed into the KLT inspired prescription (6.2), taking the modified Jacobi identity into account, yielding n(13 g 2 g 4) = 1 k 2 4 s 31 s 21 A YMs 2,2 (1, 2 g , 3 g , 4) + (s 21 + s 23 )A YMs 2,2 (1, 3 g , 2 g , 4) = − n u + s 13 n 4 = + s 13 , where by an abuse of notation we neglected factors of inverse propagators so that the graphs appear in the equation above should be understood as representing the corresponding numerators rather than the original Feynman graphs. In the following discussions we shall not distinguish numerators from Feynman graphs unless it is not apparent from the context. The other two gluon numerator at four points can be readily obtained by swapping labels (2 ↔ 3). Inserting the half-ladder numerators back into KLT relation and we find agreement with the two graviton relation (equation (4) in [14]).
Note that the above relation is not exactly the same as (5.1), but equivalent to it after using certain BCJ relations, and note particularly that the new 2 · 3 term came from the quartic graph contribution.

Off-shell continued Jacobi identity
The key point of the above calculation is the modified Jacobi identity when some of the legs becoming off-shell, e.g., n s + n t + n u + (s 12 − s 13 )n 4 = 0. This modification will lead to modified fundamental BCJ relations, to be discussed later. When considering situations for higher points, one note that the color dependency will factorize when the scalars are connected by an internal gluon line, thus the single trace part of a two-gluon partial amplitude can only contain graphs derivable from those appearing at four-point case by welding pure scalar currents to their two scalar legs. Therefore we only need to consider analytically continuing the two scalar lines of the modified Jacobi identity when two gluons are present. Careful inspection of (6.14) shows that the ( · k)( · k) part of the Jacobi sum is a pair-wise cancelation, up to terms proportional to ( 2 · k 2 ) or ( 3 · k 3 ), and therefore remains valid even when scalars become off-shell. The only modification comes from the ( 2 · 3 ) part. To completely cancel the (k 2 − k 3 ) · (k 1 − k 4 ) factor produced by n t , we see that the quartic graph needs to be multiplied by the same factor. The off-shell continued identity we need for all two gluon amplitudes is then and we will be using this identity in the following discussions.

The five-point YMs amplitudes with two gluons
Having presented the example of four points with two gluons, we further show an example of five-point amplitude with two gluons. Again, we will use the method of summing over color-ordered KK basis.
At five points the number of graphs increases considerably. Recall from [8] that there are 15 different graphs in total in the KK sector at five points if the amplitudes are to be described by cubic graphs only, 6 of them are independent when Jacobi identities are taken into account. Similarly we label the cubic graphs as n 1 , n 2 , . . ., n 15 , and we regard the quartic graphs as additional corrections n 16 , n 17 , n 18  and therefore contains quartic graphs, Plugging the above results into DDM expression yields the two graviton EYM amplitude at five points, (1, 2, 4, 3, 5) + · · · (6.38)

The five and higher point amplitude involving two gravitons
Having witness that KLT relation successfully explains the new EYM amplitude expression for two graviton scattering at four and five points, perhaps it is not much of a surprise that the explanation generalizes to higher points. Indeed, one can actually read off the n-point two gluon numerator, and the two graviton EYM amplitude is determined by the corresponding DDM expression. We shall use the algorithm introduced originally for the pure scalar scenario in [71] to systematically calculate the numerator (i.e., to systematically sum over KK basis). As we shall see, in the case when only two gluons (p, q) are involved, the numerators remain fairly simple, n(12 · · · i p g · · · j q g · · · n) = + 2(p · Y p ) (6.39) A brief review of the algorithm for numerators in the scalar scenario: For the purpose of being self-contained, we briefly review the algorithm used by the authors in [71] and [61] to calculate numerators. The idea is to divide the full S n−2 permutation sum appears in the numerator-current relation n(1 α n) = β∈S n−2 S[α T |β] J YMs (1, β, n) into BCJ sums, and proceed repeatedly if the Fundamental BCJ relation between currents admits further simplifications. For example, it was shown in [71] that the Fundamental BCJ relation between φ 3 currents yields another current, with the leg running through all insertions in the BCJ sum fixed at the off-shell continued line, so that if we divide the full permutation sum S 3 in the five-point numerator calculation into BCJ sums, after substituting these summations using Fundamental BCJ relation (6.40), the collected result is yet another BCJ sum, but only performed over permutations of the legs of fewer-point sub-currents, where we used (6.40) to replace the first and the second line of the equation above with the two graphs in (6.41). The result is another BCJ sum over currents. Repeat the substitution using Fundamental BCJ relation, and we obtain the numerator The five-point scalar Yang-Mills numerators involving two gluons The calculation explained above only complicates slightly when few gluons are present. As far as single trace contributions are concerned, all amplitudes are consisted of Jacobi satisfying cubic graphs when only one gluon participates the scattering, and the same algorithm applies. It is straightforward to see that the numerator is given by n(12 · · · p g · · · n) = f 1,2, * f * ,3, * · · · f * ,n−1,n p · x p , which when plugged into the summing expression readily reproduces the new EYM formula. In other words, (6.4) can also be understood from this point of view.
Things will become a little bit more complicated when two and more gluons are involved, since quartic vertices start to come into play, although they still remain quite manageable, in the sense that the modified Fundamental BCJ relations brought by the quartic term also permit repeated use of the relation when we carry out the summation. Explicitly, at five points the two-gluon Fundamental BCJ relations are modified as . . .
The rules to modification is as follows. Generically one only needs to replace the appropriate scalar by gluon lines in the original Fundamental BCJ relation between currents (6.40), and the right hand side of the equation is a current with the running leg fixed at the off-shell line. The only exception is when the running leg is gluonic, also that either leg 1 or leg n − 1 (legs adjacent to the off-shell line) is a gluon line. In these cases an additional current needs to be added, where a quartic vertex resides on the off-shell line connects the two gluons.
We leave the details of a proof to these relations at five-point to Appendix B because of its complicated nature. The principles are however not much different from the pure scalar scenario and is conceptually straightforward. Basically we cancel graphs related by Jacobi identities among Berends-Giele decomposed five-point current in the BCJ sum. The result after cancelation is then collected and identified to be the Berends-Giele decomposition of the right hand side of the equation. The proof for generic n points follows rather trivially from the structure of the proof, since adding more scalar lines into sub-currents at peripherals does not change Jacobi identities.
Assuming the Fundamental BCJ relations above, it is not difficult to see that the numerator is genuinely given by the formula (6.39) we claimed earlier. Consider for example the derivation that leads to numerator As was explained earlier we obtain the numerator by first dividing the full S n−2 permutation sum appears in the KLT inspired prescription (6.2) into BCJ sums, and then use the Fundamental BCJ relation between currents to fix the n − 2 legs one by one in descending order. For the most part, this procedure is not different from the derivation of a pure scalar numerator, and the result does contain a cubic half ladder graph. The only modification occurs whenever the leg we attempt to fix is gluonic, in which case an additional graph is included, where a quartic vertex connecting both gluon lines emerges. The derivation afterwards again follows that of a pure scalar numerator. In the n(123 g 4 g 5) example this leads to n(123 g 4 g 5) = k 2 1234 k 2 123 k 2 12 + k 2 1234 k 2 12 (s 31 + s 32 ) . (6.46) Note that the Mandelstam variables associated with the quartic graph was furnished by momentum kernel. Careful inspection of the derivation that leads to (6.45) shows that they should contain the inner products between gluon line carrying the smaller label and all the scalar lines which precede it. As another illustration we consider n(12 g 34 g 5), As a verification, note that applying the same rules to derive numerators with all possible combinations of gluon positions yields the same results as those listed from equation (6.32) to (6.37) previously obtained exclusively for five points.

Conclusion
In this paper, we studied the newly discovered EYM amplitude relation by gauge invariance principle, the BCFW recursion relation as well as the KLT relation respectively. It turns out that the problem of EYM amplitude expansion is also closely related to the problem of computing BCJ numerators and the boundary contribution of BCFW terms.
The major context of this paper is devoted to the principle of gauge invariance applied to the determination of EYM amplitude relations. We propose a constructive algorithm by expanding any EYM amplitude A EYM n,m as a linear sum of A EYM n+i,m−i for i = 1, 2, . . . , m with given expansion coefficients, and the contributing terms of A EYM n+i,m−i are determined by A EYM n+i−1,m−i+1 . This means that any contributing terms can be recursively determined by the very first one A EYM n+1,m−1 , while keeping the gauge invariance in each step. This leads to a compact formula (3.39) for general EYM amplitude relations with arbitrary number of gravitons. Realizing that the expansion of Einstein-Yang-Mills amplitude into Yang-Mills amplitudes shares the same kinematic coefficient as the expansion of Yang-Mills-scalar amplitude into cubic-scalar amplitudes, we copy the EYM amplitude relation to YMs amplitude relation, and generalize the later one to the expansion of pure Yang-Mills amplitude into cubic-scalar amplitudes by the help of Pfaffian expansion. With the Yang-Mills amplitude expanded recursively into the cubic graphs, we further outline the strategy of rewriting the scalar amplitudes into KK basis, manifesting the color-kinematics duality and computing the BCJ numerators of Yang-Mills amplitude.
We also study the EYM amplitude relations in the S-matrix framework, and present the proof of EYM amplitude relations with two gravitons by BCFW recursion relations. In this case, any choice of deformed momenta is not possible to avoid the boundary contributions, so we need to compare the contributions of both sides in the relations from finite poles and also the boundary. The matching of both contributions also constraints the possible form of the non-trivial relations. Besides, we examine the problem again from the perspective of KLT relations. The expansion coefficients of EYM amplitude relations are identical to the BCJ numerators of DDM basis, and by computing the BCJ relations for currents we confirm the validation of EYM amplitude relations.
Following our results, there are many interesting directions to explore further. In our paper, one of the most important results is the recursive construction (3.39) of EYM amplitude relation. We have claimed this expression by a few explicit examples plus the guidance of gauge invariance principle. For the confirmation of the claim, a rigorous derivation by other methods is favorable. In an upcoming paper, we would explore the recursive construction directly from operations on the CHY-integrand level. Furthermore, we believe that, such recursive pattern can also find its hints in the BCFW recursion relation or KLT relation investigation of EYM amplitude expansion, which worth to work on with.
Another possible work would be that, in our recursive construction, the gauge invariance is manifest for all gravitons at each step except the first one that started the recursive algorithm. As shown in [45,[47][48][49], for Yang-Mills theory, the requirement of gange invariance for (n − 1) points is sufficient to guarantee the correctness of the full amplitude. This observation seems to be also true in the EYM theory, thus finding an explicit proof along the same line as in [47][48][49] would be a thing worth to do.
A most interesting and important future direction would be the systematic study of the CHY-integrand expansion. In §2, we have laid down the general framework for the expansion, while in the whole paper we are focus only on the expansion of (reduced) Pfaffian. However, many CHY-integrands, such as (Pf (A)) 2 can be obtained from Pfaffian with proper reduction. Thus our results could be easily generalized to many other theories. Especially by similar calculations, we can check if the soft theorem can be used to uniquely determine the amplitude for some theories, such as NLSM as advertised in [47][48][49].
Finally, as a byproduct of the EYM expansion, we have outlined the strategy of computing BCJ numerators 4 from the expansion relation for general EYM amplitudes. The four-point example shows the procedure of computing the BCJ numerators as polynomial of ( · ), ( · k) and (k · k), constructed neatly from the expansion coefficients of EYM amplitudes into Yang-Mills amplitudes. This construction, when generalized to loop-level, would fascinate many important calculations involving gravitons. B Graphical proof of the two-gluon Fundamental BCJ relation between currents As a demonstration of the general idea, in this appendix we prove two of the Fundamental BCJ relations at five points involving two gluons, equations (6.43) and (6.44), following the method used in [71] (which was also briefly outlined earlier in §6.4). We shall neglect repeating a similar proof for generic n points, as it can be readily derived by induction and by attaching more external legs on the sub-currents represented by blank circles in the graphs below.

Relations with no gluon adjacent to the off-shell leg
Consider first the configuration where the leg running over all possible insertions in the BCJ sum is a gluon, and the other gluon is non-adjacent to the off-shell leg. We would like to prove that s 21 J YMs (12 g 3 g 45) + (s 2,13 )J YMs (13 g 2 g 45) + (s 2,134 )J YMs (13 g 42 g 5) = . (B.1) For this purpose we Berends-Giele decompose all three currents appear in the BCJ sum, yielding altogether nine graphs, and notice that, aside from (c1), rest of the graphs can be regrouped as BCJ sums of sub-currents. Indeed, graphs (a1)and (b1) together make up a BCJ sum of the sub-currents involving legs 1, 2 and 3, and graph (a2) is by itself a (trivial) BCJ sum of three point current. The combination of (b2) and (c2) is also a BCJ sum of the three point current, after eliminating part of the sum that carries an s 21 using U (1) decoupling identity, In the equations above we have assumed the BCJ relations between currents at four points. As for the remaining graph (c1) that does not regrouped with the others into a BCJ sum of sub-currents, we rewrite the coefficient it carries using the kinematic identity s 2,134 = k 2 5 − s 134 and then further Berends-Giele decompose the part that carries a factor s 134 , giving Relations with one gluon adjacent to the off-shell leg The proof when one gluon is adjacent to the off-shell leg follows exactly the same derivation, except that now we have a few additional quartic graphs. The relation we are aiming to prove is s 21 J YMs (12 g 34 g 5)+s 2,13 J YMs (132 g 4 g 5)+s 2,134 J YMs (134 g 2 g 5) = k 2 .

(B.15)
Note the presence of two new quartic graphs (b4) and (c4). As in the previous example we regroup graphs into BCJ sums of sub-currents. Graphs (a1) and (b1) make up a BCJ sum of the sub-currents involving legs 1, 2 and 3,  , which completes our proof.