Color-Kinematic Numerators for Fermion Compton Amplitudes
N. Emil J. Bjerrum-Bohr Gang Chen Yuchan Miao Marcos Skowronek



























hobby,calc,shapes, positioning, backgrounds,trees, decorations, decorations.markings, decorations.pathmorphing, arrows, shapes.geometric, snakes

compat=1.1.0

graviton/.style=circle, draw=green!60, fill=green!5, very thick, minimum size=7mm


codot/.style=/tikz/shape=circle,

/tikz/fill=red,/tikz/minimum size=0.1cm,/tikz/inner sep=1.8pt

sb/.style=
/tikz/shape=circle,
/tikz/fill=black,
 /tikz/minimum size=0.3cm,/tikz/inner sep=1.pt
  
 rsb/.style=
/tikz/shape=circle,
/tikz/fill=red,
 /tikz/minimum size=0.3cm,/tikz/inner sep=1.pt
  
myblob/.style=
/tikz/shape=ellipse,
/tikz/fill=red,
 /tikz/minimum width=0.5cm,
  
 myblobFF/.style=
/tikz/shape=circle,
/tikz/fill=black,
 /tikz/minimum width=0.25cm,
  
 myblob2/.style=
/tikz/shape=rectangle,
/tikz/fill=black,
 /tikz/minimum width=0.1cm,/tikz/inner sep=1.8pt
  
  fp/.style=
/tikz/shape=rectangle,
/tikz/fill=red,
 /tikz/minimum width=0.1cm,/tikz/inner sep=1.8pt
  
  myvertex/.style=
/tikz/shape=circle,
/tikz/fill=black,
 /tikz/minimum width=0.05cm,
  
box/.pic=[fill=black]  (0,0) circle (2.5pt);
				   [fill=black] (0.5,0) circle (2.5pt);
			       [line width=5pt] (0,0) – (0.5,0);


wiggle/.style=decorate, decoration=snake

 HV/.style=
/tikz/shape=circle,
/tikz/fill=rgb:black,1;white,2,
 /tikz/minimum size=0.01cm,/tikz/inner sep=1.8pt
  




1.5pt
 
theoremTheorem
definitionDefinition[section]
lemma[theorem]Lemma
conjecture[theorem]Observation
prop[theorem]Proposition
.theorem











#1#1


#1
 #1
  to1ex[#1]
  
  ·
  

η̃
η̅
m
λ̃







#1eq. (<ref>)
#1Equation (<ref>)
#1#2eqs. (<ref>) and (<ref>)
#1#2Eqs. (<ref>) and (<ref>)

#1section <ref>

⟨#|1⟨⟨#1|
|#⟩1|#1 ⟩⟩
#1[#1|
#1|#1]
⟨#|1⟩⟨⟨#1 ⟩⟩

#1.#2⟨⟨#1 #2⟩⟩
#1.#2[#1 #2]



    #1 #1#1 #1#1 #1
        
        *
          @
        #1𝖪_#1
         ^aNiels Bohr International Academy,
        Niels Bohr Institute, University of Copenhagen,
        
        Blegdamsvej 17, DK-2100 Copenhagen Ø, Denmark^bDepartment of Physics, Brown University, Providence, RI 02912, USAbjbohr@nbi.dkgang.chen@nbi.ku.dksxh962@alumni.ku.dkmarcosskowroneksantos@brown.edu
        We introduce a novel approach to compute Compton amplitudes involving a fermion pair inspired by Hopf algebra amplitude constructions. This approach features a recursive relation employing quasi-shuffle sets, directly verifiable by massive factorization properties. We derive results for minimal gauge invariant color-kinematic numerators with physical massive poles using this method. We have also deduced a graphical method for deriving numerators that simplifies the numerator generation and eliminates redundancies, thus providing several computational advantages.
        Color-Kinematic Numerators for Fermion Compton Amplitudes
            Marcos Skowronek^b
            Received xx; accepted xx
        =========================================================
        
        
        
        § INTRODUCTION
        
        The Kawai-Lewellen-Tye string theory relations <cit.> combine gravity and Yang-Mills theories <cit.> and have led to the important discovery of a color and kinematics duality <cit.>; a cornerstone for efficient computation of amplitudes and motivation for examining core principles of quantum particle scattering <cit.>.  
        
        A considerable attraction over a decade has been finding a kinematic algebra that mimics that of the color group <cit.>, and recently, there has been a breakthrough with the systematic framework developed by refs. <cit.>. This framework uses infinite-dimensional combinatorics <cit.> and a direct numerator construction approach connected to quasi-shuffle Hopf algebras. The numerator technology developed this way is universal across many theories, including those with a spinless massive field, such as the heavy-mass effective field theory refs. <cit.>, scalar Yang-Mills theory with finite mass <cit.>, the α'F^3+α'^2F^4 higher-derivative interaction theory <cit.> (see also refs. <cit.>), and the DF^2+YM theory <cit.> (see as well refs. <cit.>).
        
        This paper studies the kinematic algebra for the gauge invariant color-kinematic numerator associated with the generalized massive Compton amplitude with a fermion pair. Although this amplitude is considerably more complex than the spinless massive field case (see also refs. <cit.>), we show the realization of the algebra as a closed-form quasi-shuffle product. The amplitudes discussed are of theoretical interest as they allow for structural insight into fermion factorization limits and are useful when generating inspiration for bootstrap methods for general spin amplitudes <cit.>. We thus hope that progress in this area could potentially help develop computational methods for spinning black hole merger observables, for instance, using the modern amplitude techniques pioneered in refs. such as <cit.>.
        
        We present our work in this paper as follows. First, we introduce the kinematic Hopf algebra in section <ref> and define a new evaluation map for massive Compton amplitudes with fermion pairs. Next, we demonstrate the validity of this formalism through computation and employ factorization properties to outline a proof of the construction.
        In section <ref>, we provide an evaluation map that simplifies numerators using the on-shell condition of the external massive fermion states. This graphical method streamlines numerator generation and thereby delivers computational benefits.
        Finally, in the final section, we provide an outlook for future research.
        
        
        
        § COLOR-KINEMATIC FERMION PAIR NUMERATORS
        
        This section will deduce and demonstrate our algebraic construction for numerators. We will start by setting up a convenient and compact notation. A key element is an n-point massive fermion amplitude representation in terms of minimal color-kinematic numerators with (n-3)! elements involving nested commutators. As an illustration of this, we consider ordered cubic vertex diagrams associated with gluon lines, 1, 2, and 3 connected through a propagator to gluon legs 4,5,…,n-2 and two massive scalars n-1, n.
        
        
            [baseline=([yshift=-0.8ex]current bounding box.center)]every node=[font=]    
                (a)[];
                   [above=0.9cm of a](b)[dot];
                  [left=0.7cm of b](c);
                  [left=0.28cm of b](c23);
                  [above=0.23cm of c23](v23)[dot];
                 [above=.5cm of c](j1)p_1;
                 [right=.9cm of j1](j2)p_2;
                 [right=0.6cm of j2](j3)p_3;
               	 *(a) – [thick] (b),(b) – [thick] (j1),(v23) – [thick] (j2),(b)–[thick](j3);
                      [baseline=([yshift=-0.8ex]current bounding box.center)]every node=[font=]    
                (a)[];
                   [above=0.9cm of a](b)[dot];
                  [left=0.7cm of b](c);
                  [right=0.29cm of b](c23);
                  [above=0.23cm of c23](v23)[dot];
                 [above=.5cm of c](j1)p_1;
                 [right=.6cm of j1](j2)p_2;
                 [right=0.9cm of j2](j3)p_3;
               	 *(a) – [thick] (b),(b) – [thick] (j1),(v23) – [thick] (j2),(b)–[thick](j3);
        
        The two possible propagators respectively are p_12^2p_123^2 and p_23^2p_123^2, (defining p_ij≡ p_i+p_j and p_ij… k≡ p_i+p_j+...+p_k) and it is easy to see that the color-kinematic duality holds writing out the Jacobi identity for numerators. 
        
            A_f(1,2,3,n-1, n)    =_f([[1,2],3],n-1, n) p_12^2p_123^2+_f([1,[2,3]],n-1, n) p_23^2p_123^2
            =∑_∈ρ(3)_f(,n-1, n) d_ ,
        
        where Γ∈ρ(3)≡{[[1,2],3],[1,[2,3]]} and d_ are massless propagators associated with numerators. This notation readily extends to any number of legs as follows
        
        
            A_f(1… n-2,n-1, n)     = ∑_∈ρ(n-2)_f(,n-1, n) d_ ,
        
        
        where now Γ sums over the possible (n-3)! configurations of ordered nested numerators associated with the n-2 gluon legs.
        We can also compute graviton amplitudes using the double-copy. In this case, we write the n point amplitude as
        
        
            M_f(1… n-2,n-1, n)     = ∑_∈ρ̃(n-2)_f(,n-1, n)(, n,n-1)  d_  .
        
        We denote sums over configurations by ρ̃(n-2), where the `tilde' indicates that the sum includes all possible nested unordered commutator configurations associated with gravitons labels { 1, …  n-2 }.
        
        Using the gravity momentum S-kernel <cit.>, we know that we can always generate generic color-kinematic numerators using the following prescription  <cit.>
        
            _f([[[[⋯[1,2],3],…], n-2],n-1,n)=∑_β∈ S_n-3𝒮[1… n-2 |1β ]A(1β,n-1,n)  .
        
        In the sum, S_n-3 denotes permutations of gluon indices {2,3,…, n-2} and we work with the momentum S-kernel definition employed in <cit.>. We will, for compactness in the following, denote the left-nested commutator of gluons as [12… n-2],  i.e., [1234]≡[[[1,2],3],4]. (For more details, see also refs. <cit.>). When considering factorization, the following recursive property of the gravity momentum S-kernel is ideal 
        
        
            𝒮[1… j|1β_L jβ_R]=2(p_Θ_L(j) p_j) 𝒮[1… j-1|1β_L β_R]   .
        
        In this context, the symbol Θ_L(j) denotes the indices that come before the index j in the set {1, β_L}, with an index order that is lower than that of j. With these conventions, the expression 𝒮[1α |1β] only depends on the massless momentum and cancels out all the massless propagators in the amplitude, although not manifestly. 
        
        
        
         §.§ Numerator examples at four and five points
        
        Before diving into the kinematic algebra, we will explore a manifestly gauge-invariant double-copy numerator form of the quark-gluon Compton amplitude derived from Feynman rules. We have adopted the following convention for the three-point amplitudes
        
        
            [baseline=(current bounding box.center)]
               (a)3;
                    [right=2.cm of a](b)2;
                   [right=1.0cm of a](c);
                   [above=1.0cm of c](d)1;
                   *  
                             (b) – [fermion] (c) – [fermion] (a),
                             (c) – [gluon] (d)
                            ;
               
            = 1/2v̅ε_1 u  ,       [scale=0.8, transform shape, baseline=(current bounding box.center)]
             (c0);
               [above=1.2cm of c0] (a) 1; 
               [below left=0.55cm and 1.0cm of c0] (b) 3;
               [below right=0.55cm and 1.0cm of c0] (c) 2; 
              * 
                (c0) – [gluon] (a),
                (c0) – [gluon] (b),
                (c0) – [gluon] (c)
              ;
            
            0.3
                =(ε_1ε_3)(p_1 ε_2)+(ε_2ε_3)(p_3 ε_1)
                        +(ε_1ε_2)(p_2 ε_3) ,
        
        where ≡γ^μ_μ, for a given gluon polarization vector, v̅, and u denote  the incoming and outgoing  massive fermions states with on-shell conditions v̅(p_3-m)=0   , 
            (p_2+m) u=0,   and mass m.
        The external momenta are considered outgoing. 
        Our focus is on the first non-trivial case, the four-point quark-gluon amplitude, which is featured by two contributing color-ordered Feynman diagrams:
        
        
            [baseline=([yshift=-0.4ex]current bounding box.center)]every node=[font=] 
                 (4) 4;
                    [right=1cm of 4] (14);
                    [right=1cm of 14] (23);
                    [right=1cm of 23] (3)3;
                    [above=1cm of 14] (1) 1;
                    [above=1cm of 23] (2) 2;
                        * [edges=fermion]
                             (3) – (23) – (14) – (4)
                            ,
                                (23) – [gluon] (2),
                                (14) – [gluon] (1),
                                    ;
                
                =1 4 (p_14^2-m^2)
                (v̅ε_1 (p_14+m) ε_2u)   ,
                
               
                [baseline=([yshift=-2.1ex]current bounding box.center)]every node=[font=] 
                 (4) 4;
                    [right=1.5cm of 4] (1234);
                    [right=1.5cm of 1234] (3)3;
                    [above=1cm of 1234] (12);
                    [above right=1cm of 12] (2)2;
                    [above left=1cm of 12] (1)1;
             * 
             [edges=fermion]
             (3) – (1234) – (4)
             ,
            (1234) – [gluon] (12),
             (12) – [gluon] (1),
             (12) – [gluon] (2),
                    ;
                   =
            -1/2p_12^2[(p_2ε_1)(v̅ε_2 u)-(p_1ε_2)(v̅ε_1 u)
                          +1/2(ε_1ε_2)(v̅(p_1-p_2) u)]   .
        
        The four-point KLT momentum kernel matrix as indicated by  eq. (<ref>) provides a double-copy numerator of the form (we use 𝒮[12|12]= 
        2(p_1 p_2))
        
        
            𝒩_f([1 2], 3, 4)     = 
            2(p_1 p_2) A_f(1 2, 3, 4)=
                1/2(p_14^2-m^2)v̅(
               p_1  p_2ε_1(p_14+m)ε_2 
               -(p_1 p_4)
                [
                2(p_2ε_1)ε_2-2(p_1ε_2)ε_1+(ε_1ε_2)(p_1-p_2)
                ]
                ) u  .
        
        After simplification, this leads to a compact formulation of the double-copy numerator, expressed in terms of the field strength tensor:
        
        
            𝒩_f([12], 3, 4)=-1 p_41^2-m^2(
               v̅(p_4 F_12)  v+1 4u̅ F_1  (p_1  F_2)  u)   .
        
        In this context, we use F_i^μν to denote the Abelian field strength tensor, so that F_i^μν≡ p_i^με_i^ν-ε_i^μ p_i^ν. We represent the Dirac matrix product with a bullet symbol , and we employ the simplified notation as follows:
        
        
            v̅ X (p_1· F_2) Y u    ≡ (v̅ Xγ_μ Y  u)  (p_1· F_2)^μ,
            v̅ X F_i_1… i_r Y u    ≡  (v̅ Xγ_μγ_ν Y u)  F_i_1… i_r^μν  .
        
        Here F_i_1i_2… i_r^μν≡ (F_i_1· F_i_2⋯ F_i_r)^μν denotes the contraction of a string of multiple field strength tensors. 
        
        Next, we turn our attention to the five-point quark-gluon amplitude, which we can compute from the four diagram topologies:
        
        
            [scale=0.75, transform shape, baseline=(1.base)]
             (1) 5;
                [right=1cm of 1] (12);
                [right=1cm of 12] (123);
                [right=1cm of 123] (1234);
                [right=1cm of 1234] (5)4;
                [above=1.7cm of 12] (2) 1;
                [above=1.7cm of 123] (3) 2;
                [above=1.7cm of 1234] (4) 3;
                * [edges=fermion]
                        (5)–(1234)–(123)–(12)–(1),
                    ,
                    (12) – [gluon] (2),
                    (123) – [gluon] (3),
                    (1234) – [gluon] (4),
                ;
            
            [scale=0.75, transform shape, baseline=(1.base)]
             (1) 5;
                [right=1.33cm of 1] (12);
                [right=1.48cm of 12] (123);
                [right=1.33cm of 123] (5)4;
                [above=1.7cm of 12] (2) 1;
                [above=1cm of 123] (34);
                [above right=1cm of 34] (4) 3;
                [above left=1cm of 34] (3) 2;
                * [edges=fermion]
                        (5)–(123)–(12)–(1),
                    ,
                    (12) – [gluon] (2),
                    (123) – [gluon] (34),
                    (34) – [gluon] (3),
                    (34) – [gluon] (4),
                ;
            
            [scale=0.75, transform shape, baseline=(a.base)]
             (a) 5;
                [right=1.33cm of a] (ab);
                [right=1.48cm of ab] (c);
                [right=1.33cm of c] (d) 4;
                [above=1cm of ab] (f12);
                [above right=1cm of f12] (e) 2;
                [above left=1cm of f12] (f) 1;
                [above=1.7cm of c] (g) 3;
                * [edges=fermion]
                        (d)–(c)–(ab)–(a),
                    ,
                    (ab) – [gluon] (f12),
                    (f12) – [gluon] (f),
                    (f12) – [gluon] (e),
                    (c) – [gluon] (g),
                ;
            
            [scale=0.75, transform shape, baseline=(a.base)]
             (a) 5;
                [right=2cm of a] (b);
                [right=2cm of b] (c) 4;
                [above=0.77cm of b] (d) [myblobFF] A;
                [above left=1.5cm and 0.75cm of d] (g1) 1;
                [above=1.5cm of d] (g2) 2;
                [above right=1.5cm and 0.75cm of d] (g3) 3;
                * 
                    (c) – [fermion] (b) – [fermion] (a),
                    (b) – [gluon] (d),
                    (d) – [gluon] (g1),
                    (d) – [gluon] (g2),
                    (d) – [gluon] (g3),
                ;
        
        Building upon eq. (<ref>), we derive the five-point double-copy numerator as follows:
        
        
            𝒩_f([123], 4,5)   =𝒮[123|123]  A_f(123, 4,5)+𝒮[123|132]  A_f(132, 4,5)   ,
        
        where the elements of the KLT momentum kernel, as outlined in eq. (<ref>), are 𝒮[123|123]=4(p_1 p_2)  (p_12 p_3) and  𝒮[123|132]=4(p_1 p_2) ( p_1 p_3).
        This leads to manifestly gauge invariant and local double-copy five-point numerator
        
         
            𝒩_f([123],4,5)
               =v̅[ (p_5 F_123) + 1 4F_1(p_1  F_23) p^2_51-m^2 
               -(p_5 F_12) + 1 4F_1(p_1  F_2) p^2_51-m^22(p_12 F_3 p_5)+(p_12 F_3) p_3 p^2_512-m^2
               - (p_5 F_13) + 1 4F_1(p_1  F_3) p^2_51-m^22(p_13 F_2 p_5)+(p_13 F_2) p_2 p^2_513-m^2 ]
             u   .
        
        
        Systematically obtaining such expressions for higher point amplitudes is challenging due to the intricate process of eliminating massless poles in the double-copy numerator and rewriting in terms of strings of field strength tensors. 
        In the following section, we will work out an algebraic construction of numerators using Hopf algebra, which offers a direct and systematic universal approach.
        
        
        
         §.§ Direct algebraic construction of numerators
        
        We will now use Hopf algebra as an inspiration to simplify the construction of numerators and avoid the complexities of the direct computation approach. We will use pre-numerators, which follow the logic of <cit.>. 
        We note that pre-numerators are not uniquely determined and still have residual freedom. However, we can eliminate it by working in a framework where we equate the left-nested numerator _f([12… n-2], n-1,n) with the pre-numerator _f(12… n-2, n-1,n) and fix all numerators _f(… 1 …, n-1, n) to be zero.[We remark that working directly with crossing symmetric numerators will lead to the same result.]. We next compute the pre-numerator using the fusion product
        
        
            _f(12… n-2,n-1,  n)   =⟨⟨ T_(1)⋆ T_(2)⋆⋯⋆ T_(n-2)⟩⟩   ,
        
        where we label the generators  e.g.T_(i) and T_(τ_1),(τ_2),…, (τ_r) and the subscripts τ_i⊂{1, …, n-2} represent an ordered partition of a subset of the external gluon indices. The evaluation map denoted by ⟨⟨…⟩⟩ above provides the physical gauge-invariant object that appears in the color-kinematic numerator. We can work out the fusion products by introducing quasi-shuffle operator products given by:
        
        
            T_(1)⋆ T_(2)=-T_(12)+T_(1),(2)+T_(2),(1)
               …
               T_(τ_1),…, (τ_r)⋆ T_(j)=-∑_i=1^rT_(τ_1),…,(τ_i j),…, (τ_r)+∑_σ∈{τ_1,…,τ_r} j^r (T_(σ_1),…, (σ_r),(σ_r+1))   .
        
        leading to the following explicit result for the four-point numerator:
        
        
            _f(12,3, 4)=⟨⟨ T_(1)⋆ T_(2)⟩⟩=-⟨⟨ T_(12)⟩⟩ +⟨⟨ T_(1),(2)⟩⟩+⟨⟨ T_(2),(1)⟩⟩  .
        
        We now turn our focus to the evaluation that maps the abstract generators T_(τ_1),(τ_2),…, (τ_r) to gauge invariant functions. For our massive fermion amplitudes, it reads
        
        
            ⟨⟨ T_(1τ_1),(τ_2),…,(τ_r)⟩⟩ = ⟨⟨[baseline=([yshift=-0.8ex]current bounding box.center)]every node=[font=]    
                (a)[myblob];
                 [left=0.8cm of a] (a0)[myblob];
                 [right=0.8cm of a] (a2)[myblob];
                  [right=0.8cm of a2] (a3)[myblob];
                   [right=0.8cm of a3] (a4)[myblob];
                   [above=0.8cm of a] (b1)τ_1;
                    [above=0.8cm of a2] (b2)τ_2;
                    [above=0.8cm of a3] (b3)⋯;
                     [above=0.8cm of a4] (b4)τ_r;
                     [above=0.8cm of a0] (b0)1;
                    [above=0.8cm of a0](j1);
                 [right=0.2cm of j1](j2);
                 [right=0.6cm of j2](j3);
                 [right=0.4cm of j3](j4);
                 [right=0.8cm of j4](j5);
                   [right=0.2cm of j5](j6);
                 [right=0.6cm of j6](j7);
                  [right=0.0cm of j7](j8);
                 [right=0.8cm of j8](j9);
               	 *(a)–[very thick](a0),(a)–[very thick](a2),(a2)–[very thick](a3), (a3)–[very thick](a4),(a0) – [thick] (j1),(a) – [thick] (j2),(a)–[thick](j3),(a2) – [thick] (j4),(a2)–[thick](j5),(a4) – [thick] (j8),(a4)–[thick](j9);
                ⟩⟩
               = v̅H_1τ_1 p^2_n1-m^2p_n1τ_1-m p^2_n1τ_1-m^2  (p_Θ_L(τ_2) F_τ_2)⋯p_n1τ_1…τ_r-1-m  p^2_n1τ_1⋯τ_r-1-m^2  (p_Θ_L(τ_r) F_τ_r) u  ,
        
        
        where we have T_(1),…= T_…,(1τ_i),… =0 and 
        
            H_1τ≡ p_n F_1τ + 1 4F_1(p_1  F_τ).
        
        The symbol Θ_L(τ_r) denotes all indices in (1τ_1)∪ (τ_2)∪⋯∪ (τ_r-1) that are smaller than the first element in τ_r. We can write a closed-form color-kinematic numerator as follows: 
        
        
            _f(12… n-2, n-1, n)=∑_ ordered  partitions(-1)^n+r⟨⟨ T_(1τ_1),(τ_2),…,(τ_r)⟩⟩ ,
        
        
        where sum is over ordered partitions of the set {2,3,...,n-2} into r nonempty subsets[More examples of color-kinematic numerators and a  Mathematica package can be found at  https://github.com/AmplitudeGravity/kinematicHopfAlgebra KiHA5.0 GitHub repository <cit.>.]. We can now expand this sum as follows.
        
        
            _f(1… n-2,n-1,  n)=(-1)^n-3 v̅ H_1… n-2 p^2_n1-m^2 u
               +∑_τ_R⊂{2,…,n-2} (-1)^|τ_R| v̅[  X(1ρ(τ_R)) (p_n1ρ(τ_R)-m) (p_Θ_L(τ_R) F_τ_R)  p^2_n1ρ(τ_R)-m^2] u  ,
        
        where we have defined |τ_R| as the number of elements in the set τ_R, and ρ(τ_R) as the complement legs in {2,…, n-2}. X(1ρ(τ_R)) represents a particular Dirac matrix product of  H functions (defined above) and massive propagators for the set of gluon legs {1ρ(τ_R)} that correspond to a color-kinematic numerator that is stripped of incoming and outgoing external fermion states. From this, we introduce a factorization on an internal massive leg denoted by I and give the expression,
        
        
            _f(1… n-2,   n-1,  n) =(-1)^n-3 v̅ H_1… n-2 p^2_n1-m^2 u  
               +∑_τ_R⊂{2,…,n-2}∑_S_I(-1)^|τ_R| _f(1ρ(τ_R), I, n) v̅_I  (p_Θ_L(τ_R) F_τ_R)  p^2_n1ρ(τ_R)-m^2 u   .
        
        where S_I denotes the possible spin states of I and where we have used in the last step that v̅ X(1ρ(τ_R)) u_I ≡_f(1ρ(τ_R), I, n). 
        
        To prove the recursive form shown in eq. (<ref>), we start by considering the factorization behavior on the massive propagator that involves only one massless particle j on the right-hand side. We then take the collinear limit, in legs j and n-1,
         
        
            _f(1… n-2, n-1, n)
            ⟶[baseline=([yshift=-0.8ex]current bounding box.center)]
            every node=[font=]  
            mybox = [rectangle, draw=black, fill=red, minimum size=0.3cm]
            mydot = [circle, draw=black, fill=black, minimum size=3pt,inner sep=0pt]
             (a)[mybox];
             [left=0.8cm of a] (a0);
             [right=0.8cm of a] (cut);
                [below=0.2cm of cut] (I)I;
             [below=0.2cm of a0] (n)n;
             [right=2.4cm of a] (a1);
             [right=1.6cm of a] (j0)[mydot];
             [red, thick] ((cut) + (-0.1,-0.1)) – ((cut) + (0.1,0.1));
             [red, thick] ((cut) + (-0.1,0.1)) – ((cut) + (0.1,-0.1)colin);
             [above=0.9cm of j0] (j1)j;
             [below=0.2cm of a1] (n1)n-1;
             [above=0.9cm of a] (b)1     ρ(j);
             [above=0.8cm of a0] (b0);
             [right=0.4cm of b0] (b11);
             [right=0.8cm of b11] (b12);
             *(a)–[very thick](a0),(a)–[very thick](a1),(a) – [thick] (b11),(a)–[thick](b12),(j1)–[thick](j0);
            
               = ∑_S_I∑_β_L∈ S_n-4 𝒮[1ρ(j)|1β_L] A(1β_L, I, n) (p_Θ_L(j) p_j)   A(j,n-1, I)
               =∑_S_I _f(1ρ(j), I, n)  (p_Θ_L(j) p_j)  _f(j,n-1, I)   .
        
        We have in the above equation utilized eq. (<ref>) and  eq. (<ref>).
        Using the completeness relation for Dirac spinors and explicitly writing the expression for the fermionic propagator, the sum over intermediate states yields the following after evaluating the three-point numerator on the right side of the cut:
        
        
            _f(1… n-2, n-1, n)→ (p_Θ_L(j) p_j) v̅ X(1ρ(j)) (p_n1ρ(j)-m) _j  u   .
        
        We can now deduce a manifest gauge invariant part of the numerator that replicates this factorization behavior. 
        
        
            ^(1)(1ρ(j)|j, n-1,n)
                   = -v̅ X(1ρ(j)) (p_n1ρ(j)-m)  (p_Θ_L(j) F_j) u  p_n1ρ(j)^2-m^2  ,
        
        where ρ(j)|j indicates that the particle to the right of ρ(j) is j.  
        Subtracting all the possible one-leg contributions from the complete color-kinematic numerator, we obtain a remainder quantity,
        
        
            _f(1… n-2, n-1, n)-∑_j^(1)(1ρ(j)|j, n-1, n)   .
        
        By construction, it should also manifestly gauge invariant, free of any massive poles, and contain n-3 massless legs on the left-hand side (since the pieces we have subtracted contain the complete factorization behavior for those configurations). Next, we consider a massive cut containing two massless legs on the right-hand side
        
            [baseline=([yshift=-0.8ex]current bounding box.center)]
            every node=[font=]  
            mybox = [rectangle, draw=black, fill=red, minimum size=0.3cm]
            mydot = [circle, draw=black, fill=black, minimum size=3pt,inner sep=0pt]
             (a)[mybox];
             [left=0.8cm of a] (a0);
             [right=0.8cm of a] (cut);
                [below=0.2cm of cut] (I)I;
             [below=0.2cm of a0] (n)n;
             [right=2.4cm of a] (a1);
             [right=1.6cm of a] (j0)[mydot];
             [red, thick] ((cut) + (-0.1,-0.1)) – ((cut) + (0.1,0.1));
             [red, thick] ((cut) + (-0.1,0.1)) – ((cut) + (0.1,-0.1)colin);
             [above=0.9cm of j0] (j1);
              [left=0.2cm of j1] (j11)j_1;
             [right=0.2cm of j1] (j12)j_2;
             [below=0.2cm of a1] (n1)n-1;
             [above=0.9cm of a] (b)1     ρ(j_1j_2);
             [above=0.8cm of a0] (b0);
             [right=0.4cm of b0] (b11);
             [right=0.8cm of b11] (b12);
             *(a)–[very thick](a0),(a)–[very thick](a1),(a) – [thick] (b11),(a)–[thick](b12),(j11)–[thick](j0),(j12)–[thick](j0);
              .
        
         Then the remainder  eq. (<ref>) in this case tends to 
        
            (p_Θ_L(j_1) p_j_1) v̅ X(1ρ(j_1j_2)) (p_n1ρ(j_1j_2)-m)  X(j_1 j_2)  u   .
        
        where ρ(j_1j_2) is the complement of j_1,j_2 in the set 1,2,…, n-2. The symbol X(j_1 j_2), whose exact form is unimportant at this point, is a function of the polarizations and momenta of the external gluon legs associated with labelings j_1,j_2. 
        
        Now, subtracting the pieces that contribute to further poles in X(j_1j_2) yields a rank one tensor structure to the right of the propagator (p_n1 ρ(j_1j_2)-m). Indeed, any higher rank tensor had to be generated by a Feynman diagram containing additional fermionic propagators, and the only way of canceling the pole would be by using the Dirac matrix algebra, which in turn lowers the tensor rank. To determine the gauge invariant tensor structure, one can fit the factorization behavior in the heavy-mass regime by parameterizing p_n=m  v, p_n-1 = -m  v - q and taking m→∞. Then the corresponding remainder quantity yields <cit.>
        
            (1ρ(j_1j_2),v)× p_Θ_L(j_1) F_j_1j_2 v   .
        
        Comparing the heavy mass limit of  eq. (<ref>) with  eq. (<ref>), the piece of the numerator that contributes exclusively to the two-particle cut is found to be:
        
        
            ^(2)(1ρ(j_1j_2) |j_1 j_2, n-1, n)= v̅ X(1ρ(j_1j_2)) (p_n1ρ(j_1j_2)-m) (p_Θ_L(j_1) F_j_1j_2) u  p^2_n1ρ(j_1j_2)-m^2  .
        
        Subtracting these new contributions to the color-kinematic numerator results in an object that is free of one- and two-leg poles:
        
        
            _f(1… n-2, n-1, n)-∑_j^(1)(1ρ(j)|j, n-1, n)-∑_j_1 < j_2^(2)(1ρ(j_1j_2)|j_1 j_2, n-1, n)   .
        
        
        The procedure can be repeated recursively for cuts containing any number i<n-3 of legs on the right-hand side:
        
        
            ^(i)(1ρ(τ_R)|τ_R, n-1, n)=(-1)^iv̅ X(1ρ(τ_R)) (p_n1ρ(τ_R)-m)  (p_Θ_L(τ_R) F_τ_R) u  p^2_n1ρ(τ_R)-m^2  .
        
        For the last cut with only leg one on the left, the expression contains an additional term, 
        
        
            ^(n-3)(1|2… n-2, n-1, n)
                   =(-1)^n-3 v̅(p_n F_12… n-2) + 1/4F_1 (p_1  F_2… n-2) p^2_n1-m^2 u 
               =(-1)^n-3 v̅H_12… n-2 p^2_n1-m^2 u   .
        
        
        By repeating our previous argument, we see that subtracting all these contributions, 
        
        
            _f(1… n-2, n-1, n)-^(1)-^(2)-⋯- ^(n-3) ,
        
        leads to an object that contains no poles,  i.e., a polynomial piece. Through the use of dimensional analysis, it becomes apparent that a color-kinematic numerator must include the momenta p_i raised to an overall power of 2 and the polarization vector ε_i raised to an overall power of 3. It has been proven that no polynomial expression, which is manifestly gauge invariant, can satisfy these power requirements. This implies that the remainder must be equal to zero. Therefore, the recursive form presented in eq. (<ref>) is valid, and as a result, the evaluation map in eq. (<ref>) is correct.
        
        
        
        § AN ALTERNATIVE EVALUATION MAP FOR GENERATING SIMPLER COLOR-KINEMATIC NUMERATORS
        
        The previous section's evaluation map did not use the massive on-shell state condition. This section explores using this auxiliary constraint to provide an alternative graphical evaluation map. This map removes expression redundancies and further simplifies the intricate Dirac matrix product in the numerators. We will start by considering numerator examples at four, five, and six points, which will form a basis for universal and more straightforward graphical rules for evaluating numerator expressions. 
        
        
        
         §.§ Four-point color-kinematic numerator
        
        Our first example is the four-point color-kinematic numerator. Here, we get directly from the algebraic construction 
         
            𝒩_f(12,3,4)=    T_(12) =v̅H_12 p^2_41-m^2  u   ,
         
        and an expression, which uses the on-shell condition and H_12=p_4 F_12 + 1 4F_1(p_1  F_2) from the general definition eq. (<ref>). 
        
        
        
         §.§ Five-point color-kinematic numerator
        
        At five-point, the algebraic construction yields:
        
           
            𝒩_f(123,4,5)=    T_(123)-T_(12),(3)-T_(13),(2)
            
                =   v̅[ H_123 p^2_51-m^2  
             -H_12 p^2_51-m^2p_512-m p^2_512-m^2 (p_12 F_3) 
                                 -H_13 p^2_51-m^2p_513-m p^2_513-m^2 (p_1 F_2) ]
             u  
               =v̅[ H_123 p^2_51-m^2  -H_12 p^2_51-m^22(p_12 F_3 p_5)+(p_12 F_3) p_3 p^2_512-m^2
                                     - H_13 p^2_51-m^22(p_13 F_2 p_5)+(p_13 F_2) p_2 p^2_513-m^2 ]
             u   .
         
        In the last step, we use the Clifford algebra and the on-shell condition to reduce the rank of the two tensors to a scalar product and obtain a simpler tensor rank. The  H functions are  
        
            H_12   =p_5 F_12 + 1 4F_1(p_1  F_2)  ,
            
            	H_13   =p_5 F_13 + 1 4F_1(p_1  F_3)  , 
            
            	H_123   =p_5 F_123 + 1 4F_1(p_1  F_23)  .
        
        
        
         §.§ Six-point color-kinematic numerator
        
        Finally, we consider the six-point color-kinematic numerator. Our starting point is the algebraic construction of the numerator, which reads,
        
        
            _f(1234, 5, 6)   = 
                -T_(12),(3),(4)-T_(12),(4),(3)+T_(12),(34)
                -T_(13),(2),(4)
               
                -T_(13),(4),(2)+T_(13),(24)-T_(14),(2),(3)
                -T_(14),(3),(2)
               +T_(14),(23)
                +T_(123),(4)+T_(124),(3)
                +T_(134),(2)-T_(1234)  .
        T_(1234) is unchanged. T_(123),(4),T_(124),(3),T_(134),(2) are simplified as in the five point case. There are new non-trivial cancelations for the terms associated with the evaluations of
        -T_(12),(3),(4)-T_(12),(4),(3)+T_(12),(34), -T_(13),(2),(4)-T_(13),(4),(2)+T_(13),(24) and -T_(14),(2),(3)-T_(14),(3),(2)+T_(14),(23). 
        
        Writing out the first one of these terms, we get:
        
        
            -T_(12),(3),(4)-T_(12),(4),(3)+T_(12),(34)
               = -⟨⟨[baseline=([yshift=-0.8ex]current bounding box.center)]every node=[font=]    
                (a)[myblob];
                 [left=0.8cm of a] (a0)[myblob];
                 [right=0.8cm of a] (a2)[myblob];
                 
                   [right=0.8cm of a2] (a4)[myblob];
                   [above=0.8cm of a] (b1)2;
                    [above=0.8cm of a2] (b2)3;
                    
                     [above=0.8cm of a4] (b4)4;
                     [above=0.8cm of a0] (b0)1;
                    [above=0.8cm of a0](j1);
                 [right=0.2cm of j1](j2);
               
                 [right=1.4cm of j2](j4);
               
                   [right=0.8cm of j4](j6);
               
                  [right=0.0cm of j6](j8);
               	 *(a)–[very thick](a0),(a)–[very thick](a2),(a2)–[very thick](a3), (a3)–[very thick](a4),(a0) – [thick] (j1),(a) – [thick] (j2),(a2) – [thick] (j4),(a4) – [thick] (j8);
                ⟩⟩-⟨⟨[baseline=([yshift=-0.8ex]current bounding box.center)]every node=[font=]    
                (a)[myblob];
                 [left=0.8cm of a] (a0)[myblob];
                 [right=0.8cm of a] (a2)[myblob];
                 
                   [right=0.8cm of a2] (a4)[myblob];
                   [above=0.8cm of a] (b1)2;
                    [above=0.8cm of a2] (b2)4;
                    
                     [above=0.8cm of a4] (b4)3;
                     [above=0.8cm of a0] (b0)1;
                    [above=0.8cm of a0](j1);
                 [right=0.2cm of j1](j2);
               
                 [right=1.4cm of j2](j4);
               
                   [right=0.8cm of j4](j6);
               
                  [right=0.0cm of j6](j8);
               	 *(a)–[very thick](a0),(a)–[very thick](a2),(a2)–[very thick](a3), (a3)–[very thick](a4),(a0) – [thick] (j1),(a) – [thick] (j2),(a2) – [thick] (j4),(a4) – [thick] (j8);
                ⟩⟩+⟨⟨[baseline=([yshift=-0.8ex]current bounding box.center)]every node=[font=]    
                (a)[myblob];
                 [left=0.8cm of a] (a0)[myblob];
                 [right=0.8cm of a] (a2)[myblob];
                 
                 
                   [above=0.8cm of a] (b1)2;
                    [above=0.8cm of a2] (b2)34;
                    
                   
                     [above=0.8cm of a0] (b0)1;
                    [above=0.8cm of a0](j1);
                 [right=0.2cm of j1](j2);
               
                 [right=1.2cm of j2](j4);
                 [right=0.4cm of j4](j5);
               
               
               
               	 *(a)–[very thick](a0),(a)–[very thick](a2), (a0) – [thick] (j1),(a) – [thick] (j2),(a2) – [thick] (j4),(a2) – [thick] (j5);
                ⟩⟩
               =-v̅ H_12 (p_612-m)(p_12 F_3) (p_6123-m)(p_123 F_4) u/(p_61^2-m^2) (p_612^2-m^2) (p_6123^2-m^2)
               -v̅ H_12 (p_612-m)(p_12 F_4) (p_6124-m)(p_12 F_3) u/(p_61^2-m^2) (p_612^2-m^2) (p_6124^2-m^2)
               
            +v̅ H_12 (p_612-m)(p_12 F_3 F_4) u/(p_61^2-m^2) (p_612^2-m^2)  ,
        
        where 
        
            H_12   =p_6 F_12 + 1 4F_1(p_1  F_2)   .
        
        Terms involving two propagators yield
        
        
            v̅ H_12/p_61^2-m^2[(p_12 F_3)(p_123 F_4)/(p_612^2-m^2) +(p_12 F_4)(p_12 F_3)/(p_612^2-m^2)+ (p_612-m)(p_12 F_34)/(p_612^2-m^2) ] u
            
               =   v̅ H_12/p_61^2-m^2[(p_12 F_3)(p_123 F_4)/(p_612^2-m^2)  +(p_12 F_4)(p_12 F_3)/(p_612^2-m^2) 
               +(p_12 F_34) p_34/(p_612^2-m^2)  
                                +2(p_12 F_34 p_612)/(p_612^2-m^2) ] u    .
        
        The terms contained within the square brackets on the right-hand side we rephrase further as follows:
        
        
            2(p_12 F_3 F_4  p_612)/(p_612^2-m^2) +
                (p_12 F_3)(p_3 F_4)/(p_612^2-m^2) 
                +(p_12 F_3 F_4) p_34/(p_612^2-m^2) 
                +(p_12 F_3)(p_12 F_4)/(p_612^2-m^2) 
                +(p_12 F_4)(p_12 F_3)/(p_612^2-m^2)   ,
        
        and reexpress by employing the identity 
        
        
            (p_12 F_3)(p_12 F_4)/(p_612^2-m^2) 
                +(p_12 F_4)(p_12 F_3)/(p_612^2-m^2) =
                -2(p_12 F_3 F_4  p_12)/(p_612^2-m^2)   .
        
        It leads to the following simplified expression:
        
        
            v̅ H_12/p_61^2-m^2[2(p_12 F_3 F_4  p_6)/(p_612^2-m^2) 
                +(p_12 F_3)(p_3 F_4)/(p_612^2-m^2)
                +(p_12 F_3 F_4) p_34/(p_612^2-m^2) ] u   .
        
        From this expression, we make the following observations. It is manifest that, in the simplified form, all the mass-dependent contributions are gone in the numerator, and the massive momentum-dependent terms only appear in the scalar part. Fully accounting for the on-shell conditions allows us a much-simplified expression.
        
          It follows from above examples, that the function H_1τ in the evaluation map of the first component of the algebraic generator is unchanged. Other components are reduced to  
        
            H_jτ(p_ρ,p_ω)   =2 p_ρ F_jτ p_ω+
            ∑_σ_1∈τ(p_ρ  F_jσ_1)( p_Θ'_L(σ_2) F_σ_2)   ,  if  j>1  .
        
        The expression involves the summation of all possible orderings of τ, denoted as σ_1σ_2≡τ, here, Θ'_L(σ_2) represents all the indices in jσ_1 that come before the index of σ_2 in the canonical ordering. When σ_2 is an empty set, Θ'_L includes every index in jτ. In this case, F_σ_2 reduces to an identity matrix. For example 
        
            H_j(p_ρ,p_ω)   =2 p_ρ F_jτ p_ω+ (p_ρ  F_j) p_j ,
            
            	H_jk(p_ρ,p_ω)   =2 p_ρ F_jk p_ω+ (p_ρ  F_jk) p_jk+ (p_ρ  F_j)(p_j  F_k)  .
        
        
         §.§ General color-kinematic numerator
        
        Armed with the understanding from the above examples, we will now introduce a more general and universal formalism. 
        We have the following simplified evaluation map for numerators:
        
        
            ⟨⟨ T_(1τ_1),(τ_2),…, (τ_r)⟩⟩≡⟨⟨[baseline=([yshift=-0.8ex]current bounding box.center)]every node=[font=]    
                (a)[myblob];
                 [left=0.8cm of a] (a0)[myblob];
                 [right=0.8cm of a] (a2)[myblob];
                  [right=0.8cm of a2] (a3)[myblob];
                   [right=0.8cm of a3] (a4)[myblob];
                   [above=0.8cm of a] (b1)τ_1;
                    [above=0.8cm of a2] (b2)τ_2;
                    [above=0.8cm of a3] (b3)⋯;
                     [above=0.8cm of a4] (b4)τ_r;
                     [above=0.8cm of a0] (b0)1;
                    [above=0.8cm of a0](j1);
                 [right=0.2cm of j1](j2);
                 [right=0.6cm of j2](j3);
                 [right=0.4cm of j3](j4);
                 [right=0.8cm of j4](j5);
                   [right=0.2cm of j5](j6);
                 [right=0.6cm of j6](j7);
                  [right=0.0cm of j7](j8);
                 [right=0.8cm of j8](j9);
               	 *(a)–[very thick](a0),(a)–[very thick](a2),(a2)–[very thick](a3), (a3)–[very thick](a4),(a0) – [thick] (j1),(a) – [thick] (j2),(a)–[thick](j3),(a2) – [thick] (j4),(a2)–[thick](j5),(a4) – [thick] (j8),(a4)–[thick](j9);
                ⟩⟩
               =v̅H_1τ_1 p^2_n1-m^2H_τ_2(p_Θ_L(τ_2),p_Θ_R(τ_2)) p^2_n1τ_1-m^2 ⋯H_τ_r(p_Θ_L(τ_r),p_Θ_R(τ_r)) p^2_n1τ_1…τ_r-1-m^2 u  .
        
        We obtain this map by analyzing the color-kinematic numerator pattern observed in the above examples provided by eq. (<ref>) using the additional on-shell conditions. 
        We will also use the convenient musical diagrams <cit.> to ease the computation of numerators. The musical diagrams  take the form
        
        
            [baseline=([yshift=-0.8ex]current bounding box.center)]every node=[font=]    
                (l1)[]; 
                 [below=0.6cm of l1](ls)[]; 
                 [above=0.9cm of l1](lg)[]⋮; 
                  [right=1.0cm of lg](lg3)[]⋮; 
                 [right=3.5cm of lg](lg2)[]⋮; 
                 [left=0.5cm of ls](rms)[](ητ_1);
                 [right=5.cm of ls](rs); 
                 [right=4.3cm of ls](vphi)[myblob2]𝐧; 
                 [left=0.5cm of l1](rm1)[](τ_2);
                 [right=5.cm of l1](r1)[]; 
                 [right=0.5cm of ls](v1)[sb]1;
                  [right=1.9cm of ls](v11)[sb];
                  [right=2.9cm of ls](v12)[sb];
                 [right=1.cm of l1](v3)[sb];
                 [right=3.5cm of l1](v4)[sb];
               	 *(l1)- -[thick](v3)- -[scalar,thick](v4)–[thick](r1),(ls)–[thick](v1)–[thick](v11)–[scalar,thick](v12)–[thick](vphi)–[thick](rs);
        
        In a given diagram, the reference leg indices are located at the bottom line, and we sum all possible positions of the non-reference legs. For a given diagram, each horizontal line contains an argument of an  H function required for the corresponding numerator. (We use the notation Θ_L/R(τ_i) that refers to the set of all the left-lower/right-lower indices of τ_i. 
        This definition of Θ_L(τ_i) is consistent with the Θ_L function in the definition of gravity KLT momentum kernel used in eq. (<ref>)). 
        
        Below, we illustrate using musical diagrams to evaluate abstract generator terms. As an example, we can consider ⟨⟨ T_(16),(235),(4)⟩⟩, which appears in the eight-point color-kinematic numerator.
        
        
            [baseline=([yshift=-0.8ex]current bounding box.center)]every node=[font=]    
                (l1)[]; 
                 [below=0.4cm of l1](ls)[]; 
                 [above=0.2cm of l1](ls2)[]; 
                 [left=0.5cm of ls2](rms2)[](τ_2); 
                 [right=1.cm of ls2](v22)[sb]2; 
                 [right=1.8cm of ls2](v23)[sb]3; 
                 [right=3.2cm of ls2](v25)[sb]5; 
                 [right=5.0cm of ls2](rs2); 
                 [above=0.8cm of l1](ls3)[]; 
                  [left=0.5cm of ls3](rms3)[](τ_3); 
                  [right=2.5cm of ls3](v34)[sb]4; 
                  [right=5.0cm of ls3](rs3); 
                 [left=0.5cm of ls](rms)[](ητ_1); 
                 [right=5.0cm of ls](rs); 
                 [right=4.5cm of ls](vphi)[myblob2]8; 
               
                 [right=5.0cm of l1](r1)[]; 
                 [right=0.5cm of ls](v1)[sb]1; 
                 [right=3.8cm of ls](v4)[sb]6; 
               	 *
               	 (ls)–[thick](v1)–[thick](v4)–[thick](vphi)–[thick](rs),(ls2)–[thick](v22)–[thick](v23)–[thick](v25)–[thick](rs2),(ls3)–[thick](v34)–[thick](rs3);
        
        Reading the arguments of the  H functions (by following the indices of the various horizontal lines), we arrive at the contributions
        
            ⟨⟨ T_(16),(235),(4)⟩⟩   =
              v̅H_16 p^2_81-m^2 H_235(p_1,p_68) p^2_816-m^2H_4(p_123,p_568) p^2_816235-m^2 u  ,
        
        where the H_235 function reads:  
        
        
            H_235(p_1,p_68)   =2(p_1 F_235 p_68)+(p_1 F_25) (p_2 F_3)
               +(p_1 F_235) p_235
            +(p_1 F_23) (p_23 F_5)+(p_1 F_2) (p_2 F_35)   ,
        
        as well as the H_16 function,
        
        
            H_16= (p_8 F_16) + 1 4F_1(p_1  F_6)   ,
        
        and H_4,
        
        
            H_4(p_123,p_568) =2 (p_123 F_4 p_568)+
             (p_123  F_4) p_4  .
        
        
        We conclude this section by presenting the following compact expression for numerators based on the above rules for evaluation: 
        
        
            _f(12… n-2, n-1, n)=(-1)^n-3 v̅[  H_1… n-2 p_n1^2-m^2] u 
               +∑_τ_R⊂{2,…,n-2}(-1)^|τ_R| v̅[
                X(1ρ(τ_R))  H_τ_R(p_Θ_L(τ_R),p_Θ_R(τ_R)) p^2_n1τ_L-m^2] u  .
        
        We have explicitly checked the compact numerator expressions provided by this formula up to eight points. We also note that in all cases the fermion color-kinematic numerator tends to the scalar Yang-Mills case in <cit.> after the replacement v̅γ^μ u→ p_n^μ and v̅γ^μ_1γ^μ_2⋯γ^μ_r u→ 0 as expected for consistency.
        
        
        
        § CONCLUSION AND OUTLOOK
        
        Using a universal method, we have presented a new and compact construction for a color-kinematic numerator associated with the generalized massive Compton amplitude with a fermion pair. We have verified evaluation maps explicitly until eight points and also presented a more straightforward graphical method for computing such numerators.
        
        As the formulation relies on the Lorentz group representations of the matter fields, one can apply similar techniques to miscellaneous Compton scattering amplitudes, for instance, those with two massive vectors or other massive particles. An exciting track is to investigate if it is possible to shed new light on gauge and gravity interactions in high-spin matter fields at finite spin order. Studying the classical limit for such amplitudes and resolving questions related to consistent factorization limits could potentially aid the development of more efficient computational methods for observing spinning black hole mergers. It is beyond this paper's scope, so we leave these and related questions for future studies.
        
        
        
        § ACKNOWLEDGEMENTS
        
        The work of N.E.J.B.-B. and G.C. was supported by the DFF grant 1026-00077B and partially by the Carlsberg Foundation. G.C. has received funding from the Marie Sklodowska-Curie grant agreement No. 847523 "INTERACTIONS" under the European Union Horizon 2020 research and innovation program.
        
        JHEP
