Three-loop helicity amplitudes for four-quark scattering in massless QCD

We compute the three-loop corrections to the helicity amplitudes for $q\bar{q}\to Q\bar{Q}$ scattering in massless QCD. In the Lorentz decomposition of the scattering amplitude we avoid evanescent Lorentz structures and map the corresponding form factors directly to the physical helicity amplitudes. We reduce the amplitudes to master integrals and express them in terms of harmonic polylogarithms. The renormalised amplitudes exhibit infrared divergences of dipole and quadrupole type, as predicted by previous work on the infrared structure of multileg scattering amplitudes. We derive the finite remainders and present explicit results for all relevant partonic channels, both for equal and different quark flavours.


Introduction
The success of the collider particle physics program, whose main player today is the Large Hadron Collider (LHC) at CERN, relies heavily on our ability to model with high precision and accuracy the scattering of high energetic protons in Quantum Chromodynamics (QCD). Thanks to asymptotic freedom and the factorization properties of QCD, this intrinsically non-perturbative problem can be treated with perturbative methods, supplemented by non-perturbative information about the distribution of partons in the proton. Within this picture, an important role is played by higher order perturbative QCD calculations, which allow for a reliable and precise description of a wide range of collider processes and observables.
Thanks to a concerted effort in the high-energy community over the last few years, it is currently possible to compute predictions for many interesting reactions to second order in the strong coupling expansion, i.e. to what is usually referred to as next-to-next-to-leading order (NNLO). This has required, on the one hand, major advances in computational techniques for multi-loop scattering amplitudes , which, notably, have recently made it possible to compute various 2 → 3 processes up to two loops in QCD . On the other hand, the use of these amplitudes to perform phenomenological studies for the relevant processes at NNLO [58][59][60][61] has required the development of so-called subtraction or slicing frameworks [62][63][64][65][66][67][68][69][70][71] to properly deal with the intricate IR divergences that appear in QCD reactions.
Beyond NNLO, predictions at third order in the perturbative couplings, i.e. at N 3 LO, are known only for a handful of important LHC processes [72][73][74][75][76][77][78][79][80]. In particular, N 3 LO results are currently available only for reactions that require at most three-point threeloop integrals. Given the remarkable success of the experimental program at the LHC, it is desirable to extend these calculations to more complex processes. A particularly interesting one is di-jet production. In fact, jets are ubiquitous at hadron colliders, so understanding their dynamics is of great interest. Moreover, di-jet production is the first massless 2 → 2 process that has a non-trivial colour structure. This makes it an ideal ground for studying the structure of perturbative QCD. For example, it is by now well known that when four or more coloured partons interact, starting at the three-loop order, non-trivial colour correlations can affect the pattern of IR divergences, generating new structures [81] beyond the standard dipole formula [82][83][84][85][86][87][88]. Also, the non-trivial colour structure may create subtle violations of the factorization framework that is at the very core of theoretical predictions at hadronic colliders [89][90][91][92]. This makes jet production at hadronic colliders an extremely interesting process to investigate at higher orders.
A key ingredient for the study of jet production at N 3 LO is provided by the virtual three-loop corrections to the scattering amplitudes for the production of two jets in massless QCD. Modulo crossings, there are three main partonic channels that need to be computed: four-gluon scattering, the scattering of two quarks and two gluons, and the scattering of four quarks. All ingredients necessary for the calculation of the two-loop QCD corrections to these processes have been known for a long time [93][94][95][96][97], which have made it possible to compute the relevant scattering amplitudes [98][99][100][101]. Also, in view of extending these calculations to three loops, results for the two loop helicity amplitudes up to order 2 have been obtained [102]. For what concerns the three loop results, instead, the relevant master integrals have been computed in ref. [103], and have then been used to obtain the first three loop results for 2 → 2 scattering amplitudes in supersymmetric theories [104,105]. More recently also the first three loop corrections to the production of two photons in full QCD have been obtained [106].
In this paper, we move one step further and consider one of the three classes of partonic processes listed above, namely the scattering of four massless quarks. This particular process is interesting not only because it allows us, for the first time, to check the full structure of IR divergences at three loops in QCD, but also because it involves two external spinor structures. In fact, this property makes the use of the standard form factors method for the calculation of the helicity amplitudes particularly cumbersome, due to the fact that the γ-algebra does not close in d space-time dimensions. For the calculation of the helicity amplitudes we then make use of a different approach, recently described in refs [107,108], which allows us to calculate the helicity amplitudes in a simpler way, corresponding to the 't Hooft-Veltman scheme (tHV) [109] for the processes considered here. 1 In doing this, we also expose some subtleties in the usual approach to compute helicity amplitudes in tHV with the standard form factor method. The rest of the paper is organised as follows. We start in section 2 by establishing the notation for the calculation of the fundamental partonic channel qq → QQ, from which all other channels can by obtained by crossing. We continue in section 3, where we describe the colour and tensor decomposition of the scattering amplitude, and show how to efficiently compute the helicity amplitudes without considering evanescent Lorentz structures. In section 4 we provide details on our computational set up, and in section 5 we discuss the renormalisation and the infrared structure of the threeloop scattering amplitudes. In section 6 we discuss our final results for the main partonic channel. In section 7 we then explain how to obtain all other partonic channels from our calculation, both for quarks of equal and different flavour. Finally we conclude in section 8. In appendices A and B we provide some details on the structure of infrared divergences up at three loops, in particular focusing on the explicit derivation of the quadrupole terms, which appear for the first time with the scattering of at least four coloured partons at three loops. In appendix C, we review the analytical continuation of the amplitudes to different regions of phase space.

Kinematics and Notation
We consider the massless quark-antiquark scattering process where q and Q are in general differently-flavoured quarks and p 2 1 = p 2 2 = p 2 3 = p 2 4 = 0. From the master process in eq. (2.1) one can obtain all 2 → 2 quark scattering amplitudes with arbitrary initial states, including equal-flavour scattering q = Q. We discuss in detail how this can be achieved in section 7. The minus signs appearing in the final-state momenta imply that all momenta are taken to be incoming, This choice is convenient when performing the required crossings to obtain all the relevant partonic channels. The kinematics of the process eq. (2.1) can be parametrized in terms of Mandelstam invariants, defined as We also find it convenient to define the dimensionless ratios The variables y, z are not needed for the computation of process (2.1), but will be convenient to describe all the other processes which can be derived from it using crossing symmetry. In terms of these variables the physical scattering region is given by s > 0 , t, u < 0 which imply 0 < x < 1 . For our calculation, we work in QCD with n f massless quarks and N c colours. We denote the generators of the fundamental representation of SU (N c ) by (T a ) ij , with Tr[T a T b ] = 1 2 δ ab . Also, we indicate with C F and C A the quadratic Casimir constants. They are defined through (T a ) ij (T a ) jk = C F δ ik and f acd f bcd = C A δ ab where f abc are the structure constants. These definitions imply for the algebra of SU (N c ).

Lorentz Decomposition
We compute the scattering amplitude for the process (2.1) up to three loops in QCD, i.e. up to order O(α 4 s,b ) , where α s,b is the bare strong coupling constant. We find it convenient to write the amplitude as where we have made the external quarks colour indices i j explicit. The colour space for our process can be spanned using the two tensor structures 2 By decomposingĀ i 1 i 2 i 3 i 4 with respect to the C j , j = 1, 2 basis, we can write it as a vector in colour space:Ā where we indicate all colour-space vectors in boldface. It is important to notice that the components of the colour vectorĀ are independently gauge invariant, since a gauge transformation cannot mix different colour contributions to the amplitude. From now on we will work in the vector notation defined above, unless explicitly stated. 2 We point out that C1 and C2 are not orthogonal.
Turning to the spin structure of the process, we further decompose the scattering amplitude in terms of in terms of a basis of Lorentz structures ("tensors") T ī where the F i are scalar form factors that only depends on the Mandelstam invariants and N L is the number of elements of the basis of Lorentz structures. In our calculation, we employ dimensional regularization to deal with ultraviolet (UV) and infrared (IR) divergences. This makes the decomposition eq. (3.4) subtle. Indeed, if one works in Conventional Dimensional Regularisation (CDR), one finds that, since the γ-algebra in d dimensions does not close, the number N L of independent structures depends on the loop order [100]. However, since we are ultimately interested in computing helicity amplitudes in four dimensions, we find it convenient to work in a scheme, where we can ignore evanescent Lorentz structures right from the start. In this approach, it is possible to show that the number of Lorentz structures which are physically relevant is the same at any number of loops (N L = N ), and it equals the number of independent helicity amplitudes [107,108]. Specifically, we consider all internal momenta and polarizations in d dimensions, but restrict momenta and polarizations of the external quarks to a 4-dimensional subspace. A convenient choice for the two independent Lorentz structures describing our process is [108]: These two Lorentz structures are then sufficient at any loop order. In order to isolate the form factors from the rest of the amplitude, we define tensor projectors P i satisfying where the dot products indicates the sum over the polarisations of the external quarks, P i · T j = pol P i T j . It then follows from eq. (3.4) that P i ·Ā = F i . For the choice eq. (3.5), the explicit form of the projectors is (3.7) We recall here that the main advantage of working with scalar form factors F i is that, by construction, they only contain scalar integrals, because all the Lorentz tensor structure has been factorised out by the basis tensors T i .

Helicity Amplitudes
Ultimately, we are interested in computing the helicity amplitudes A λ 1 λ 2 λ 3 λ 4 , for the process in eq. (2.1), where we indicate with λ j the helicity of the (anti)particle with momentum p j . Since the quarks are massless, helicity is conserved along the quark lines and there are only four different possibilities that we need to consider Moreover, the symmetries of the process allow us to compute only the first two helicity amplitudes (+, −, +, −), (+, −, −, +) and then obtain the other two by acting on the result with a parity transformation, which flips the signs of the external helicities: In what follows, we will focus on the two independent configurations (λ 1 , λ 2 , λ 3 , λ 4 ) = (+, −, +, −), (+, −, −, +). We adopt the following definition for helicity states of spin-1 2 fermions where we use the well known spinor-helicity formalism [21], and indicate with ± the projection of the (anti-)particle spin along its four-momentum. Using this notation and starting from the general structure of the amplitude given in eqs (3.4) and (3.5), we obtain where In eqs. (3.12) we have introduced an explicit label for the process that we are considering, which will turn out to be useful later on when we describe the other partonic channels.
Since we expect that final analytic results for the helicity amplitudes should display the maximum degree of simplicity, in what follows we focus directly on the two linear combinations H 1 and H 2 . We write their expansion in terms of the bare coupling α s,b as (3.14) In the next section, we discuss the computation of H 1 and H 2 up to three loops in QCD.

Computation
We perform our calculations in dimensional regularization with d = 4−2 dimensions for all internal momenta and gluon fields. UV and IR singularities will then manifest themselves as poles in the dimensional regulator . In order to compute the helicity amplitudes H 1 and H 2 , we begin by producing all relevant Feynman diagrams for the process in eq. (2.1) with QGRAF [112]. Only 1 diagram contributes at tree level, 9 diagrams at one loop, 158 diagrams at two loops and 3584 at three loops. We give a few representative samples of the three-loop diagrams in figure 1. We use FORM [113] to apply the tensor projectors of eq. (3.7) to the diagrams, perform the Dirac traces and the colour algebra. The latter can be boiled down to a repeated application of the identities After the colour algebra has been performed, the quark colour indices can only appear via the two independent rank-4 tensor structures defined in (3.2). They appear in the amplitude accompanied by coefficients of the type Terms in the amplitude with different a and b are separately gauge invariant so there cannot be any gauge cancellations among them. Because of this, we compute them separately, which allows us to deal with smaller expressions. After performing the colour and Dirac algebras we can express the helicity amplitudes as linear combinations of scalar Feynman integrals with rational coefficients depending on the Mandelstam invariants s, t and the dimensional regulator . At L loops, we write the integrals appearing in the amplitudes as where γ E ≈ 0.5772 is the Euler constant and µ 0 is the dimensional regularisation scale.
Here the factor e L γ E is purely conventional and it is chosen for later convenience, while the factor µ 2L 0 ensures that the integrals have integer mass dimension. In general, for a given process with E independent external momenta and L loops one needs L(L + 1)/2 + LE independent denominators to describe all possible scalar products of loop momenta with loop or external momenta. In our case E = 3, and therefore we need 4 denominators at one loop, 9 denominators at two loops and 15 at three loops. A specific complete set of Family PL NPL {D i } at a given loop order is usually referred to as an "integral family". Having in mind the calculation of the three-loop scattering amplitudes, it is useful to organise the relevant integrals in as few integral families as possible, up to permutation of the external momenta. While one family is sufficient at one loop, we need one planar family and one non-planar one at two loops and one planar and two non-planar ones at three loops. We report them in Tabs (1), (2) and (3) . There, we indicate the loop momenta with k 1 , k 2 and k 3 . We name PL the families corresponding to planar graphs and NPL, NPL1, NPL2 the ones corresponding to non-planar graphs. The different number of loop momenta for the planar topologies eliminates any ambiguity in the naming convention.
As it is well known, the integrals in eq. (4.3) are not all independent and, instead, various types of relations can be established among them, most notably by the algorithmic exploitation of symmetry relations and integration-by-parts identities [1,2]. The latter in particular allow one to reduce the number of independent integrals by potentially several orders of magnitudes and to express the amplitude in terms of a relatively small number of so-called master integrals. While the reduction to master integrals for 2 → 2 massless scattering up to two loops can be performed very easily with automated tools, going one order higher involves a considerable increase in complexity. In practice, at three loops we proceed as follows. Once the amplitude has been expressed in terms of scalar integrals, we first use Reduze 2 [114,115] to find trivially vanishing integrals, non-trivial symmetry relations among the various integrals and corresponding ones obtained by permutations of the external momenta. This step already reduces the size of the expressions significantly as well as the number of integrals within them. We then use Finred, an in-house Table 3. Planar and non-planar 3-loop integral families.
implementation of the Laporta algorithm [116] augmented by the use of finite-field arithmetics [9][10][11]117] and syzygy-based techniques [4-7, 118, 119], in order to solve the system of integration-by-parts identities satisfied by the remaining integrals. This allows us to express the three-loop helicity amplitudes for this process in terms of the basis of master integrals computed in ref. [103]. In this way we find that the physical three-loop scattering amplitudes can be expressed in terms of the same 486 master integrals necessary for the calculation of diphoton production at three loops [106].

Ultraviolet and Infrared Subtraction
The result of the computation described in the previous section are the divergent helicity amplitudes for the process (2.1) to O(α 4 s,b ). In the following, we describe the UV renormalisation and IR subtraction of the divergent amplitude for this process.

Ultraviolet Renormalisation
We work in massless QCD with an arbitrary number of fermion flavours n f . We adopt the standard MS renormalisation scheme, where the bare coupling α s,b is written in terms of the renormalised coupling α s (µ) in the following way , µ is the renormalisation scale (for the rest of the paper we set µ 0 = µ) and The β-function coefficients are defined through where in this equation α s ≡ α s (µ). To the relevant order, they read By inserting eq. (5.1) in the α s expansion for the helicity amplitudes (3.14), we obtain the renormalised helicity amplitudes These are free from UV poles, but they still contain poles in of IR origin. We discuss how to subtract them in the next subsection.

Infrared Subtraction
While the structure of IR singularities of scattering amplitudes in massless QCD up to twoloop order has been known for a long time [120], its generalisation to three-and higher-loop order has been understood only more recently, in particular in the case where four or more coloured partons participate to the scattering process [81,85,87]. In particular, it has been shown that IR singularities are in one-to-one correspondence to the UV poles of operator matrix elements in SCET [85,87]. Therefore, UV renormalisation in SCET corresponds to IR subtraction in QCD and one can write the finite remainder of the scattering amplitude by means of a multiplicative colour-space operator Z as where {p} stands for the dependence on the external kinematics. We point out that, since we are working with vectors in colour-space as defined in Section 3, the Z colour operator can be represented as a 2 by 2 matrix which mixes the colour structures defined in eq. (3.2). Solving a renormalisation group equation one finds that Z can be rewritten as with P the path-ordering symbol, i.e. operators are ordered from left to right with decreasing values of µ . Following the notation of [81], where Γ was first computed up to three loops, the anomalous dimension operator for 4 coloured external particles is written as Above, Γ dipole represents the well known dipole colour correlations between two coloured external legs, namely with T a i the i-th particle SU (N c ) generator and from now on we use the shorthand α s = α s (µ) to indicate the renormalised coupling at scale µ. The constants γ cusp [121][122][123][124][125][126][127] and γ i [128][129][130][131] are given in Appendix B. It is also useful to define the expansions The operator ∆ 4 appears instead for the first time at three loops and contains quadrupole colour correlations among all four external legs. It can also be expanded in α s as In terms of these quantities one can organise the IR poles of the helicity amplitudes as where the IR subtraction operators read 20) with together with the quadrupole contribution ∆ 4 , which is non-diagonal in colour space and reads explicitly .

(5.25)
We stress here that, even if their x dependence has been left implicit for clarity, D 1 and D 2 are non trivial functions of the kinematics. More explicitly, the abbreviations in (5.25) read where D 1 and D 2 are expressed in terms of generalized polylogarithms of argument x with letters 0 and 1. In eqs. (5.27,5.28) we have suppressed the x dependence and used a compact index notation for the harmonic polylogarithms similar to [25,132], and we use this notation also to present our explicit analytical results for the finite remainders. We note here that both D 1 and D 2 are of uniform transcendental of weight five and are independent of the matter content of the theory, i.e. they do not depend explicitly on n f . More details on the derivation of these formulas are provided in Appendix A.

Results
After the ultraviolet and infrared pole subtraction described in the previous section, we arrive at the main result of this paper, the fully analytical expressions for the finite remainders of the helicity amplitudes for process (2.1). As previously mentioned, the helicity amplitudes for other 2 → 2 quark processes with different initial states, including the equal-flavour case q = Q, can be obtained by a combination of analytical continuation and momenta renaming from the ones for our main process (2.1). We discuss this in more detail in Section 7. We provide all finite remainders in electronic format as ancillary files attached to the arXiv submission of this manuscript.

Checks
We have performed various checks on our results. First of all, we have verified that the IR poles of our scattering amplitudes follow the pattern predicted in refs [81,85,87] up to three loops. We have then checked the finite part of our one loop amplitudes for all different partonic channels against the automated one-loop generator OpenLoops [133,134]. Finally, we have checked our one-and two-loop amplitudes through to order 4 and 2 , respectively, against the results presented in refs [100,102]. In order to successfully perform this check, one has to pay particular attention when comparing the amplitudes before IR subtraction, due to a subtlety in the dimensional regularisation scheme used in refs [100,102]. This is due to the fact that the tensor structures used to decompose the scattering amplitude in those references contain an explicit dependence on the dimensional-regulation parameter , even if the external states are taken to be in four dimensions, as in the 't Hooft-Veltman prescription. While ignoring this dependence does not change the finite remainder of the scattering amplitudes after UV and IR poles have been subtracted, it does change the bare results. We illustrate this point explicitly for the one loop case. In refs [100,102], the following four tensors are used to decompose the amplitude at one loop in CDR 3 3) where we stress that this decomposition is loop dependent and is only sufficient up to one-loop order, L ≤ 1Ā Importantly, in the 't Hooft-Veltman scheme the vector indices of the γ matrices in eqs. (6.1)-(6.4) that are not explicitly contracted with four-dimensional vector fields are in general to be taken in d dimensions. While this makes no difference for the first two tensors, the second two (6.3,6.4) depend explicitly on d and are responsible for an ambiguity in the way the helicity amplitudes are defined. Let us consider for example the fourth tensors, eq. (6.4). If we consider the γ matrices to carry d dimensional vector indices we can split the four dimensional part from the -dependent one by writing γ µ d = γ µ 4 + γ µ −2 and then use the equation as done in ref. [135] to extract the (−2 )-dimensional dependence of the γ-strings. All traces can then be evaluated as usual using the Clifford algebra relation {γ µ , γ ν } = 2g µν for d-dimensional indices µ, ν. In this case, the dimensional splitting procedure simply amounts to dependent coefficients of the 4-dimensional γ-strings. For instance, taking the first string of T 4 we find where Similarly, by repeating the same exercise for both helicities as we did in eqs. (3.12), but this time starting from the tensor decomposition in eq. (6.5), we find It is instructive to compare these formulas to the corresponding ones obtained in our approach, see eq. (3.13). Our expressions for the helicity amplitudes, despite not displaying any explicit dependence on the parameter , are exact in the 't Hooft-Veltman scheme. Notice, in particular, that the two tensors onto which we decompose the amplitude, T 1 and T 2 in eq. (3.5), have been chosen such that it makes no practical difference whether the γ algebra to fix the helicity amplitudes is performed 4 or in d = 4 − 2 dimensions. Upon substituting the form factors provided in refs [100,102] in eqs. (6.9) and (6.10) (and in the corresponding generalisations for the two-loop corrections), we find perfect agreement up to weight six with our results for the bare helicity amplitudes. We stress, nevertheless, that the results for the bare helicity amplitudes as provided in refs [100,102] are obtained by setting = 0 in the coefficients of the form factors of (6.9) and (6.10) before substituting the results for the form factors. This amounts to having assumed that the γ matrices in eq. (6.5) are purely four-dimensional. This produces a difference for the bare amplitudes with respect to ours of order O( ) at one loop and O(1/ ) at two loops. 4 However, it is easy to see that, as long as this choice is made consistently to all orders, one obtains the same results for the finite remainders in d = 4. In fact, one can imagine to first subtract UV and IR poles at the level of the individual form factors F j and, only afterwards, substitute the finite form factors in eqs. (6.9) and (6.10), and fix = 0. It is then obvious that the finite remainder cannot depend on the -suppressed contributions in eqs. (6.9) and (6.10). We have verified the last statement directly, finding perfect agreement with refs [100,102] at the level of the finite remainders.

Numerical Evaluation
We present numerical results for the finite form factors defined in (5.16), (5.17) and (5.18) calculated in the physical region 0 < x < 1 for the process in eq. (2.1), qq → QQ. To evaluate our results numerically we made use of the Mathematica package PolyLogTools [33], which in turn uses the Ginac library [136][137][138]. For the various parameters we use the following values: Here, the colour index [1] is related to the colour structure δ i 1 i 4 δ i 2 i 3 while the index [2] refers to the coefficient of δ i 1 i 2 δ i 3 i 4 . Lastly, the index (L) refers to the number of loops of the corresponding amplitude.

Crossed Channels and Equal Flavour Amplitudes
In the previous sections, we presented the computation of the helicity amplitudes for the scattering of four quarks of different flavour as in eq. (2.1). We discuss here how to use this result to derive the helicity amplitudes for all other 2-to-2 quark scattering processes, both for different and equal flavours. We start by listing all of these processes.
In the following we will consider only two helicity amplitudes for each process, since, as already discussed, the other two can be obtained by a parity transformation 5 which flips the signs of helicities. We start with the case q = Q. From the helicity amplitudes (3.12), one can find those for (7.2) and (7.3) by performing appropriate crossings of particles from the initial to the final state. This can be achieved by analytically continuing the results for (7.1) to the appropriate kinematical region (see figure 5 for reference) and then renaming the momenta of the particles involved in the crossing to match the definitions of the processes above. Amplitudes for the process (7.4) are then simply found by acting with a charge conjugation transformation 6 on (7.2). Details on how to perform the analytic continuation of the transcendental functions appearing in the scattering amplitudes are described for example in [139]. However, for convenience of the reader, we outline the method in Appendix C. In the following we will abandon the colour-vector notation for helicity amplitudes as the crossing of particles also changes the indices of the colour basis (3.2). Instead, we adopt the more conventional notation This way the crossing acts on the colour structure of the amplitude by simply exchanging indices in C 1 and C 2 . With this notation and following the procedure described above, we find 7 for processes (7.2) and (7.4) We now turn to the processes with q = Q. To obtain the amplitude for (7.5), we can use the fact that the Feynman diagrams contributing to this process are exactly the sum of the ones of processes (7.1) and (7.3). Accounting for a relative minus sign between the two contributions to the amplitude due to the exchange of two identical fermions, we write, similarly to [100],Ā qq→qq λ 1 λ 2 λ 3 λ 4 =Ā qq→QQ λ 1 λ 2 λ 3 λ 4 −Ā qQ→Qq λ 1 λ 2 λ 3 λ 4 . where in the second equation we usedĀ qQ→Qq +−−+ = 0. In the same way, we can find the amplitudes for processes (7.6) and (7.7) as a sum of those for process (7.2) and the ones obtained by crossing p 3 ↔ p 4 in (7.2). Explicitly, we obtain where we stress the exchange of helicity labels in the second term of the r.h.s. This implies where we have used thatĀ qQ→qQ +−+− = 0.
We point out that the formulas of this section are valid for the bare amplitudes as well as for the UV renormalised and the finite remainders after IR subtraction. We thus arrive at the following representation for the finite remainders of the full helicity amplitudes: where Φ process λ 1 λ 2 λ 3 λ 4 are the corresponding spinor phases, see e.g. eq. (3.12). We list the spinor phases and provide explicit analytical results for the form factors H process (L) λ 1 λ 2 λ 3 λ 4 , fin in the ancillary files.

Conclusions
In this paper we presented the analytic calculation of the three-loop helicity amplitudes for the process qq → QQ and all other related 2 → 2 quark scattering processes in massless QCD. Our results were obtained through a tensor decomposition of the scattering amplitude in the 't Hooft Veltman scheme, which allowed us to work with a minimal number of independent basis tensors, exactly matching the number of independent helicity amplitudes. This method is particularly convenient when applied to processes with multiple external fermion lines, where the number of basis tensors in conventional dimensional regularisation is known to increase with the number of loops. We employed modern integration-by-parts reductions based on finite field arithmetic and syzygy techniques to perform the demanding mapping of the three-loop integrals to a basis of master integrals. For the master integrals, analytical solutions in terms of harmonic polylogarithms were already available [103]. We point out that the finite remainders at three loops are remarkably compact with file sizes of ∼ 1 MB for each partonic channel. We provide the analytical results as ancillary files attached to the arXiv submission of this manuscript. This is the first three-loop calculation in full QCD for a process that involves the scattering of four coloured particles. In particular, our results confirm, for the first time, the structure of the colour quadrupole contribution to the IR poles of QCD scattering amplitudes with non-trivial colour flow. We performed various checks on our results, both comparing them with previous analytical calculations at one and two loops, and numerically with OpenLoops 2. The first test highlighted a subtlety in the definition of the bare helicity amplitudes in 't Hooft-Veltman scheme for processes which involve the scattering of more than one pair of fermions. In this case, using the standard projectors defined in conventional dimensional regularisation, different results can be obtained depending on whether the γ algebra to fix the helicity amplitudes is performed in 4 or in d = 4 − 2 dimensions. While no ambiguity arises at the level of finite remainders in four dimensions, extra care should be paid when comparing divergent results.
The four-quark scattering processes presented here are arguably the least complex among the ones involving four coloured particles. However, the calculation of the three-loop amplitudes for qq → gg, gg → gg and the processes related by crossings require no new concepts and can be performed following the steps described in this paper.

Acknowledgments
LT is grateful to T. Peraro for collaboration on a related project on the tensor decomposition in the 't Hooft-Veltman scheme used in this paper. We thank F. Buccioni for his assistance in comparing our one-loop results to OpenLoops. The research of FC is partially supported by the ERC Starting Grant 804394 HipQCD. AvM is supported in part by the National Science Foundation through Grant 2013859. LT and GG are supported by the Royal Society through Grant URF/R1/191125.

A Details on the IR Structure
In this appendix we provide details on the calculation of the IR subtraction operators given in eqs. (5.19), (5.20) and (5.21). We recall that one can write the finite remainder of the scattering amplitude by means of a multiplicative colour-space operator Z which we can apply similarly to the form factors in eq. (5.8): For the case of four-quark scattering considered in this paper, the Z colour operator can be represented as a 2-by-2 matrix, which mixes the colour structures defined in eq. (3.2). It is possible to write the Z operator as [85,87] where P is the path-ordering symbol, i.e. operators are ordered from left to right with decreasing values of µ . To proceed, it is useful to decompose the anomalous dimension operator according to into a dipole and a quadrupole contribution. The dipole matrix was given in eq. (5.12). We repeat its form here for convenience where α s = α s (µ), γ cusp (α s ) is the cusp anomalous dimension and γ i the anomalous dimension of the i-th external particle, which only depends on whether i is a quark or a gluon. The perturbative expansion for these quantities are given up to order O(α 3 s ) in Appendix B. In eq. (A.5), the T a i colour operators are related to the SU (N c ) generator associated with the external particle i. Specifically, following the conventions of [140] they are defined as for a initial(final) state quark(anti-quark). (A.8) We now consider how (A.5) can be written for our process (2.1) as a matrix in colour space. The operator T a 1 T a 2 acts on our colour space basis tensors as follows We change to component notation with respect to our C i basis by identifying and obtain the matrix representation Similarly, we find for the first sum in eq. (A.5) the matrix representation 1≤i<j≤4 T a i T a j γ cusp (α s ) log where t = −s−u is implied. In these equations, a small positive imaginary part is associated with the Mandelstam invariants to fix branch cut ambiguities. The dipole anomalous matrix that we have just discussed is sufficient to reproduce the IR poles of our scattering amplitude up to two loops. Starting at three loops, however, new four-parton colour correlations contained in the matrix ∆ 4 spoil the simple dipole picture. To compute ∆ 4 explicitly, we start by writing  such that, at the three-loop level, we are only interested in the first term of the expansion. By analytically continuing eq. (7) of Ref. [81] to the kinematical region with s > 0, t, u < 0, we find with x = −t/s given in eq. (2.4), the constant C in eq. (5.26), and the functions D 1 , D 2 in eqs. (5.27), (5.28). We note that all of the four partons of the kinematically dependent parts of eq. (A.17) are correlated through colour, while the kinematically independent part has a three parton colour correlation. We also note that, at this level, the quadrupole corrections to the anomalous dimension matrix can be generated solely by gluon interactions, and hence it is expected to be universal for all gauge theories involving the gluon in the particle content. Two of the representative diagrams that give rise to the quadrupole corrections are those in panels (e) and (f) of figure 1. We can highlight contributions to the quadrupole divergence, and in particular to the colour structures in the first and second line of eq. (A.17), by drawing the representative diagrams in figure 4, where the black lines represent a tree level diagram and the red lines the gluons responsible for the soft quadrupole corrections. We stress here that in reference [81] the authors give results for ∆ 4 prior to imposing momentum conservation. In this formulation, there exists a kinematical region where ∆ is purely real, which corresponds to s ij < 0 for all i, j. Consequently, we first analytically continue to the region in which s = s 12 > 0 and then impose momentum conservation.
Just like for the dipole case it is useful to represent the quadrupole operator as a matrix acting in colour space. This is immediate to do for the first line of ∆ With the dipole and quadrupole contributions at hand, it is easy to compute the Z opperator using eq. (A.3) and to perform the exponentiation as a series expansion in α s . To do so, we follow ref. [87] and introduce the short-hand T a i T a j γ cusp (α s ) = −4C F γ cusp (α s ), (A. 18) and the α s expansions (A. 19) In terms of these quantities, it is easy to find log Z = α s 4π and of its inverse Z −1 , The explicit form of the Z i and I 1 coefficients are presented in eqs. (5.19-5.24).

B Anomalous Dimension Coefficients
In this appendix, we give the universal cusp anomalous dimension γ cusp and the quark anomalous dimension γ q = γq up to O(α 3 s ) relevant for the three loop calculation of this paper. We define the α s expansion coefficients according to We then have [122,123]  and for the quark anomalous dimension [129,130] γ q 0 = −3C F , where we have set T F = 1 2 .

C Analytic Continuation
In this appendix, we briefly outline the analytical continuation of the scattering amplitudes needed to obtain the crossed channels presented in section 7. The kinematic regions and the connecting paths used in the following are shown in figure 5.
Process qQ → qQ: To derive the partonic channel (7.2) from our calculation of process (7.1), one must effectively continue our results from region I to region II, where p 1 and p 3 are incoming, see figure 5. After the analytic continuation, one can then rename p 2 ↔ p 3 . The connecting path crosses two branch cuts, one at t = 0 and one at s = 0. After crossing each branch cut, we wish to maintain the amplitude written in terms of explicitly real HPLs. This is achieved by the change of variable x → −x a − iδ for the first branch cut and x a → −1/y − iδ for the second, where the variables used have been defined in eqs. (2.4). We stress here, that the sign of the infinitesimal imaginary parts are Process qQ → qQ: To obtain process (7.3) from process (7.1), we move from region I to region III, and then exchange p 1 with p 4 . This time, we first cross the branch cut at u = 0, and then the one at s = 0. It is convenient to perform a preliminary change of variables x = 1 − x prior to crossing the cuts. We can then use exactly the same changes of variables as for the previous case: to cross the u = 0 branch cut we use x → −x b − iδ, while for the s = 0 branch cut we use x b → −1/z − iδ, such that 0 < z < 1 in region III. We can then perform the renaming p 2 ↔ p 4 , z = x| p 2 ↔p 4