Pseudo-scalar Form Factors at Three Loops in QCD

The coupling of a pseudo-scalar Higgs boson to gluons is mediated through a heavy quark loop. In the limit of large quark mass, it is described by an effective Lagrangian that only admits light degrees of freedom. In this effective theory, we compute the three-loop massless QCD corrections to the form factor that describes the coupling of a pseudo-scalar Higgs boson to gluons. Due to the axial anomaly, the pseudo-scalar operator for the gluonic field strength mixes with the divergence of the axial vector current. Working in dimensional regularization and using the 't~Hooft-Veltman prescription for the axial vector current, we compute the three-loop pseudo-scalar form factors for massless quarks and gluons. Using the universal infrared factorization properties, we independently derive the three-loop operator mixing and finite operator renormalisation from the renormalisation group equation for the form factors, thereby confirming recent results in the operator product expansion. The finite part of the three-loop form factor is an important ingredient to the precise prediction of the pseudo-scalar Higgs boson production cross section at hadron colliders. We discuss potential applications and derive the hard matching coefficient in soft-collinear effective theory.


Introduction
Form factors are the matrix elements of local composite operators between physical states. In the calculation of scattering cross sections, they provide the purely virtual corrections. For example, in the context of hard scattering processes such as Drell-Yan [1,2] and the Higgs boson production in gluon fusion [3][4][5][6][7][8][9][10][11][12][13][14][15], the form factors corresponding to the vector current operator ψγ µ ψ and the gluonic operator G a µν G a,µν contribute, respectively. Here ψ is the fermionic field operator and G a µν is the field tensor of the non-Abelian gauge field A a µ corresponding to the gauge group SU(N). In Quantum Chromodynamics (QCD) the form factors can be computed order by order in the strong coupling constant using perturbation theory. Beyond leading order, the ultraviolet (UV) renormalisation of the form factors involves the renormalisation of the composite operator itself, besides the standard procedure for coupling constant and external fields.
The resulting UV finite form factors still contain divergences of infrared (IR) origin, namely, soft and collinear divergences due to the presence of massless gluons and quarks/antiquarks in the theory. The inclusive hard scattering cross sections require, in addition to the form factor, the real-emission partonic subprocesses as well as suitable mass factorisation kernels for incoming partons. The soft divergences in the form factor resulting from the gluons cancel against those present in the real emission processes and the mass factorisation kernels remove the remaining collinear divergences rendering the hadronic inclusive cross section IR finite. While these IR divergences cancel among various parts in the perturbative computations, they can give rise to logarithms involving physical scales and kinematic scaling variables of the processes under study. In kinematical regions where these logarithms become large, they may affect the convergence and reliability of the perturbation series expansion in powers of the coupling constant. The solution for this problem goes back to the pioneering work by Sudakov [16] on the asymptotic behaviour of the form factor in Quantum Electrodynamics: all leading logarithms can be summed up to all orders in perturbation theory. Later on, this resummation was extended to non-leading logarithms [17] and systematised for non-Abelian gauge theories [18]. Ever since, form factors have been central to understand the underlying structure of amplitudes in gauge theories.
The infrared origin of universal logarithmic corrections to form factors [19] and scattering amplitudes results in a close interplay between resummation and infrared pole structure. Working in dimensional regularisation in d = 4 + dimensions, these poles appear as inverse powers in the Laurent expansion in . In a seminal paper, Catani [20] proposed a universal formula for the IR pole structure of massless two-loop QCD amplitudes of arbitrary multiplicity (valid through to double pole terms). This formula was later on justified systematically from infrared factorization [21], thereby also revealing the structure of the single poles in terms of the anomalous dimensions for the soft radiation. In [13], it was shown that the single pole term in quark and gluon form factors up to two loop level can be shown to decompose into UV (γ I , I = q, g) and universal collinear (B I ), color singlet soft (f I ) anomalous dimensions, later on observed to hold even at three loop level in [22]. An all loop conjecture for the pole structure of the on-shell multi-loop multi-leg amplitudes in SU(N) gauge theory with n f light flavors in terms of cusp (A I ), collinear (B I ) and soft anomalous dimensions (Γ IJ , f I -colour non-singlet as well as singlet) was proposed by Becher and Neubert [23] and Gardi and Magnea [24], generalising the earlier results [20,21]. The validity of this conjecture beyond three loops depends on the presence/absence of non-trivial colour correlations and crossing ratios involving kinematical invariants [25] and there exists no all-order proof at present. The three-loop expressions for cusp, collinear and colour singlet soft anomalous dimensions were extracted [26,27] from the three loop flavour singlet [28] and non-singlet [29] splitting functions, thereby also predicting [22] the full pole structure of the three-loop form factors.
The three-loop quark and gluon form factors through to finite terms were computed in [30][31][32][33] and subsequently extended to higher powers in the expansion [34]. These results were enabled by modern techniques for multi-loop calculations in quantum field theory, in particular integral reduction methods. These are based on integration-by-parts (IBP) [35,36] and Lorentz invariance (LI) [37] identities which reduce the set of thousands of multi-loop integrals to the one with few integrals, so called master integrals (MIs) in dimensional regularisation. To solve these large systems of IBP and LI identities, the Laporta algorithm [38], which is based on lexicographic ordering of the integrals, is the main tool of choice. It has been implemented in several computer algebra codes [39][40][41][42][43][44]. The MIs relevant to the form factors are single-scale three-loop vertex functions, for which analytical expressions were derived in Refs. [31,[45][46][47][48][49].
Recently, some of us have applied these state-of-the-art methods to accomplish the task of computing spin-2 quark and gluon form factors up to three loops [50] level in SU(N) gauge theory with n f light flavours. These form factors are ingredients to the precise description of production cross sections for graviton production, that are predicted in extensions of the Standard Model. In addition, the spin-2 form factors relate to operators with higher tensorial structure and thus provide the opportunity to test the versatility and robustness of calculational techniques for the vertex functions at three loop level. The results [50] also confirm the universality of the UV and IR structure of the gauge theory amplitudes in dimensional regularisation.
In the present work, we derive the three-loop corrections to the quark and gluon form factors for pseudo-scalar operators. These operators appear frequently in effective field theory descriptions of extensions of the Standard Model. Most notably, a pseudo-scalar state coupling to massive fermions is an inherent prediction of any two-Higgs doublet model [51][52][53][54][55][56][57][58]. In the limit of infinite fermion mass, this gives rise to the operator insertions considered here. The recent discovery of a Standard-Model-like Higgs boson at the LHC [59,60] has not only revived the interest in such Higgs bosons but also prompted the study of the properties of the discovered boson to identify either with lightest scalar or pseudoscalar Higgs bosons of extended models. Such a study requires precise predictions for their production cross sections. In the context of a CP-even scalar Higgs boson, results for the inclusive production cross section in the gluon fusion are available up to N 3 LO QCD [10][11][12]15], based on an effective scalar coupling that results from integration of massive quark loops that mediate the coupling of the Higgs boson to gluons [61][62][63]. On the other hand for the CP-odd pseudo scalar, only NNLO QCD results [12,[64][65][66][67] in the effective theory [68] are known. The exact quark mass dependence for scalar and pseudo-scalar production is known to NLO QCD [6,69], and is usually included through a re-weighting of the effective theory results. Soft gluon resummation of the gluon fusion cross section has been performed to N 3 LL for the scalar case [70][71][72][73][74][75][76][77][78] and to NNLL for the pseudo-scalar case [79]. A generic threshold resummation formula valid to N 3 LL accuracy for colour-neutral final states was derived in [78], requiring only the virtual three-loop amplitudes as process-dependent input. The numerical impact of soft gluon resummation in scalar and pseudo-scalar Higgs boson production and its combination with mass corrections is reviewed comprehensively in [80]. The three-loop corrections to the pseudo-scalar form factors computed in this article are an important ingredient to the N 3 LO and N 3 LL gluon fusion cross sections [81] for pseudoscalar Higgs boson production, thereby enabling predictions at the same level of precision that is attained in the scalar case.
The framework of the calculation is outlined in Section 2, where we describe the effective theory [68]. Due to the pseudo-scalar coupling, one is left with two effective operators with same quantum number and mass dimensions, which mix under renormalisation. Since these operators contain the Levi-Civita tensor as well as γ 5 , the computation of the matrix elements requires additional care in 4+ dimensions where neither Levi-Civita tensor nor γ 5 can be defined unambiguously. We use the prescription by 't Hooft and Veltman [82,83] to define γ 5 . We describe the calculation in Section 3, putting particular emphasis on the UV renormalisation. Exploiting the universal IR pole structure of the form factors, we determine the UV renormalisation constants and mixing of the effective operators up to three loop level. We also show that the finite renormalisation constant, known up to three loops [83], required to preserve one loop nature of the chiral anomaly, is consistent with anomalous dimensions of the overall renormalisation constants. As a first application of our form factors, we compute the hard matching functions for N 3 LL resummation in soft-collinear effective theory (SCET) in Section 4. Section 5 summarises our results and contains an outlook on future applications to precision phenomenology of pseudo-scalar Higgs production.

The Effective Lagrangian
A pseudo-scalar Higgs boson couples to gluons only indirectly through a virtual heavy quark loop. This loop can be integrated out in the limit of infinite quark mass. The resulting effective Lagrangian [68] encapsulates the interaction between a pseudo-scalar Φ A and QCD particles and reads: where the operators are defined as The Wilson coefficients C G and C J are obtained by integrating out the heavy quark loop, and C G does not receive any QCD corrections beyond one loop due to the Adler-Bardeen theorem [84], while C J starts only at second order in the strong coupling constant. Expanded in a s ≡ g 2 s /(16π 2 ) = α s /(4π), they read In the above expressions, G µν a and ψ represent gluonic field strength tensor and light quark fields, respectively and G F is the Fermi constant and cotβ is the mixing angle in a generic Two-Higgs-Doublet model. a s ≡ a s µ 2 R is the strong coupling constant renormalised at the scale µ R which is related to the unrenormalised one,â s ≡ĝ 2 s /(16π 2 ) througĥ and µ is the scale introduced to keep the strong coupling constant dimensionless in d = 4 + space-time dimensions. The renormalisation constant Z as [85] is given by up to O(a 3 s ). β i are the coefficients of the QCD β functions which are given by [85] with the SU(N) QCD color factors n f is the number of active light quark flavors.

Treatment of γ 5 in Dimensional Regularization
Higher order calculations of chiral quantities in dimensional regularization face the problem of defining a generalization of the inherently four-dimensional objects γ 5 and ε µνρσ to values of d = 4. In this article, we have followed the most practical and self-consistent definition of γ 5 for multiloop calculations in dimensional regularization which was introduced by 't Hooft and Veltman through [82] Here, ε µνρσ is the Levi-Civita tensor which is contracted as and all the Lorentz indices are considered to be d-dimensional [83]. In this scheme, a finite renormalisation of the axial vector current is required in order to fulfill chiral Ward identities and the Adler-Bardeen theorem. We discuss this in detail in Section 3.2 below.

Pseudo-scalar Quark and Gluon Form Factors
The quark and gluon form factors describe the QCD loop corrections to the transition matrix element from a color-neutral operator O to an on-shell quark-antiquark pair or to two gluons. For the pseudo-scalar interaction, we need to consider the two operators where, n = 0, 1, 2, 3, . . . . In the above expressions |M λ,(n) β is the O(â n s ) contribution to the unrenormalised matrix element for the transition from the bare operator [O λ ] B (λ = G, J) to a quark-antiquark pair (β = q) or to two gluons (β = g). The expansion of these quantities in powers ofâ s is performed through where, Q 2 = −2 p 1 .p 2 and p i s (p 2 i = 0) are the momenta of the external quarks and gluons. Note that |M G,(n) q and |M J,(n) g start from n = 1 i.e. from one loop level.

Calculation of the Unrenormalised Form Factors
The calculation of the unrenormalised pseudo-scalar form factors up to three loops follows closely the steps used in the derivation of the three-loop scalar and vector form factors [31,33]. The Feynman diagrams for all transition matrix elements (Eq. (3.1), Eq. (3.2)) are generated using QGRAF [86]. The numbers of diagrams contributing to three loop amplitudes are 1586 for |M where all the external particles are considered to be on-shell. The raw output of QGRAF is converted to a format suitable for further manipulation. A set of in-house routines written in the symbolic manipulating program FORM [87] is utilized to perform the simplification of the matrix elements involving Lorentz and color indices. Contributions arising from ghost loops are taken into account as well since we use Feynman gauge for internal gluons. For the external on-shell gluons, we ensure the summing over only transverse polarization states by employing an axial polarization sum: where p i is the i th -gluon momentum, q i is the corresponding reference momentum which is an arbitrary light like 4-vector and s stands for spin (polarization) of gluons. We choose q 1 = p 2 and q 2 = p 1 for our calculation. Finally, traces over the Dirac matrices are carried out in d dimensions.
The expressions involve thousands of three-loop scalar integrals. However, they are expressible in terms of a much smaller set of scalar integrals, called master integrals (MIs), by use of integration-by-parts (IBP) [35,36] and Lorentz invariance (LI) [37] identities. These identities follow from the Poincare invariance of the integrands, they result in a large linear system of equations for the integrals relevant to given external kinematics at a fixed loop-order. The LI identities are not linearly independent from the IBP identities [88], their inclusion does however help to accelerate the solution of the system of equations. By employing lexicographic ordering of these integrals (Laporta algorithm, [38]), a reduction to MIs is accomplished. Several implementations of the Laporta algorithm exist in the literature: AIR [39], FIRE [40], Reduze2 [41,42] and LiteRed [43,44]. In the context of the present calculation, we used LiteRed [43,44] to perform the reductions of all the integrals to MIs.
Each three-loop Feynman integral is expressed in terms of a list of propagators involving loop momenta that can be attributed to one of the following three sets (auxiliary topologies, [31]) (3.5) In the above sets To accomplish this, we have used the package Reduze2 [41,42]. In each set in Eq. (3.5), D s are linearly independent and form a complete basis in a sense that any Lorentz-invariant scalar product involving loop momenta and external momenta can be expressed uniquely in terms of D s from that set. As a result, we can express the unrenormalised form factors in terms of 22 topologically different master integrals (MIs) which can be broadly classified into three different types: genuine three-loop integrals with vertex functions (A t,i ), three-loop propagator integrals (B t,i ) and integrals which are product of one-and two-loop integrals (C t,i ). These integrals were computed analytically as Laurent series in in [45][46][47][48][49] and are collected in the appendix of [31]. Inserting those, we obtain the final expressions for the unrenormalised (bare) form factors that are listed in Appendix A.

UV Renormalisation
To obtain ultraviolet-finite expressions for the form factors, a renormalisation of the coupling constant and of the operators is required. The UV renormalisation of the operators [O G ] B and [O J ] B involves some non-trivial prescriptions. These are in part related to the formalism used for the γ 5 matrix, section 2.2 above.
This formalism fails to preserve the anti-commutativity of γ 5 with γ µ in d dimensions. In addition, the standard properties of the axial current and Ward identities, which are valid in a basic regularization scheme like the one of Pauli-Villars, are violated as well. As a consequence, one fails to restore the correct renormalised axial current, which is defined as [83,89] in dimensional regularization. To rectify this, one needs to introduce a finite renormalisation constant Z s 5 [84,90] in addition to the standard overall ultraviolet renormalisation constant Z s M S within the M S-scheme: By evaluating the appropriate Feynman diagrams explicitly, Z s M S can be computed, however the finite renormalisation constant is not fixed through this calculation. To determine Z s 5 one has to demand the conservation of the one loop character [91] of the operator relation of the axial anomaly in dimensional regularization: The bare operator [O J ] B is renormalised multiplicatively exactly in the same way as the axial current J µ 5 through whereas the other one [O G ] B mixes under the renormalisation through with the corresponding renormalisation constants Z GG and Z GJ . The above two equations can be combined to express them through the matrix equation In the above expressions Z JG = 0 to all orders in perturbation theory , We determine the above-mentioned renormalisation constants Z s M S , Z GG , Z GJ up to O a 3 s from our calculation of the bare on-shell pseudo-scalar form factors described in the previous subsection. This procedure provides a completely independent approach to their original computation, which was done in the operator product expansion [92]. Our approach to compute those Z ij is based on the infrared evolution equation for the form factor, and will be detailed in Section 3.3 below. Moreover, we can fix Z s 5 up to O(a 2 s ) by demanding the operator relation of the axial anomaly (Eq. (3.8)). Using these overall operator renormalisation constants along with strong coupling constant renormalisation through Z as , Eq. (2.5), we obtain the UV finite on-shell quark and gluon form factors.
To define the UV renormalised form factors, we introduce a quantity S λ β , constructed out of bare matrix elements, through (3.14) Expanding the quantities appearing on the right hand side of the above equation in powers of a s : we can write Then the UV renormalised form factors corresponding to O G are defined as Similarly, for defining the UV finite form factors for the other operator O J we introduce Expanding Z s M S and |M λ β in powers of a s , following Eq. (3.15), we get (3.20) With these we define the UV renormalised form factors corresponding to O J through The finite renormalisation constant Z s 5 is multiplied in Eq. (3.19) to restore the axial anomaly equation in dimensional regularisation. We determine all required renormalisation constants from consistency conditions on the universal structure of the infrared poles of the renormalised form factors in the next section, and use these constants to derive the UV-finite form factors in Section 3.4.

Infrared Singularities and Universal Pole Structure
The renormalised form factors are ultraviolet-finite, but still contain divergences of infrared origin. In the calculation of physical quantities (which fulfill certain infrared-safety criteria [93]), these infrared singularities are cancelled by contributions from real radiation processes that yield the same observable final state, and by mass factorization contributions associated with initial-state partons. The pole structures of these infrared divergences arising in QCD form factors exhibit some universal behaviour. The very first successful proposal along this direction was presented by Catani [20] (see also [21]) for one and twoloop QCD amplitudes using the universal subtraction operators. The factorization of the single pole in quark and gluon form factors in terms of soft and collinear anomalous dimensions was first revealed in [13] up to two loop level whose validity at three loop was later established in the article [22]. The proposal by Catani was generalized beyond two loops by Becher and Neubert [23] and by Gardi and Magnea [24]. Below, we outline this behaviour in the context of pseudo-scalar form factors up to three loop level, following closely the notation used in [72].
The unrenormalised form factors F λ β (â s , Q 2 , µ 2 , ) satisfy the so-called KG-differential equation [16][17][18]94] which is dictated by the factorization property, gauge and renormalisation group (RG) invariances: where all poles in the dimensional regulator are contained in the Q 2 independent function K λ β and the finite terms in → 0 are encapsulated in G λ β . RG invariance of the form factor implies β,i on the right hand side are the i-loop cusp anomalous dimensions. It is straightforward to solve for K λ β in Eq. (3.24) in powers of bare strong coupling constantâ s by performing the following expansion The solutions K λ β,i ( ) consist of simple poles in with the coefficients consisting of A λ β,i and β i . These can be found in [72,73]. On the other hand, the renormalisation group can be solved. The solution contains two parts, one is dependent on µ 2 R whereas the other part depends only the boundary point µ 2 R = Q 2 . The µ 2 R dependent part can eventually be expressed in terms of A λ β : The boundary term can be expanded in powers of a s as The solutions of K λ β and G λ β enable us to solve the KG equation (Eq. (3.23)) and thereby facilitate to obtain the ln F λ β (â s , Q 2 , µ 2 , ) in terms of A λ β,i , G λ β,i and β i which is given by [72] ln All these form factors are observed to satisfy [13,22] the following decomposition in terms of collinear (B λ β ), soft (f λ β ) and UV (γ λ β ) anomalous dimensions: where the constants C λ β,i are given by [73] C λ β, In the above expressions, X λ β,i with X = A, B, f and γ λ β,i are defined through Within this framework, we will now determine this universal structure of IR singularities of the pseudo-scalar form factors. This prescription will be used subsequently to determine the overall operator renormalisation constants. We begin with the discussion of form factors corresponding to O J . The results of the form factors F J β for β = q, g, which have been computed up to three loop level in this article are being used to extract the unknown factors, γ J β,i and g J,k β,i , by employing the KG equation. Since the F J β satisfy KG equation, we can obtain the solutions Eq. (3.28) along with Eq. (3.29) and Eq. (3.30) to examine our results against the well known decomposition of the form factors in terms of the quantities X J β . These are universal, and appear also in the vector and scalar quark and gluon form factors [13,22]. They are known [13,28,[95][96][97] up to three loop level in the literature. Using these in the above decomposition, we obtain γ J β,i . The other process dependent constants, namely, g J,k β,i can be obtained by comparing the coefficients of k in Eq. (3.29) at every order inâ s . We can get the quantities γ J g,i and g J,k g,i up to two loop level, since this process starts at one loop. From gluon form factors we get Similarly, from the quark form factors we obtain Note that γ J q,i = γ J g,i which is expected since these are the UV anomalous dimensions associated with the same operator [O J ] B . The γ J β,i are further used to obtain the overall operator renormalisation constant Z s M S through the RGE: The general solution of the RGE is obtained as By substituting the results of γ J β,i in the above solution we get Z s M S up to O(a 3 s ): which agrees completely with the known result in [83]. In order to restore the axial anomaly equation in dimensional regularization (see Section 3.2 above), we must multiply the Z s M S [O J ] B by a finite renormalisation constant Z s 5 , which reads [83] Following the computation of the operator mixing constants below, we will be able to verify explicitly that this expression yields the correct expression for the axial anomaly.  (3.28) in conjunction with Eq. (3.29) and Eq. (3.30) to compare our results against the universal decomposition of the form factors in terms of the constants X G β . Upon substituting the existing results of the quantities A G β,i , B G β,i and f G β,i up to three loops, which are obtained in case of quark and gluon form factors, we determine the anomalous dimensions γ G β,i and the constants g G,k β,i . However, it is only possible to get the factors γ G q,i and g G,k q,i up to two loops because of the absence of a tree level amplitude in the quark initiated process for the operator O G . Since [F G β ] R are UV finite, the anomalous dimensions γ G β,i must be equal to the anomalous dimension corresponding to the renormalisation constant Z GG . This fact is being used to determine the overall renormalisation constants Z GG and Z GJ up to three loop level where these quantities are parameterized in terms of the newly introduced anomalous dimensions γ ij through the matrix equation The general solution of the RGE up to a 3 s is obtained as where, γ ij is expanded in powers of a s as γ ij = ∞ n=1 a n s γ ij,n . (3.41) Demanding the vanishing of γ G β,i , we get In addition to the demand of vanishing γ G β,i , it is required to use the results of γ JJ and γ JG , which are implied by the definition, Eq. (3.39), up to O(a 2 s ) to determine the abovementioned γ GG and γ GJ up to the given order. This is a consequence of the fact that the operators mix under UV renormalisation. Following Eq. (3.39) along with Eq. (3.13), Eq. (3.37) and Eq. (3.38), we obtain As it happens, we note that γ JJ 's are -dependent and in fact, this plays a crucial role in determining the other quantities. Our results are in accordance with the existing ones, γ GG and γ GJ , which are available up to O(a 2 s ) [83] and O(a 3 s ) [92], respectively. In addition to the existing ones, here we compute the new result of γ GG at O(a 3 s ). It was observed through explicit computation in the article [83] that holds true up to two loop level but there was no statement on the validity of this relation beyond that order. In [92], it was demonstrated in the operator product expansion that the relation holds even at three loop. Here, through explicit calculation, we arrive at the same conclusion that the relation is still valid at three loop level which can be seen if we look at the γ GG,3 in Eq. (3.42) which is equal to the β 2 . Before ending the discussion of γ ij , we examine our results against the axial anomaly relation. The renormalisation group invariance of the anomaly equation (Eq. (3.8)), see [83], gives Through our calculation up to three loop level we find that our results are in complete agreement with the above anomaly equation through in the limit of → 0. This serves as one of the most crucial checks on our computation. Additionally, if we conjecture the above relations to hold beyond three loops (which could be doubted in light of recent findings [25]), then we can even predict the -independent part of the γ JJ at O(a 3 s ): The results of γ ij uniquely specify Z ij , through Eq. (3.40). We summarize the resulting expressions of Z ij below: (3.49) Z GG and Z GJ are in agreement with the results already available in the literature up to O(a 2 s ) [83] and O(a 3 s ) [92], where a completely different approach and methodology was used.

Results of UV Renormalised Form Factors
Using the renormalisation constants obtained in the previous section, we get all the UV renormalised form factors F λ β R , defined in Eq. (3.17) and Eq. (3.21), up to three loops.
In this section we present the results for the choice of the scales µ 2

Universal Behaviour of Leading Transcendentality Contribution
In [32], the form factor of a scalar composite operator belonging to the stress-energy tensor super-multiplet of conserved currents of N = 4 super Yang-Mills (SYM) with gauge group SU(N) was studied to three-loop level. Since the theory is UV finite in d = 4 space-time dimensions, it is an ideal framework to study the IR structures of amplitudes in perturbation theory. In this theory, one observes that scattering amplitudes can be expressed as a linear combinations of polylogarithmic functions of uniform degree 2l, where l is the order of the loop, with constant coefficients. In other words, the scattering amplitudes in N = 4 SYM exhibit uniform transcendentality, in contrast to QCD loop amplitudes, which receive contributions from all degrees of transcendentality up to 2l. The three-loop QCD quark and gluon form factors [31] display an interesting relation to the SYM form factor. Upon replacement [98] of the color factors C A = C F = N and T f n f = N/2, the leading transcendental (LT) parts of the quark and gluon form factors in QCD not only coincide with each other but also become identical, up to a normalization factor of 2 l , to the form factors of scalar composite operator computed in N = 4 SYM [32].
This correspondence between the QCD form factors and that of the N = 4 SYM can be motivated by the leading transcendentality principle [98][99][100] which relates anomalous dimensions of the twist two operators in N = 4 SYM to the LT terms of such operators computed in QCD. Examining the diagonal pseudo-scalar form factors F G g and F J q , we find a similar behaviour: the LT terms of these form factors with replacement C A = C F = N and T f n f = N/2 are not only identical to each other but also coincide with the LT terms of the QCD form factors [31] with the same replacement as well as with the LT terms of the scalar form factors in N = 4 SYM [32], up to a normalization factor of 2 l . This observation holds true for the finite terms in , and could equally be validated for higher-order terms up to transcendentality 8 (which is the highest order for which all three-loop master integrals are available [101]). In addition to checking the diagonal form factors, we also examined the off-diagonal ones namely, F G q , F J g , where we find that the LT terms these two form factors are identical to each other after the replacement of colour factors. However, the LT terms of these do not coincide with those of the diagonal ones.

Hard Matching Coefficients in SCET
Soft-collinear effective theory (SCET, [102][103][104][105][106][107][108]) is a systematic expansion of the full QCD theory in terms of particle modes with different infrared scaling behaviour. It provides a framework to perform threshold resummation. In the effective theory, the infrared poles of the full high energy QCD theory manifest themselves as ultraviolet poles [109][110][111], which then can be resummed by employing the renormalisation group evolution from larger scales to the smaller ones. To ensure matching of SCET and full QCD, one computes the matrix elements in both theories and adjusts the Wilson coefficients of SCET accordingly. For the on-shell matching of these two theories, the matching coefficients relevant to pseudo-scalar production in gluon fusion can be obtained directly from the gluon form factors.
The UV renormalised form factors in QCD contain infrared (IR) divergences. Since the IR poles in QCD turn into UV ones in SCET, we can remove the IR divergences with the help of a renormalisation constant Z A,h g , which essentially absorbs all residual IR poles and produces finite results. The result is the matching coefficient C A,eff g , which is defined through the following factorisation relation: where, the UV renormalised form factor F A g R , is defined as The parameter µ h is the newly introduced mass scale at which the above factorisation is carried out. For the UV renormalised form factors [F A g ] R in Eq. (4.1), we fixed the other scales as µ 2 R = µ 2 F = µ 2 h . Upon expanding the Z A,h g and C A,eff g in powers of a s as and utilising the above Eq. (4.1), we compute the Z A,h g,i as well as C A,eff g,i up to three loops (i = 3). Demanding the cancellation of the residual IR poles of F A g R against the poles of (Z A,h g,i ) −1 , we compute Z A,h g,i which comes out to be

+
In the above expressions, L = ln Q 2 /µ 2 h = ln −q 2 /µ 2 h . These matching coefficients allow to perform the matching of the SCET-based resummation onto the full QCD calculation up to three-loop order.
Before ending the discussion of this section, we demonstrate the universal factorisation property fulfilled by the anomalous dimension of the Z A,h g which is defined through the RG equation a i s (µ 2 h )γ A,h g,i (Q 2 , µ 2 h ) . (4.7) By explicitly evaluating the γ A,h g,i using the results of Z A,h g (Eq. (4.4)) up to three loops (i = 3), we find that these satisfy the following decomposition in terms of the universal factors A g,i , B g,i and f g.i : which is in complete agreement with the existing results [112] upon identifying γ V = B g.i + 1 2 f g,i . (4.10)

Conclusions
In this paper, we derived the three-loop massless QCD corrections to the quark and gluon form factors of pseudo-scalar operators. Working in dimensional regularisation, we used the 't Hooft-Veltman prescription for γ 5 and the Levi-Civita tensor, which requires nontrivial finite renormalisation to maintain the symmetries of the theory. By exploiting the universal behaviour of the infrared pole structure at three loops in QCD, we were able to independently determine the renormalisation constants and operator mixing, in agreement with earlier results that were obtained in a completely different approach [83,92]. The three-loop corrections to the pseudo-scalar form factors are an important ingredient to precision Higgs phenomenology. They will ultimately allow to bring the gluon fusion cross section for pseudo-scalar Higgs production to the same level of accuracy that has been accomplished most recently for scalar Higgs production with fixed order N 3 LO [15] and soft-gluon resummation at N 3 LL [75,77,78,80].
With our new results, the soft-gluon resummation for pseudo-scalar Higgs production [79,80] can be extended imminently to N 3 LL accuracy [81], given the established formalisms at this order [75,78]. With the derivation of the three-loop pseudo-scalar form factors presented here, all ingredients to this calculation are now available. Another imminent application is the threshold approximation to the N 3 LO cross section [81]. By exploiting the universal infrared structure [78], one can use the result of an explicit computation of the threshold contribution to the N 3 LO cross section for scalar Higgs production [113] to derive threshold results for other processes essentially through the ratios of the respective form factors (which is no longer possible beyond threshold [15,114], where the corrections become process-specific), as was done for the Drell-Yan process [115] and for Higgs production from bottom quark annihilation [97]. These will be investigated in a separate publication.