Two-Loop anomalous dimensions of QCD operators up to dimension-sixteen and Higgs EFT amplitudes

We consider two-loop renormalization of high-dimensional Lorentz scalar operators in the gluonic sector of QCD. These operators appear also in the Higgs effective theory obtained by integrating out the top quark loop in the gluon fusion process. We first discuss the classification of operators and how to construct a good set of basis using both off-shell field theory method and on-shell form factor formalism. To study loop corrections, we apply efficient unitarity-IBP strategy and compute the two-loop minimal form factors of length-3 operators up to dimension sixteen. From the UV divergences of form factor results, we extract the renormalization matrices and analyze the operator mixing behavior in detail. The form factors we compute are also equivalent to Higgs plus three-gluon amplitudes that capture high-order top mass corrections in Higgs EFT. We obtain the analytic finite remainder functions which exhibit several universal transcendentality structures.


Introduction
Understanding hadron spectrum is one central problem in QCD. The hadron states can be understood as gauge invariant composite operators which are composed of gluon and quark fields. The problem of computing the hadron spectrum is then reduced to calculating the (anomalous) dimension of composite operators. Due the non-perturbative nature of confinement, an analytic derivation of anomalous dimensions remains a dream; 1 on the other hand, at high energy scale the asymptotic freedom ensures that a perturbative expansion still applies. A good knowledge of such perturbative information is helpful to understand the RG flow of the spectrum, and should also provide an important probe to the full spectrum. One goal of this paper is to provide a working framework that can be efficiently used to compute the anomalous dimension of high-dimensional operators as well as at high loop orders. To be concrete, we will focus on gauge invariant and Lorentz invariant local operators O(x) where all elementary fields are located at a common point in spacetime.
As another motivation, the local operators we consider are also related to the Higgs effective action, which describes the Higgs production via gluon fusion process at LHC. The Higgs particle has no direct interaction with gluons but through Yukawa coupling with quarks. The coupling is proportional to the mass of quarks, which is dominated by the heaviest top quark [2,3]. To simplify the computation, a useful approximation is to use an effective field theory (EFT) which describes the interaction between Higgs and gluons by integrating out heavy top quark [4][5][6][7][8][9][10]. The EFT Lagrangian can be schematically given as: whereĈ i is the Wilson coefficient, H is the Higgs field, and O ∆ 0 ;i are the effective operators of canonical dimension ∆ 0 . For the Higgs plus one jet production, the contribution of higher dimension operators can be important when the Higgs transverse momentum is comparable to the top mass. The two-loop Higgs plus three-parton amplitudes with the leading operator O 4;0 = Tr(F µν F µν ) were computed in [11], and similar two-loop amplitudes with dimension-6 operators were computed in [12,13]. The two-loop amplitudes with higher dimension operators may be used to improve the precision for the cross section of Higgs plus a jet production at N 2 LO, which is so far known in the infinite top mass limit [14][15][16][17][18][19][20]. At NLO QCD accuracy, the full top mass effect can be taken into account by integrating the top quark loop directly [21][22][23]. See also [24] for a recent extensive review about related studies on Higgs amplitudes and their phenomenological applications.
To study the operator spectrum and the corresponding Higgs amplitudes, we consider the form factor which is defined as a matrix element between an operator O(x) and n on-shell states (see e.g. [25] for an introduction): F O,n = d 4 x e −iq·x p 1 , . . . , p n |O(x)|0 . (

1.2)
Such form factor is equivalent to a Higgs plus n-parton amplitude in the Higgs EFT (1.1), where q 2 = m 2 H . In this following, we will often refer Higgs amplitudes as form factors.
The concrete goal of this paper is to clarify the operators of general high dimensions and to compute their anomalous dimensions and related Higgs amplitudes up to two loops. Since we discuss multiple high dimension operators, the first problem to solve is to classify the operators and construct a convenient operator basis (at classical level). The operators at a given canonical dimension are generally not independent with each other, since they may be related to each other through equations of motion or Bianchi identities. We will apply both (off-shell) field theory method and (on-shell) minimal form factor method to construct the operator basis. We point out that rather than counting the dimension of the basis in this paper, a specific emphasis is put on how to construct a good set of basis operators explicitly that will be convenient for the high loop computations.
To compute the loop form factors for the basis operators, we apply the unitarity-IBP strategy that combines unitarity cut [26][27][28] and integration by parts (IBP) methods [29,30]. This strategy has been applied to computing form factors (and Higgs amplitudes) in [12,13,31] and for pure gluon amplitudes in [32,33]. Similar strategy has also been used in the numerical unitarity approach [34,35], in which unitarity and IBP are used together with numerical variables to avoid large intermediate expressions. The idea of using cuts to simplify IBP has also been used in e.g. [36][37][38][39].
Given the form factor results, one can extract the ultraviolet (UV) divergences by subtracting the universal infrared (IR) divergences. The basis operators of same classical dimensions can mix with each other through renormalization at quantum level (see e.g. [40]). This is captured by the renormalization matrix Z defined in O R,i = Z j i O B,j . With the explicit two-loop results, we can study various aspects of the operator mixing in detail. By diagonalizing the renormalization matrix, one can also compute anomalous dimensions.
As mentioned above, the form factors we consider also correspond to the Higgs amplitudes. In particular, high dimension operators corresponds to the high-order top mass corrections in (1.1). The central quantity of the Higgs amplitudes is the finite remainder function. We obtain the analytic expressions, and find that that maximal transcendental part is universal, namely, independent of the operators. This provides further support of the maximal transcendentality principle introduced in study of anomalous dimensions in [41,42] and then observed for the form factor of tr(F 2 ) in [43] and later also for other operators [12,13,31,[44][45][46][47][48][49][50][51][52]. Besides, the lower transcendentality degree-3 and degree-2 parts also contain universal building blocks, which was observed in the results with dimension-6 operators [13].
This paper is structured as follows. In Section 2, we construct the operator basis using both the field theory method and on-shell minimal form factor methods. Section 3 describes the strategy of computing form factors based on the unitarity-IBP method. In Section 4, the renormalization of the form factors are discussed in detail, from which we extract the renormalization matrix of the operators and compute anomalous dimensions. In Section 5, we consider the full Higgs amplitudes and in particular we focus on the finite remainder functions and discuss their properties. Section 6 provides a summary and outlook. Several technical details are included in the Appendix A-E. For convenience, we also summarize the structure of the paper as following graph: In this section we consider the construction of operator basis. We will first consider the field theory method and then apply the on-shell form factor method. Besides counting the number of basis, a central goal is to explain how to construct a convenient set of basis operators that will facilitate the high loop computations. We will provide explicit basis for length-3 operators up to dimension 16, and in later sections we will compute their anomalous dimension and related Higgs EFT amplitudes.

Operator setup
We consider local gauge invariant scalar operators in pure Yang-Mills theory composed of field strength F µν and covariant derivatives D µ . The field strength carries an adjoint color index as F µν = F a µν T a , where T a are the adjoint generators of gauge group and satisfy The covariant derivative acts in the standard way as

2)
A gauge invariant scalar operator takes the following form c(a 1 , ..., a n ) D µ 11 ...D µ 1m 1 F ν 1 ρ 1 a 1 · · · D µ n1 ...D µnm n F νnρn an X(η, ǫ) , (2.3) where c(a 1 , ..., a n ) are color factors, such as given in terms of products of Tr(..T a i ..T a j ..). And to form a scalar operator, all Lorentz indices {µ i , ν i , ρ i } are contracted in pairs by metric η µν or by antisymmetric tensor ǫ µνρσ , which are contained in the function X(η, ǫ). In this paper, for simplicity we will consider the parity even operators where X contains only η's. For convenience of the upcoming discussions, we define following useful quantities for the operators: (2.4) Since we consider Lorentz scalar operators, this dimension is always an even integer number, starting with dim=4. The canonical dimension typically receives quantum corrections at loop level, and the correction is called the anomalous dimension γ(O).
• Length of an operator: len(O) = (# of F 's) . (2.5) We will also call an F in an operator together with all the derivatives in front of it, i.e. (D . . . DF ), as one block. Obviously, the number of blocks is equal to the length of the operator.
• Descendants. If an operator can be written as a total derivative of lower dimensional operator, it is called a descendant. Since we always consider Lorentz invariant and gauge invariant operators, the overall derivatives are just covariant derivatives and always appear in pairs so a descendant must take the form like Below we also summarize the color factors for length-2 to length-4 operators: 1. The length-2 case has a unique color factor: c(a, b) = Tr(T a T b ) = δ ab . Tr(T a 1 T a σ(2) T a σ(3) T a σ(4) ), Tr(T a 1 T aσ (2) )Tr(T aσ (3) T aσ (4) ) , σ ∈ S 3 ,σ ∈ Z 3 .
Our goal is to find a set of independent operator basis, which requires no equivalent relation holds among the basis we choose. Two operators are said to be equivalent if their difference is proportional to the equation of motion (EoM) or Bianchi identity (BI): BI : Below we will use two different ways to do this classification: (1) field theory method, and (2) on-shell spinor helicity method. For the convenience of notation, we will often use integer numbers to represent Lorentz indices and abbreviate product D i D j D k .... to D ijk... . For example, (2.12)

Field theory method
The choice of the operator basis is not unique. To make a preference for the choice, we impose two requirements on basis operators: • We consider operators with smallest length first and then increase the length one by one. If two operators at length l differ by an operator of high length (L > l): we will say that O l and O ′ l are equivalent at length l level, and only one of them may be kept in the operator basis. In the following, we will start with the length-2 operators and then consider length-3 operators and finally briefly discuss length-4 operator.
• At a given dimension, we prefer to choose descendants as basis operators, so that the total derivatives of lower dimensional operators will contribute automatically as the basis operators in higher dimensional cases. Since the total derivatives do not change the anomalous dimension of the operator, this choice has an important advantage, in the sense that the renormalization property of a descendant is determined by its lower dimensional counterpart. We point out that the descendants can have non-trivial mixing with other operators, therefore it is important to keep them in the basis in our computation, see Section 4 for detailed examples.
At a given length and a given dimension, the numbers of D and F are fixed, and different operators are differed by their ways of Lorentz contraction. A covariant derivatives D contracts either with another D or with an F . We call operators with no DD contraction as primitive operators. The full set of operator can be generated by adding pairs of identical Ds to the primitive operators.
Denote the number of Ds in an primitive operator asn d . At a given length l, it is easy to seen d ranges from 0 to 2l and must be even. Starting from arbitrary operator, one can get a primitive one by taking off all the identical DD pairs from it, like: We say the operator O 1 belongs to primitive class ⌊O 0 ⌋.

Length-2 case
In the length-2 case we illustrate the above classification in detail. Since there is only one color structure (2.7), we will ignore the color factor in the discussion for simplicity. First let us classify primitive operators. Number of Ds in a primitive operator might take valuesn d = 0, 2, 4, respectively corresponding to (2.14) Here we have used the EoM constraint (2.10), such that D i should not contract with F ij that contains the same index. The aboven d = 2 operator can be reduced to a non-primitive operator using Bianchi identity: Then d = 4 operator turns out to be a higher length-3 one and therefore can be dropped out at length-2 level: In summary, F 12 F 12 is the only primitive operator for length-2 cases. Second let us construct non-primitive operators through DD pair insertion. One can check that a DD pair acting on the same F is always equivalent to a higher length operator: According to our requirement mentioned around (2.13), at a given length, we only need to consider the cases that pairs of identical D i s being inserted into different blocks. Thus the non-primitive length-2 operator obtained from DD pair insertion must be in the form of Finally we can see such operators are equivalent to descendants: Recursively, an arbitrary length-2 operator, say with dimension ∆ 0 , is equivalent to a total derivative of F 12 F 12 like

Length-3 case
For the length-3 case, the classification of primitive operators is similar to the length-2 case, and we provide the derivation in Appendix A. We directly cite the result here: there are two length-3 primitive operators They can be rewritten as Second, based on that, construction of basis operators at a given dimension ∆ 0 is just to correctly count inequivalent ways of DD pair insertion into blocks of O P1 and O P2 , which require n = ∆ 0 −8 2 and m = ∆ 0 −6 2 pairs of identical Ds respectively. Let us discuss the two cases in more detail.
• We first consider inserting pairs of identical Ds into primitive operator O P1 , which has blocks 1,2,3, namely D 1 F 23 , D 4 F 23 , F 14 . As explained above, D 2 acting on a single F vanishes up to higher length, so a pair of identical Ds can only be inserted into two different blocks, which means there are three choices: 1,2 or 1,3 or 2,3.
Taking n = 2 pairs of Ds as an example, all the six possible insertions are listed below: where the three slots represent the positions of block D 1 F 23 , D 4 F 23 , F 14 . Alternatively, they can be obtained from the three-part partitions of number n = 2: 0 + 2 + 0, 2 + 0 + 0, 0 + 0 + 2, 1 + 1 + 0, 0 + 1 + 1, 1 + 0 + 1 . (2.20) The three parts in a sum correspond to numbers of D-pairs inserted into three block pairs (1,2) + (1,3) + (2,3). For a general n, the number of inequivalent insertions are given by partition counting • The O P2 case is a little bit different, since the three blocks F 12 , F 13 , F 23 have Z 3 symmetry under index renaming. As a result, three insertions given by the first line of (2.19) actually produce the same operator. The same is true for the second line, so there are only two inequivalent ways to insert m = 2 pairs of identical Ds into three blocks of O P2 .
Generally, we should equalize different D-insertions if they are related by Z 3 cyclic symmetry. A special insertion appears when number m is multiple of 3, because the trisection partition produces a D-insertion which is invariant under Z 3 symmetry. In summary, for m = 0 (mod 3), number of inequivalent D-insertions equals to one third of the partition number, i.e., 1 3
Third, operators listed in Appendix B are not the final basis choice. Since we prefer to choose descendants as the basis elements, we need to further solve descendant relations and replace some operators by descendants.
Finally, it is convenient to symmetrize the operators according to the color factors. As discussed in (2.8) -(2.9), for length-3 operators, one can take either single trace products or {f abc , d abc } as color basis. In this paper we choose the later. To get an operator with color factor d abc or f abc from a single traced one, we need to (anti)symmetrize the kinematic part.
To illustrate the above procedure, let us take the dimension 10 case as a concrete example. To start with, there are five length-3 operators obtained by add DD pairs in the two primary classes: As mentioned before, we choose length-3 color factor to be {f abc , d abc } instead of single trace basis {Tr(T a T b T c ), Tr(T a T c T b )}. To achieve that we need to symmetrize or antisymmetrize the original ones (2.22), and the resultant new basis are: (2.23) Here equivalence relation A ∼ B means the difference between A and B is a higher length operator. Furthermore, we prefer to choose descendants as basis operators. At dimension 10, there are two descendants from total derivatives of lower dimensional operators O P1 and O P2 , and they are related to operators in (2.23) as:

Length-4 case
Above construction should be straightforward to generalize to higher length operators. Here we will not go into details but only give a brief account of the length-4 case. For length-4 case there are seven types of primitive operators as shown in Table 1.
Here we provide the basis of dimension-8 operators as a concrete example. At dimension 8, all length-4 operators are primitive ones, and they correspond to the first two rows of Adding length-2 and length-3 operators, at dimension-8, the operator basis contains 1+2+8 = 11 operators in total. Unlike length-2 and length-3 cases however, length-4 operators from different primitive classes might not be independent, so naively adding pairs of D into the primitive operators may create an overcomplete basis set. For example, consider an identity involving three operators: These operators belong to three different primitive classes at the fifth row of Table 1, i.e., (abdc), (dbca), (abcd). Therefore the linear independence of primitive operators doesn't guarantee the independence of non-primitive ones.
Before entering next subsection, let us summarize the strategy of field theory classification: 2. After primitive operators being classified, one can then generate other (non-primitive) operators by enumerating inequivalent ways of DD pair insertion into primitive ones.
3. While independent operators obtained from inserting DD pairs into primitive ones already form a set of basis, they are not a good choice since we require descendants to be included. One can apply identities between descendants and above basis and solve for part of them in terms of descendants.
4. For length-3 case, we also organize the basis into f abc and d abc sectors, to manifest the symmetry properties.
5. Finally, to obtain a full basis at a given dimension, one needs to sum operators of all possible length. For example, for dimension-8 case, operators up to length-4 all contribute.

On-shell spinor helicity method
An alternative way to do the classification is to make use of on-shell technique and read properties of operators from their form factors. In other words, we follow the dictionary between operators and their tree-level minimal form factors and translate all the operator information into spinor helicity formalism for mathematical clearness.

Operator-spinor dictionary
A minimal tree form factor means the number of external gluons is equal to the length of the operator. One can establish a dictionary from an operator to its tree-level minimal form factor [53][54][55]: especially each single D and F contained by the operator are mapped to certain spinor structures, as shown in Table 2.
The map of D results from spinor represention of momentum Pα α =λα p λ α p . As for field strength F , first one takes decomposition F µν → F ααββ = ǫ αβfαβ + ǫαβf αβ (2.28) to obtain self-dual and anti-self-dual components Then one makes use of LSZ reduction formula to get their final matrix elements We summarize the correspondence between operators and on-shell spinors in Table 2, and the example on reconstructing operators from spinor-helicity formalism will be given in upcoming context, see (2.42). The correspondence listed in Table 2 is not limited within pure Yang-Mills theory, and the result can be generalized when fermions enter in. The above on-shell language has several advantages: 1. Equivalent relations between operators take much simpler forms. Equation of motion holds automatically, and Bianchi identities are translated into Schouten identities: 2. Contribution from higher length operators vanish automatically since they have vanishing form factor now. Therefore two operators equivalent up to higher length have identical tree level minimal form factor. For example, following three operators are equivalent at the level of length 2: and they have the same form factor under arbitrary helicity setting, like s 12 12 2 for 1 − 2 − and 0 for 1 − 2 + .
3. In field theory classification we treat DF contraction and DD contraction differently. In on-shell language, DD contraction only contributes to scalar factor like s ij , so the primitive class just corresponds to minimal form factors with all the s ij stripped off: primitive class → ignoring scalar factors .
For example: Below we will discuss the construction of operator basis based on the spinor-helicity formalism. As we will see, all the standard steps in field theory method have on-shell analogies. Primitive operator classifiction becomes enumerating spinor structures, and counting DD pair insertion becomes enumerating scalar factors. Besides, the minimal form factor picture also allows a further choice of the helicity sector which will further simplify the kinematic structure.

Length 2
For length-2 operators, the tree-level minimal form factor involves two scattered gluons labeled by 1 and 2. The only possible spinor bracket structures are 12 and [12], since forms like 1|P i |2] must vanish because P i is either p 1 or p 2 .
For an operator with minimal dimension, the only possible expressions of the form factor are 12 2 and [12] 2 , corresponding to projected operator tr(f αβ f αβ ) and tr(fαβfαβ). Their sum is tr(F 2 ), and difference is ǫ µνρσ tr(F µν F ρσ ). In this paper we will not consider operators with odd P -parity, so only the former is kept.
For an even operator with general dimension ∆, its minimal form factor is (s 12 ) . Taking dimension 6 as an example, from s 12 12 2 one can read two holomorphic operators, which are self-dual components of the first two operators in (2.32), and they are both equivalent to the descendant 1 4 ∂ 2 tr(F 2 ), which is chosen as the only independent length-2 operator at dimension 6: (2.33)

Length 3
As for length-3 cases, we introduce a new concept helicity sector to help the discussion. The map from operator equivalent class to tree-level minimal form factor is not injective. For example, for two inequivalent operators O ′′ 8;1 and O ′′ 8;2 given in (B.2):  Table 3. One can see these two operators are undistinguishable under helicity (−, −, −).
A natural solution is to create two new operators using O ′′ 8;1 and O ′′ 8;2 so that one has nonzero tree-level minimal form factor only under (−, −, +) and its h.c., and the other is non-zero only under (−, −, −) and its h.c. We say these two new operators belong to helicity sector α and β respectively: α-sector : F (2. 35) In following context, we require each length-3 basis operator belongs to either α-sector or β-sector.
The nature of "helicity sector" for an operator is its holomorphic structure. Taking decomposition (2.28) of a length-3 operator might create four possible components f f f ,fff , f ff,ff f . For operator with even P -parity, conjugate components always appear in pairs so there are only two inequivalent structures: We can see the former can only be detected by helicity (−, −, +) and (+, +, −), while the later only by helicity (−, −, −) and (+, +, +). Let us redo the classification for length-3. First we need to enumerate all the possible spinor structures for 3-gluon form factors, which is an analogy to enumerating primitive operators in field theory method. There are two types of spinor structures The second step is to add Mandelstam scalar factors to A 1 and A 2 until certain dimension is reached, and this is an analogy to inserting DD pairs to primitive operators in field theory method. Below we take dimension 10 case as an example to explain the details.
1. Let us enumerate scalar factors for spinor structures A 1 and A 2 . For A 1 , dimension of scalar factor is ∆ − 8, so at ∆ = 10 there are three choices: s 12 , s 13 , s 23 . Notice that helicity setting −, −, + has broken the total bosonic symmetry, while the exchange invariance between gluon 1 and 2 is still maintained. So 1,2-flipping property of kinematic part should be compatible with that of color factor. As a result, s 12 and s 13 + s 23 correspond to f -sector while s 13 − s 23 corresponds to d-sector.
For A 2 , scalar factor has dimension ∆ − 6 and must be cyclic symmetric as A 2 . At ∆ = 10 there are two choices: s 2 12 + s 2 23 + s 2 13 and s 12 s 23 + s 12 s 13 + s 23 s 13 , which both belong to f -sector.
(2.40) 2. After obtaining above spinor-helicity forms, one can apply the dictionary in Table 2 to read out the operators. Taking s 12 A 1 as an example, we write the bracket form back to spinors. For an operator of primitive class ⌊O P1 ⌋ under helicity setting (1 − 2 − 3 + ), one can lock the external gluon label i to block position i, so we have is not the only way to group spinors into counterparts of f and D. For example, (2.41) can also be rearranged to which is the f ff-component of operator Tr(D 12 F 34 D 13 F 45 F 25 ) equivalent to − 1 2 O ′ 10;1 : From above example we can see that the reconstruction of operators from on shell minimal form factor is not an injective mapping, but the different operators obtained from the same form factor are equivalent (up to possible overall factor) at the level of length-3, so one can choose an arbitrary operator as the representative one and here we choose (2.42) for s 12 A 1 .
For an A 2 -involved form factor we do the similar thing. A little bit difference comes from the fact that under helicity (−, −, −) the form factor has full bosonic symmetry so gluon i might come from any one of the three blocks. So first we write the cyclic symmetric scalar back to an asymmetric one (but maintaining (1, 2)-flipping symmetry), e.g. write s 12 s 23 + s 13 s 23 + s 12 s 13 back to s 13 s 23 , (s 12 + s 13 )s 23 , (s 12 − s 13 )s 23 , and then transform the spinor into operator according to dictionary, with gluon label i locked to Table 4. Operator-form factor dictionary of basis operators at dimension 10. minimal form factor projected operator full operator   Table 4 is mapped to three inequivalent operators.

Basis operator
In total, there are five different operators in the third column, they are equivalent to the basis operators obtained from field theory method introduced in Section 2.2, see (2.23).
3. Now we need to solve the descendant relations as constraints. One can easily find operator relations (2.24) between {O ′ 10;i } and descendants of O P1 and O P2 from scalar factor relations s 123 = s 12 + s 13 + s 23 . We choose to replace O ′ 10;2 and O ′ 10;4 with these two descendants and the basis set becomes (2.25).
4. Finally we need to recombine five operators in (2.25) so that each basis operator only belongs to one helicity sector. The final operator basis is summarized in Table 5.
Let us explain the operator notation in Table 5. In following context, we will label length-3 basis operators as O ∆;α/β;f /d;i , where ∆ is the operator dimension, α/β labels the helicity sector, f /d labels the color factor, and i is the numbering. We also give the final basis operators for dimension 6 and 8 in Table 6, and final basis for dim 12-16 are given in Table 6. Final basis operators at dimension 6 and 8 cases. Table 7. Color-ordered tree-level minimal form factors of length-4 operators at dimension 8. The color factors are tr(1234) and tr(12)tr(34) for single and double trace operators, respectively.

Basis operator
Appendix C. These operators will be the starting point for the computation of subsequent sections.

length 4
For length-4 case, we take dimension 8 as an example. There are eight independent operators in total as shown in (2.26). When considering color ordered form factors, we take coefficients of tr(1234) for single trace operators and take coefficients of tr(12)tr(34) for double trace operators. In order to make each operator belong to only one helicity sector, i.e. having color-ordered form factors that are non-zero under only one type of helicity configuration, we recombine above eight operators as follows: The correspondence between these newly obtained operators and their color-ordered form factors are summarized in Table 7.
Notice the last column of Table 7 corresponds to uniform helicity, so a total form factor has full bosonic symmetry, which means its kinematic part and color factor have the same symmetric property under label permutation. For example, kinematic part of single-trace operators are invariant under Z 4 cyclic permutation, and those of double trace operators are invariant under 1 ↔ 2, 3 ↔ 4 permutation.
Form factors under one plus three minus or three plus one minus are not written in the table. The results are all zero, because these eight operators do not have f f ff orfff f components under decomposition (2.28).
For lengh-2 and length-3 operators, the countings of operators in these two different approaches agree with each other. However, this is not the case for higher length operators, where the evanescent operators appear. Such evanescent operators do not show up in the 4-dim spinor approach, on the other hand the field theory approach is valid for generic dimensions and can captures these operators. We leave the discussion to the future work.

Two-loop form factor computation via unitarity
In this section, we compute the one and two-loop form factors of the high dimensional operators discussed in the last section. Our computation is based on the on-shell unitarity methods [26][27][28], where the cut integrands are constructed by sewing tree-level components. Furthermore, we combine the unitary method together with the integration by parts (IBP) reduction [29,30]. This "unitarity-IBP" strategy not only makes the computation very efficient, but also provides important internal consistency checks for the results. Below we first outline the main strategy of the computation and then apply it to the concrete form factor computations.
The work flow of our calculation can be illustrated as follows: where I i are IBP master integrals. In the beginning, a particular cut channel (or cut configuration) is chosen and one can calculate the cut integrand through tree-level data. In order to avoid the issue of rational terms, here it is essential to use D-dimensional cut instead of four-dimensional cut. The resulting cut integrand contains all the integrals whose topologies are permitted by the chosen cut. As for the integral reduction, we use IBP method combined with on-shell conditions for the cut propagators. Because terms proportional to cut inverse propagators vanish under cut condition, the expressions of IBP relations can be sharply shortened and therefore computing efficiency is improved. After cut-constrained IBP reduction, one obtains coefficients c i of all the cut-permitted master integrals. Finally, by repeating the process for different cut channels, coefficients of all the master integrals are probed. See also [31] for discussion. In this paper we will mostly focus on the three-point form factors of length-three operators up to two-loop level. In these cases only planar integrals appear (and therefore the planar cuts are sufficient). This can be understood from a simple color analysis of Feynman diagrams that up to two-loop level the minimal form factors of length-3 operators do not contain subleading (1) (2) Figure 1. (1) The 2-loop non-planar topology has vanishing color factor. (2) Nonplanar topology contributing to leading color begins to appear at 3-loop.
(1) color factors. The only two-loop non-planar topology is shown in Fig. 1(a), but its color factor is zero, for length-3 operators in both f abc -sector and d abc -sector. As a comparison, at three loops nonplanar topology (even at leading N c color) will contribute, as shown in Fig. 1(b), and therefore nonplanar cut is necessary. Since the one-loop case is quite simple, below we will focus on the two-loop computation.
The complete set of two-loop master integrals for minimal length-3 form factors are given in Fig. 2. With color decomposition, the two-loop color-ordered form factors, associated with color factor tr(T a 1 T a 2 T a 3 ), can be written as a sum of master integrals I i as where master integrals {I i } strictly correspond to the topology and labeling given in Fig. 2.
The master coefficients c i are what we want to obtain using unitarity-IBP method. Before considering that, let us discuss one important feature of the master integrals.
One can see that (5), (6) and (7), (8) in Fig. 2 are pairs of 'mirror' topologies. In color-ordered form factors, they should be considered to be independent because they are inequivalent planar diagrams and therefore probed by different planar cuts. On the other hand, they are closely related to each other: graphically, (5) and (6) are related by label flipping 1 ↔ 2, while (7) and (8) are related by flipping 3 ↔ 1. From the planar color point of view, they are related by reversing color orientation, which is equivalent to a "C-parity transformation" (see e.g. [56]), so the kinematic parts of a fixed color order tr(123) and the Figure 3. Complete set of cuts fully probing contributions from all the master integrals  (5) and (6) as well as (7) and (8) are related with each other as: Notice also that I 3 and its two cyclic partners share a degenerate expression, but here we treat them as distinct ones and sum cyclic permutations together. A spanning set of planar cuts fully probing these master integrals are shown in Fig. 3. As already mentioned, a particular cut can probe only a subset of master integrals. Among master integral coefficients, c 1 , c 2 , c 3 , c 4 are probed respectively by cuts (c), (b), (a), (d) in Fig. 3. To probe c 5 one should apply s 123 -triple-cut (a), which also probes the coefficient of integral I 6 | (p 3 →p 1 →p 2 →p 3 ) . To probe c 7 and c 9 one can apply s 12 -triple-cut (b), or s 312triple-cut. Notice the coefficients of I 8 | (1→3→2→1) and I 9 | (1→2→3→1) can also be probed by cut (b). Since different cut channels can probe same or symmetry-related master integrals, this provides strong consistency checks for the results.
Below we provide some more details of the calculation by considering a particular cut channel. Taking cut (b) in Fig. 3 as an example, this cut allows us to determine the coefficients of master integrals as shown in Fig. 4.
The cut integrand is obtained by sewing a planar four-gluon tree form factor together with a planar five-gluon tree amplitude. Since we consider D-dimensional cuts, the tree results are computed via Feynman rules. The sewing process involves the helicity sum of cut states: dPS helicities of ǫ l 1 ,l 2 ,l 3 which can be performed using where q µ i is an arbitrary reference momentum. From (3.3), one obtains a cut integrand as a function of {ǫ i , p i , l a }.
Before performing integral reduction using IBP, it is convenient to further isolate external polarizations. This can be achieved by choosing a set of gauge invariant basis {B α } = {B α ({ǫ i , p i })} and expanding loop integrand as (see e.g. [11,32,57]): so that all the polarization information are absorbed into basis {B α }, leaving external and loop momenta as the only variables of coefficients {f α n }. For the three-gluon setting, basis {B α } has rank-4 and can be chosen as where A i and C ij are defined as To get the coefficients f α n (p i , l a ), we only need to project the left of (3.5) by dual basis {B α }: Mutually dual basis {B α } and {B α } satisfy: See also [13] for further discussion of basis for the cases with external fermions. It is then straightforward to perform IBP reduction for the scalar function {f α n } under cut conditions, using public packages, (see e.g. [58][59][60]). After doing so we obtain the coefficients for five masters that are probed by cut channel (b), such as shown in Fig. 4. Finally, by considering other cuts, all the master coefficients can be obtained and the full form factors are given as (3.1).
The master integrals listed in Fig. 2 are known in terms of 2d harmonic polylogarithms [61,62]. Substituting the expressions, bare form factors are obtained in explicit functional form.

Anomalous dimensions of high dimensional operators
The form factors obtained in the last section contain the information of anomalous dimensions of high dimensional operators and also provide the Higgs-gluons amplitudes in the EFT. In this section, we will focus on the anomalous dimensions, which can be extracted from the UV divergences of bare form factors.
As a brief outline, in Section 4.1 we first give a review of IR and UV subtraction, and the effect of mixing between operators with different lengths is also analyzed. Dilatation operators for length-two and three operators up to dimension 16 and their eigenvalues up to O(α 2 s ) are given in Section 4.2. In Section 4.3 we take dimension 8 case as an example to cure the incompleteness of length truncation by considering the length-crossing mixing between length-4 and length-3 operators, and make correction to two-loop anomalous dimension derived in Section 4.2.

Subtraction of IR and UV Divergence
Bare form factors contain both IR and UV divergences. The general strategy is that by subtracting the universal IR divergences, one can extract the UV information unambiguously. The renormalization with multi-operators in the basis generally exhibit non-trivial operator mixing effect. We will first consider the simple case with an eigenstate operator that does not mix with other operators. Then we discuss the more general cases that involve the operator mixing effects.

Eigen operator case
For an eigen-operator, the renormalization constant is a (non-matrix) single number, The simplest example is the dimension-4 operator tr(F 2 ). The analysis of UV divergence structure is relatively simple in this case, see e.g. [11]. Below we give a brief review which will also help to set up the notation.
Perturbative expansion of renormalization factors of gauge coupling constant and operator O are: where we take renormalization under MS scheme and S ǫ = (4πe −γ E ) ǫ . The normalization of beta functions are taken to be The perturbative expansion of anomalous dimension γ is given as: 4) and the 1-loop and 2-loop order anomalous dimension can be read out from (4.1) and (4.2): 2ǫ , which can be totally determined by 1-loop data, and resultantly the intrinsic contribution for two-loop level is of the ǫ −1 part.
The coupling taken by tree-level E-point form factor of length L operator O is g E−L . Perturbative expansion of form factors can be written either in bare quantities or in renormalized quantities: O,B } can be obtained by comparing above two expressions order by order. Plugging in gauge coupling renormalization (4.1), one can get renormalization of form factors up to two-loop level as: Here δ = E − L accounts for the difference between number of external gluons and length of the operator.
To determine the UV divergences, one needs to subtract the IR divergences. This can be achieved thanks to the universality of the IR divergences, in the sense that they are independent of the type of operators but only depend on the data of external particles. The IR subtraction formula of renormalized E-gluon amplitudes up to 2-loop order is known [63,64] (see also [11]): where Operator mixing structure Next we consider the more general cases where different operators can mix with each other through loop corrections. In such case, a multi-operator generalization of renormalization multiplier Z O is needed. The main picture of the above discussion still applies, except that the renormalization constant should be taken as a matrix : Operator mixing can (and only can) take place among operators with the same canonical dimension. Also, operators with different lengths but same dimension can mix with each other. To put the discussion in a general context, we will use O L i to denotes a length-L operator labeled by i. We denote the mixing of length-L operator into length-L ′ operator at ℓ-loop level as Z (ℓ) L→L ′ , and the perturbative expansion of operator renormalization can be given as: 3 (4.14) The upper bound of k for L → (L − k) case is based on the fact that the length of terminal operator L − k should be no shorter than 2, and the loop order of a L → (L − k) process is no less than k + 1. Let us give a few general remarks regarding (4.14).
1. For length L = 2, there is only one independent operator at dimension 2∆, i.e. the descendant ∂ ∆−2 trF 2 . Since the length-2 operator is an eigenstate of anomalous dimension matrix, the length-2 operator will not mix into higher length operators along RG flow, therefore,  For form factors with operator mixing, the previous renormalization formulae (4.7)-(4.8) are generalized to In the case of minimal form factors of length-3 operators (with E = L i = 3), they reduce to Here we label the only length-2 operator as O 0 .
In the operator mixing case, we generalize the anomalous dimension (4.4) by introducing the dilation operator as 20) and the eigenvalues of the dilatation operator are anomalous dimensions. The expansion of dilatation operator may contain terms with odd power of coupling g ∼ α , D (2) can be read out from (4.1) and (4.14): L→(L−1) , (4.22) The requirement that D (ℓ) is regular as ǫ → 0 predicts several analytical structures of renormalization matrices: 1. Two-loop order length-changed mixing Z L→L is totally determined by the 1-loop data: These provide important consistency checks of the calculation.

Anomalous dimension matrices and eigenvalues
Now we apply the above strategy to compute the results of dilation operators of length-3 operators up to canonical dimension 16 and up to 2 loop order. Their eigenvalues, i.e. anomalous dimensions, are also computed up to O(α 2 s ). The renormalization of the dimension-four operator is well known, see e.g. [11], and one-loop renormalization for dimension 6 and 8 operators were considered in [65][66][67][68]. The two-loop renormalization for dimension-6 operators (also with quark operators) were obtained recently in [12,13]. Other results for operators with dimension 8-16 are given for the first time to our knowledge. These results have passed various non-trivial consistency checks, which are listed in the end of this subsection.
The basis operators we choose are given in Appendix C, labeled as O ∆,α/β,f /d,i which means it has dimension ∆, color factor f abc /d abc , numbering i and belongs to helicity sector α/β as introduced in (2.35). As we consider operators only up to length-3, the dilatation operator (4.23) in this section is defined as a truncated version as: 3→2 , (4.25) Note that Z 2→3 in D ( 3 2 ) vanishes as mentioned in (4.15). Correction with high length operators will be discussed in Section 4.3 and Appendix D.

Dimension 4
At dimension 4 there is only one independent operator O 4 = Tr(F 2 ), which has been thoroughly studied. Here we briefly review the result and also help to clarify the notations. The renormalization constants up to two loops are: 4 The double pole term of Z (2) is determined by the one-loop result as expected. According to (4.5) the anomalous dimension at 1-loop and 2-loop level can be read out:

Dimension 6
At dimension 6 there is one length-2 descendent operator ∂ 2 tr(F 2 ) and one length-3 operator tr(F 3 ) belonging to β-helicity sector according to (2.35): Two-loop minimal form factor of tr(F 3 ) was calculated in [12]. The renormalization matrix at one and two-loop level are: (4.29) At one-loop level there is no mixing, and as expected Z 3→2 so they are associated with α 3/2 s and contribute to D (3/2) , while the diagonal elements are associated with α 2 s and contribute to D (2) . The dilation matrix is straightforward to obtain using (4.25), and it reads: where for the convenience of notation, we introduce newly normalized 't Hooft couplingλ and gauge couplingĝ:λ By diagonalizing the matrix, one obtains the anomalous dimension (eigenvalues) as: (4.32)

Dimension 8
There are two length-3 basis operators at dimension 8, which are given in Table 6   We arrange the operators into a vector {O 8;0 , O 8;α;f ;1 , O 8;β;f ;1 }. According to this order, the one-loop renormalization matrix is: (4.40) At two-loop level, the Z (2) matrix is: (4.41) Using (4.25), the dilation operator is given as Note that the off-diagonal elements of the first column belong to Z

Dimension 10
There are five length-3 basis operators at dimension 10, as shown in Table 13. Together with O 10;0 = 1 4 ∂ 6 O 4 , they can be classified into three sectors: We add the lines in the matrix to separate length-2 operator and also different helicity sectors.
We can see that the operator mixing starts to appear even at one-loop order, but only within the same helicity sector.
Mixing between different helicity sectors does not happen at one-loop level, which can be understood from unitarity cut. A double cut divides one-loop minimal form factor F Still, the off-diagonal elements of the first column of (4.46) belong to Z 3→2 , so they are associated withλ 2 /ĝ, while rest of the elements are associated withλ 2 . And we can see that operators of different helicity sectors can indeed mix with each other at two loops.
The dilation operator matrix can be obtained from (4.25), and for f abc sector it is

Dimension 12
There are 10 length-3 basis operators at dimension 12, as shown in Table 14  The intrinsic two-loop renormalization matrices of f abc and d abc sector are (4.52) The dilation operator matrix is straightforward to obtain using (4.25). The anomalous dimensions are given by the eigenvalues, which read

Dimension 14
There are 15 length-3 basis operators at dimension 14, as shown in Table 15   (4.57) The intrinsic two-loop renormalization matrices of f abc and d abc sector are:  The dilation operator matrix is straightforward to obtain using (4.25) and we will not provide here. The anomalous dimensions are given by the eigenvalues, which read At operator dimension 14, the two-loop anomalous dimensions in d abc sector start to involve irrational numbers.

Dimension 16
There are 22 length-3 basis operators at dimension 16, as shown in Table 16 We (4.65) The intrinsic two-loop renormalization matrices of f abc and d abc sector are where two block matrices M , N are The dilation operator matrix is straightforward to obtain using (4.25) and we will not provide here. The anomalous dimensions are given by the eigenvalues. We summarize the anomalous dimensions in Table 8. Irrational numbers also appear in the two-loop anomalous dimensions in f abc sector at dimension 16.

Checks and analysis
Some consistency checks for our calculation have been mentioned above, and here we make a summary: 3. At a given dimension, mixing from descendent operators to non-descendent operators never takes place, such as length-2 to higher length operators in (4.15). Our results satisfy all these requirements. Some further consistency checks will be also mentioned for the computation of finite remainder function in next section. Let us make a few comments on the anomalous dimensions and dilatation matrix.

As explained in
• In Table 8, the irrational number appears in the dimension 14 and 16 cases. As eigenvalues of dilatation operators, anomalous dimensions can be obtained straightforwardly by solving characteristic equation. Alternatively, one can get their series expansions inλ up to arbitrary finite order through perturbation method introduced in quantum mechanics, which is equivalent to treat dilatation operator as a Hamiltonian of a finite system, see e.g. [69]. From perturbative calculation, one can find that whether irrational numbers appear in perturbative expansions is determined by characteristic equation of Table 9. Largest prime number in the denominator of dilatation matrix elements. the first non-degenerate order, which is the one-loop order here. As a result, if dilatation operator D (1) is of lower-triangular, then eigenvalues and eigenvectors at O(λ) must be rational, which guarantees the series expansions up to arbitrary order to be rational. On the contrary, if Z (1) contains sufficiently many non-vanishing upper-right elements, then the O(λ) characteristic equation might be complicated and irrational solution may appear. We find more and more non-vanishing upper-right elements of Z (1) emerge as operator dimension increases, so irrational numbers are expected to appear when operator dimension is high enough, which is just the case we find.
• There are large integer numbers appearing in the denominator of the dilatation matrix elements. They can be decomposed in terms of prime integer factors. We find that these prime numbers are relatively small: at one-loop no larger than (∆/2) and at two-loop no larger than (∆/2) 3 , as shown in Table 9. The pattern of denominators indicates that the elements are related to Harmonic numbers: S It would be interesting to see if the matrix elements of dilatation operator could be organizable in certain generic analytic form involving harmonic numbers.
• It is important to keep in mind that the operators we consider contain only length-2 and length-3 cases. A full basis should contain length-4 and higher length operators. By including these operators, the eigenvalues will be changed. We will investigate the effect by adding higher length operators in next subsection. To study the analytic structure of anomalous dimensions, it would be meaningful only after including the full basis operators, which we hope to discuss in the future work.

Correction from higher length operators
So far we consider only operators up to length-3. As discussed in Section 2, to simplify operator classification, we introduce an extra equivalence relation apart from E.o.M and Bianchi identity, by identifying two length-3 operators if their difference turns out to be higher length operators: At two-loop order, the length-3 operators can also mix with length-4 operators, which can modify the dilatation operator and also the anomalous dimensions computed in previous subsection. In this subsection we consider the basis operators of dimension 8 by adding length-4 ones. This allows us to compare the correct anomalous dimensions with the former result. We will find that the correction from length-change mixing between length-3 and length-4 is small at order O(10 −1 ). Let us first consider the basis of dimension-8 operators. For simplicity, we will also consider the large N c limit, thus we can ignore the double-trace length-4 operators. Besides the previously considered three lower length operators  Note that we have included a coupling g in the definition of length-4 operators. As we mentioned in Footnote 3, this will change the formula for renormalization constant and the corresponding form factors. We leave the detailed discussion to Appendix D. From the form factor results, one can extract the dilatation operators up to two loops as where the order of operators corresponds to (4.71). The last four columns contain undetermined z ij terms contributed from uncalculated Z (2) 4→4 elements. Compute the eigenvalues up to O(λ 2 ), one can find three eigenvalues (out of seven) are independent of z ij which are: Consider another length-3 operator (belonging to α-sector) defined as  (4.43), the first two eigenvalues stay unchanged as expected, since they have already appeared at lower dimensions and should not be affected by operators intrinsically emerging at dimension 8. Furthermore, the third eigenvalue is corrected only at O(λ 2 ) order from 269 18 to 235 18 with amount of O(10 −1 ). We expect this to be a general feature, as the mixing between different length operators are suppressed physically. As for our consideration, the results in Section 4.2 by truncating higher length operators are expected to provide a good approximation for the anomalous dimension of length-3 operators.

Finite remainder functions
As mentioned before, the form factors can be understood as the Higgs to three-gluon amplitudes, in which the high dimension operators correspond to the interaction vertices in the Higgs EFT. In this section, we compute the finite remainder functions of form factors, which contain the essential information of the Higgs amplitudes. O,fin is defined as in (4.9) -(4.10). We introduce the finite remainder function R One can further decompose the two-loop remainder according to their trancendentality degree as:

Transcendentality structure of remainder
In this subsection, we discuss the two-loop remainders according their transcendentality degrees. Explicit results of two-loop finite remainders are given in the ancillary file submitted together with this paper. As an example, the result of O 8;α;f ;1 is explicitly given in Appendix E.

Universal building blocks
For two-loop remainders under matched helicities, we find the transcendentality degree-4 part of two-loop minimal form factors (under match helicity) always share a universal expression: which is expected and also appears in previous computations of lower dimension operators [12,13,31,49,51,52]. 5 This implies the two-loop minimal form factor of a length-3 operator with arbitrary dimension in pure Yang-Mills theory always obeys the maximal transcendentality principle.
For two-loop remainders under unmatched helicities, there are no degree 4 or 3 parts, in accord with the vanishing of ǫ −4 , ǫ −3 poles in bare form factors. Finite remainders and poles of the same degree originate from the same term in the master integral coefficients, so they usually vanish simultaneously. The absence of degree 4 and 3 poles at two-loop level can be traced back to the absence of one-loop divergence. As mentioned in section 4.2, under unmatched helicity the tree-level form factor is zero and the one-loop form factor only has rational term, so divergence subtraction formula from (4.8) and (4.10) becomes which explicitly shows the leading singularity is of O(ǫ −2 ) from I (1) (ǫ)F B , and no term can contribute to ǫ −3 , ǫ −4 .
Apart from maximal transcendental universality, degree-3 and degree-2 parts also signify some universal structure, in the sense that complicated transcendental functions can always be absorbed into a set of universal building blocks, and no other polylogarithm functions like Li 2 , Li 3 are left outside these basis functions. 6 Building blocks of degree-3 part are six functions {T 3 [σ(x), σ(y), σ(z)]|σ ∈ S 3 } together with π 2 log and ζ 3 , where T 3 (u, v, w) is given as Similar function has appeared in the N = 4 form factors [45,47,49]. Building blocks of degree-2 part are three functions {T 2 [σ(x), σ(y)]|σ ∈ Z 3 } together with log 2 and π 2 , where T 2 (u, v) is given as (see also [13]) When expanding the form factor remainders in these building blocks, the coefficients in front of them are just rational functions of u, v, w, see examples in Appendix E. 5 The computation here in QCD using Catani IR subtraction scheme, and the expression is slightly different (as purely a scheme change) from the N = 4 results which are based on the BDS subtraction scheme [70]. 6 When quark is added, T3, T2 no longer compose complete basis for polylogarithms, see [13].
Comment on log(q 2 ) terms As mentioned in (5.4), remainder functions also contain two type of terms that are linear and quadratic in log(q 2 ) respectively. All these terms have simple universal structures which we explain below. Coefficient of log 2 (−q 2 ) term of O i equals to the ǫ −2 residue of j (Z (2) ) j i F O j . This is because every master integral contains an overall factor (−q 2 ) −2ǫ , whose ǫ-expansion produces 1 and log 2 (−q 2 )ǫ 2 at the same time. So leading singularity 1/ǫ 2 can be recovered back to 1/ǫ 2 (1 + log 2 (−q 2 )ǫ 2 ) and this is the only origin of log 2 (−q 2 ) term. Plugging in (4.24), we can write log 2 (−q 2 ) term as is the tree-level scalar factor of O j , as introduced in Table 10. Since matrix elements of Z (1) is nonzero only for operators within same α or β sector, the log 2 (−q 2 ) terms vanish under unmatched helicity. Terms linear in log(−q 2 ) can be determined by data of one-loop finite remainders and renormalization matrices Z (1) , Z (2) . Coefficients of linear log(−q 2 ) terms contain degree 3,2,1,0 parts. The generic expression of this part is: is the one-loop finite remainder, whose log(−q 2 ) coefficient is j (ǫZ (1) ) j i f (0),± O j . The degree 3 and 2 parts of log(−q 2 ) coefficient are universal. Both (5.10) and (5.11) provide nice consistency checks for our calculation.

Cancellation of spurious poles
While the transcendentality functions in the remainder take very simple universal structures as discussed in last subsection, their rational coefficients depend on the dimension of operators and are the main complication of the remainders. In particular, they contain spurious poles which cancel only after summarizing the results.
Spurious poles exist at transcendentality degree 3,2,1,0. They contain high order poles of u n , v n , w n with n > 1, as well as linear combination poles u + v and u + w. In Table 11 we summarize the leading spurious poles in the form factor remainders.
A non-trivial feature is that the cancellation of these spurious poles takes place across different transcendentality degrees. As a concrete non-trivial example, we consider the remainder function of O 8;α;f ;1 under match helicity (−, −, +). Explicit expressions of two-loop remainder of O 8;α;f ;1 (with degree 3,2,1,0 parts) can be found in Appendix E. We first summarize the property of pole structures:  In the following we discuss these spurious poles separately and show that they explicitly cancel in the full remainder.

1/u m -pole
To analyze the 1 u m -poles, one needs to consider the limit u → 0. Polylogarithm functions can be expanded in this limit, for example: To cancel the O(u −5 ) term, one needs to consider the contribution from the degree-1 part. Concretely, one can extract the residue for 1/u 5 pole terms from various degree parts as: One can check that after expanding the polylogarithm functions to appropriate orders in u, all u −k poles vanish except for the physical 1/u pole. Explicitly, residues of sub-leading poles contained by different transcendentality degree parts are: All these 1/u m poles do not cancel within single transcendentality degree, but only after the sum of different degree parts.

1/v m -pole
For 1/v m -poles the analysis is similar, one can take limit v → 0 and expand polylogarithm functions in v. After doing so, 1 v 5 pole in degree-2 part vanishes. Residues of sub-leading poles contained by four parts become: These 1 v m poles cancel after summing over degree 2,1,0 parts. Poles of 1 w m are related by symmetry and are cancelled in the same way.

1/(u + v) m -pole
Similarly, take limit v → −u and expand polylogarithm functions in (u + v). After doing so, leading poles 1 (u+v) 2 in degree 1 and 0 parts cancel within their own degree, and residues of sub-leading 1/(u + v) pole cancel when summing over degree 1 and 0 parts: , deg-1 : 5 18u . (5.19) 6 Summary and outlook In this paper we study the two-loop renormalization of high dimensional QCD local operators and related Higgs amplitudes in Higgs EFT. We have discussed the classification of operators and the construction of their basis, for which we apply both field theory method and on-shell spinor helicity formalism. Given the operator basis, we explain the unitarity-IBP method to compute their form factors. Explicitly the two-loop minimal form factors of length-3 operators up to dimension 16 are computed. The on-shell methods have played important roles for both operator construction and loop computation. Based on the form factor results, we extract the renormalization constant matrix up to two loop order. At one loop, the operator mixing starts at dimension 10, while at two loops it begins from dimension 6. These show that the operator mixing is a general feature for QCD high dimension operators. By diagonalizing the dilatation matrix, we compute the anomalous dimensions. To study the correction by high length operators, we consider the dimension 8 case and include the length-4 operators to form a complete set of basis. We find that the anomalous dimension only receives a small correction. This implies that the mixing with high length operators are small, and considering only length-3 operators is expected to be a good approximation.
The form factor results also provide the Higgs to three-gluon amplitudes in Higgs EFT, where the high dimension operators correspond to high order of top mass corrections. The non-trivial information of the amplitudes is captured by the finite remainder function. We find that for all length-3 operators up to dimension 16, the finite remainders always share the same maximal transcendentality part. This implies that the maximal transcendentality principle is independent of the dimension of operators. Patterns for the lower transcendental parts are also analyzed.
There are several straightforward generalizations based on this paper. In this paper, we focus on the pure gluon sector of QCD. It is also important to consider the operators with quark fields, similar to the dimension-6 operators considered to two-loop order in [13]. It would be also interesting to consider general parity odd operators, such examples contain the Weinberg's operator [71], see e.g. [72] for a recent study. There have been many studies of constructing operator basis in standard model EFT (see e.g. [73][74][75][76][77][78][79][80][81]) and consider their renormalization using both conventional and on-shell formalism (see e.g. [82][83][84][85][86][87][88][89][90]). While the state-of-the-art computation is mostly at one-loop order, it would be worth extending the method developed in this paper to to general two-loop renormalization in SMEFT (see also [91]). values 0, 2, 4, 6. We will show there are only two independent primitive classes at length-3: First, forn d = 6, the three field strengths must take forms like F ij , F kl , F mn . There is at least one block contains nonzero covariant derivatives, so one can rewrite this block using Bianchi identity. In this way one Lorentz index will be moved from a D to the F , which meansn d is reduced from 6 to 4. For example: Second, forn d = 4, the three field strengths must take forms like F ij , F jk , F mn . Among them F mn is the special one because it shares no Lorentz index with the other two. We consider following two cases: 1 If there is at least one covariant derivative, say D i or D k , acting on F mn , we can make use of Bianchi identity to move the Lorentz index from this D to the F so thatn d is reduced from 4 to 2. For example: 2 If there is no covariant derivative acting on F mn , then other two blocks must both contain nonzero Ds. We can rewrite one of those two blocks using Bianchi identity, For example: The first term on the r.h.s is already reduced ton d = 2. The second term still haŝ n d = 4, but now F 23 becomes the special F mn , and there are nonzero Ds acting on it, so this actually belongs to the previous case.
In summary, every primitive operator withn d = 4 can be reduced ton d = 2. Third, forn d = 2, the three field strengths may take forms like (F ij F kl , F kl ) or (F ij , F jk , F kl ). We consider various different cases as follows: a For configuration (F ij , F kl , F kl ), the only permitted operator is D 3 F 12 D 4 F 12 F 34 . Other seemingly qualified operators like D 34 F 12 F 12 F 34 and F 12 D 34 F 12 F 34 actually have higher length and therefore should not be taken into account, e.g., b For configuration (F ij , F jk , F kl ), we can see F jk is the special one because it shares two Lorentz indices with the other two F s.
b.1) If there is only one covariant derivative, say D i or D l , acting on F jk , we can rewrite this block using Bianchi identity, relabel the Einstein indices and finally obtain an operator belonging to case a. For example: b.2) If there are two Ds acting on F jk , the operator turns out to be a higher length one, so it should be dropped out: Forn d = 2 operators, we prefer configuration (F ij , F kl , F kl ) instead of (F ij , F jk , F kl ), because the former one has a simpler expression under decomposition F ααββ = ǫ αβfαβ + ǫαβf αβ . The components from two F kl 's must be both self-dual or both anti-self-dual because contraction between f orf and antisymmetric tensor gives zero. If we probe the minimal form factor under configuration (−, −, +) or (+, +, −), the particle with the opposite helicity must be emitted from F ij . That is why we choose the only independentn d = 2 primitive operator from case a.
Above analysis shows there are only two independent primitive class as shown in (A.1). Based on this, one can construct basis of length-3 operators for any given dimension ∆ 0 systematically by inserting pairs of identical D i s into the primitive ones until ∆ 0 is reached.

B Preliminary operator basis O ′′
In this appendix, we provide length-3 basis operators constructed by inserting DD pairs to primary operators, as described in Section 2.2. The underlined number (n + m) represents that there are n operators created from primitive operator O P1 = Tr(D 1 F 23 D 4 F 23 F 14 ) and m ones from the primitive operator O P2 = Tr(F 12 F 13 F 23 ).  As discussed in Section 2, the above operators are not the final basis choice, we still need to solve descendant relations as constraints and symmetrize the operators which do not have symmetric properties. As shown in Appendix C, the final basis operators we choose are linear combinations of these preliminary ones.

C Operator basis up to dimension 16
In Section 2.2 and 2.3 we've explained how to find length-3 basis operators from both field theory method and on-shell method by showing an example at dimension 10, and these two strategies apply for general operator dimensions.
As mentioned in (2.24) there is a freedom in choosing which non-descendant operators to be replaced by the descendant ones. We require that the scalar factors f (0),± O (see Table 10, here + for α-sector, − for β-sector) of the chosen operators should have the following forms: Here u = s 12 s 123 , v = s 23 s 123 , w = s 13 s 123 , and ∆ 0 is the dimension of the operator. Under this constraint we can reduce the number of candidate operators and do not violate the completeness. However, there is still some choice freedom. For example, the number of independent non-descendent operators for dimension-12 (d abc , −−+)-sector is 1, but there are two non-descendent scalar factors satisfying the stated forms: u(w − v)s 2 123 and w 2 − v 2 s 2 123 . Here during computation we choose the former. Generally speaking, the choice freedom can not be avoided, and the different operator choices correspond to different tree-level scalar factor.
The chosen basis operators are listed in following tables, labeled as O ∆ 0 ,α/β,f /d,i , where ∆ 0 is the dimension operator, α/β denotes helicity sector (introduced in (2.35) ), f /d denotes color factor f abc /d abc , and i denotes numbering. All the new basis operators are linear combinations of old ones given in (B.1)-(B.6) labeled as O ′′ ∆ 0 ,i . In the following tables, we show the tree level minimal form factor of every basis operator in spinor-helicity formalism. As for length-3 operators, spinor structure of minimal form factor is universal within each given helicity sector. Concretely with tree-level minimal form factor A 2 . Basis operators at dimension 8, 10, 12, 14, 16 are given in Table 12, 13, 14, 15, 16.  Since descendants are always kept into basis choice, the remainder list of basis operators at a certain dimension cover all the remainders of lower dimensional basis.

D Dilatation operator under new operator definition
In this appendix, we provide some details regarding the length-4 operators in Section 4.3.

Dilatation operator with length-4
Back to perturbative expansion of dilatation operator, by adding length-4 operators using (4.22)-(4.23), we modify the length-3 formula (4.25) as 4→3 , 4→2 . (D.1) Above formulae correspond to operator definition which does not include gauge coupling as mentioned in Footnote 3. Again notice there is no mixing from length-2 operator to higher length ones as mentioned in (4.15). For dimension-8 case, there is no operator with length higher than 4, so including length-4 operators makes the operator basis complete and (D.1) provide the full result. Renormalization matrix is obtained from UV subtraction of loop form factors. Using 4→4 Z 3→3 , Z 3→2 . We also need to consider non-minimal form factors: 4→4 one needs to compute two-loop 4-gluon form factors of length-4 operators, while to determine Z (2) 4→4 one needs to take three-loop calculation, and these two types of data will not be discussed in this paper. Nevertheless, to determine the two-loop anomalous dimension of length-3 operators, the information of Z Operators with coupling in the definition As we mentioned in Footnote 3, if we multiply all the basis length-4 operators by an overall coupling g, the perturbative order of length-changed matrix elements will change and running effect of gauge coupling will enter in, and as a result perturbative expansion (D.1) will change to: 3→2 , (D.6) The basis (4.71) we choose in Section 4.3 correspond to this operator redefinition. Notice perturbative expansion (D.2)-(D.5) will not change because the order change of Z-matrix elements will be compensated for the order change of form factors.