Soft Theorem to Three Loops in QCD and N = 4 Super Yang-Mills Theory

: The soft theorem states that scattering amplitude in gauge theory with a soft gauge-boson emission can be factorized into a hard scattering amplitude and a soft factor. In this paper, we present calculations of the soft factor for processes involving two hard colored partons, up to three loops in QCD. To accomplish this, we developed a systematic method for recursively calculating relevant Feynman integrals using the Feynman-Parameter representation. Our results constitute an important ingredient for the subtraction of infrared singularities at N 4 LO in perturbative QCD. Using the principle of leading transcendentality between QCD and N = 4 super Yang-Mills theory, we determine the soft factor in the latter case to three loops with full-color dependence. As a by-product, we also obtain the finite constant f (3)2 in the Bern-Dixon-Smirnov ansatz analytically, which was previously known numerically only.

A remarkable property of gauge theories is that scattering amplitude containing a soft gauge boson can be factorized into a universal soft factor S µ (k) and a hard scattering amplitude with the soft gauge boson removed, This is known as the soft theorem [1][2][3].Note that for non-abelian gauge theory such as QCD, the soft factor is an operator acting on the color space of n-point amplitude.The soft theorem has found many applications in high energy physics both phenomenologically and theoretically.
In this paper we present a calculation of the soft theorem for soft gluon radiation from two hard partons to three loops in the perturbative QCD.Our main motivation is from the intimate relation between soft theorem and infrared behavior of QCD in higher order perturbation theory [4][5][6][7][8].In particular, the precision knowledge of the soft theorem in QCD is essential for constructing infrared subtraction terms in fixed order calculation [9][10][11][12][13][14][15][16][17][18].It also contributes to the calculation of various soft function in Soft-Collinear Effective Theory (SCET) [19][20][21][22][23].In this paper, we focus on the calculation of the soft factor with two hard-scattering partons and a single soft gluon emission.While it is not the most general soft factor beyond one-loop, it suffices for applications to some of the most important processes, such as the Drell-Yan process, e + e − to dijet, and 1+1 jet production in deepinelastic scattering.The one-loop contribution of it was calculated more than two decades ago [24][25][26][27][28][29][30].The two-loop soft factor was initially extracted from the soft limit of splitting amplitude up to O(ϵ 0 ) in dimensional regularization parameter [31], and was obtained through to O(ϵ 2 ) and beyond, either by direct calculation in SCET [32], or by extracting from amplitude [33].The two-loop soft factor constitutes an essential contribution to the total cross section of Higgs boson production at N 3 LO in the threshold limit [34][35][36][37].The two-loop soft factor is also an important ingredient for constructing infrared subtraction for generic perturbative QCD calculation at N 3 LO.Besides the two-loop soft factor for single gluon emission, also relevant are the one-loop double-parton soft emission [38][39][40], and the tree-level triple-parton soft emission [41][42][43].In addition, starting from two loops, a nontrivial color structure that correlates more than two partons first arises and has been computed in [44].To further push the theoretical accuracy towards the N 4 LO frontier for scattering cross section, in this paper we perform the calculation of the single soft-emission soft factor with two hard partons at three loops through O(ϵ 2 ).
To facilitate the calculation, we have developed a systematic approach to calculate single-scale soft integrals using Feynman parameter representation.Our main idea is to introduce an auxiliary scale in the parametric representation and directly construct differential equations in the parametric representation with respect to this scale.A parametric integral is nothing but a multi-fold integral.Thus, an auxiliary scale can trivially be introduced by leaving one fold of integration untouched.The obtained integrals can be calculated by using the standard differential-equation method [45,46].The boundary conditions of the differential-equation system can again be expressed in terms of parametric integrals, which can be calculated by further applying this method.Thus, this method allows us to calculate Feynman integrals recursively until the boundary conditions can be trivially determined.
Besides phenomenological interests, the soft factor is also useful in determining quantities of theoretical interests.For instance, since the soft factor can be understood as the soft limit of the corresponding full amplitude, it shares the same iterative structure of the full amplitude in the maximally supersymmetric N = 4 Yang-Mills theory (MSYM) [47,48].It was conjectured by Bern, Dixon, and Smirnov in ref. [48] that the planar maximally helicity violating (MHV) amplitudes in MSYM can be obtained iteratively.Specifically, the l-loop planar MHV n-point amplitude in MSYM is determined by the one-loop amplitude up to some kinematic-independent constants, which are known to three loops numerically [49].Assuming the principle of transcendentality [50], we obtain the soft factor in MSYM from reading off the leading transcendental part of the QCD results.We obtain the analytic expression for the three-loop constant f (3) 2 in the BDS ansatz, which agrees well with the previously numerically determined one [49].In addition, we also predict the full-color dependence for the soft function in MSYM at three loops, which provides a test to the three-loop non-planar form factor of 1 → 3 decay in MSYM [51], once the relevant master integrals there are computed.
The rest of this paper is organized as follows: in sec.2, we describe the method to calculate the soft factor based on an effective theory.The result is expressed in terms of single-scale soft master integrals.In sec.3, we develop a systematic method to calculate these master integrals recursively based on the differential-equation method, with all the boundary integrals evaluated to gamma functions.The final results in both QCD and SYM are presented in sec.4.

Calculation of QCD soft theorem to three loops
In this section we introduce the method for constructing the integrand for loop-level soft factor in QCD.Our approach is based on SCET, where the soft factor can be expressed as a transition matrix element of soft Wilson lines from vacuum to single gluon state.We use this definition to construct the integrand through three loops.

Soft theorem from Soft-Collinear Effective Theory
The soft theorem with an outgoing soft gluon can be well described by Wilson lines in Soft-Collinear Effective Theory (SCET) [19][20][21][22][23] (See also, for example, ref. [52] for a discussion from QCD factorization).That is, where Y k (x) is a semi-infinity Wilson line standing for the color source of an external hard parton.In this paper, we restrict ourselves to the case of two Wilson lines with m = 2.For an outgoing Wilson line, it starts from the origin and extends to null infinity, where the subscript 's' in A a s refers to the soft gluon field.Similarly, an incoming Wilson line is defined as In the above equation, the P refers to path ordering where we define A(x + sn k ) = A a (x + sn k )T a k , and T a k is the color-charge operator defined in the color space formalism [27].For an outgoing quark (incoming anti-quark), (T a k ) ij = (t a ) ij , for an outgoing anti-quark (incoming quark), (T a k ) ij = − (t a ) ji , for a gluon, (T a k ) bc = −if abc , where f abc are structure constants, and t a are the Gell-Mann matrices with the normalization Tr[t a t b ] = 1 2 δ ab .According to Lorentz and color structures, the soft factor J µ with two Wilson lines can be decomposed into the following form up to three loops: where and their coefficients C 12 , D 12 only receive contributions starting from three-loop order.
We stress that all these scalar factors don't depend on the particular representation of the Wilson lines, a form of Casimir scaling.The form of eq.(2.5) is constructed from scaling symmetry and dimensional analysis.It has been checked by an explicit computation which will be described in detail below.A related quantity, the l-loop Eikonal function can be derived from the soft factor, where in SU(N c ) group for quarks and gluons respectively.The eikonal function directly contributed to the soft-virtual cross section of Drell-Yan or Higgs production, see e.g.[35,36].Here and in the following, we always expand the quantity in with α s = g 2 s /(4π), for example, (2.9) The S (0) 12 is the well-known tree-level Eikonal function, (2.10) We are mainly interested in the higher-order corrections for scalar form factors in eq.(2.5) and the tree-level Eikonal function in eq.(2.10).The soft factor in general depends on the direction of the Wilson lines (incoming or outgoing).However, due to a rescaling symmetry n µ 1 → λ 1 n µ 1 and n µ 2 → λ 2 n µ 2 for arbitrary λ 1 and λ 2 , this dependence can be fully encoded in terms of a factor S ϵ : where ϵ = (4 − d)/2 is the dimensional regulator, λ AB in the phase factor e −iλ AB π is 1 if A and B are both incoming or outgoing, and λ AB = 0 for other cases (see for example [30]).By factoring out the dependence on S ϵ for eq.(2.5) and eq.(2.7) at each order, the remaining contributions are not sensitive to the direction of Wilson lines, In this paper, we determine r 12 and b 12 , c 12 , d 12 in above equation to three-loop order.

Construction of loop integrand for Soft theorem
To construct the loop integrand for the soft theorem, we first derive the effective Feynman rules for soft Wilson lines as shown in eq.(2.2) and eq.( 2.3).It can be conveniently done by expanding the Wilson lines order by order in g s .We get the following eikonal Feynman rules up to three gluon emissions (We checked explicitly that the Feynman rules with more gluon emissions are not needed for the computation of three-loop soft theorem), where the sign of Feynman's prescription iδ k 0 + stems from the path along the light cone of outgoing or incoming Wilson lines.The δ k is 1 for an outgoing Wilson line and δ k = −1 for an incoming Wilson line.The color-charge operator T a k is the same operator as defined in eq.(2.2) and eq.(2.3).In the above equations, plus permutations indicates a summation over all external gluon indices (simultaneous permutation of µ i , a i , k i ).Each vertex for the Wilson line above can be understood as a sum of conventional diagrams with soft gluon emissions.For example, the vertex with two soft gluons in eq.(2.13) for an outgoing Wilson line can be understood as the sum of the following two diagrams: Here the effective coupling between a soft gluon and an eikonal line is ig s T a k n µ k .And the propagator of the eikonal line with momentum k is i n k •k+i0 + .We generated in QGRAF [53] all relevant Feynman diagrams, which implement particle interactions from standard QCD and interactions due to the above effective vertices with up to three-gluon emissions.In figure 1, we show some sample Feynman diagrams.The amplitude is invariant under the rescaling of n 1 , n 2 , such that the only scale is µ 2 S (0) 12 .Therefore, the soft factor only receives contributions from one-particle-irreducible (1PI) diagrams with a number 550 at three loops.Subsequently, an in-house Mathematica code was used to substitute the Feynman rules into Feynman diagrams, and FORM [54][55][56] and Color.h[57] were used to evaluate Dirac and color algebra.To verify the (generalized) Casimir scaling principle, we use the effective Feynman rules as shown in eq.(2.13) for Wilson lines in both fundamental and adjoint representations.Regarding the topology classification, the package Apart [58] was first used to eliminate the linear dependence of propagators largely due to the multiple linear propagators from the effective vertices.After the partial fraction, we found 780 topologies that were then reduced to 160 topologies by applying a self-written code.The code implemented a simple algorithm that tries to find a loop momentum transformation relating the two topologies with each other by carefully searching all possible loop momentum transformations.We noted that a similar algorithm was also implemented in the public package Reduze 2 [59]1 .By appending some proper propagators stemming from irreducible numerators, the 160 topologies can be further mapped into 25 integral families.The definition of these integral families can be found in Appendix.A. The integration-by-parts (IBP) [60] reductions were done by Kira [61] equipped with FireFly [62], which implements the Laporta algorithm [63] as well as finite fields and function reconstruction techniques [64,65].After IBP reductions, we found 52 master integrals and only 49 of them appeared in the amplitude, and the master integrals were found to appear only in six integral families.
To check the gauge invariance of the amplitude, we use the Feynman gauge as well as the light cone gauge for the polarization summation of the external gluon as shown in eq.(2.7), For internal gluons, we use the R ξ gauge but with the amplitude truncating at We found the amplitude is indeed gauge invariant provided that six extra relations exist within the 49 master integrals.We verified explicitly these six extra relations by computing all 49 master integrals as shown in section 3.

Calculation of master integrals
In this section, we present the details of our approach to compute the single-scale soft master integrals.Our approach is based on Feynman parameter representation.Using the differential equation with respect to an auxiliary scale appearing in the intermediate step and reduction of integrals in Feynman parameter representation, we manage to compute all master integrals iteratively, with the final boundary conditions coming from simple Gamma functions.

Differential equations
We calculate the master integrals by using the differential-equation method [45,46].For soft integrals, the scale dependence is trivial.To get a nontrivial differential-equation system, we need to introduce an auxiliary scale.While there are several widely used methods to achieve this, such as the Drinfeld-associator method [66] and the auxiliary-mass-flow method [67][68][69][70], these methods are less appropriate for the calculation in this work.Because an extra mass scale on either an external line or an internal line may highly increase the complexity of the IBP reduction.A better choice is to introduce a scale that is essential to the integral to be calculated.It is evident that an extra scale can be introduced by leaving one fold of integration untouched, or equivalently, by inserting a delta function.For phase-space integrals, this can be done in the momentum space (see e.g.refs.[71,72]).While for normal loop integrals, an extra delta function in the momentum space may even complicate the calculation.A better choice is to insert a delta function in the parametric representation.
The obtained integrals can further be reduced by using the method developed in refs.[73][74][75].
We consider the calculation of the parametric integrals of the following form Here the integration measure is dΠ , with E (n) (x) a positive definite homogeneous function of x of degree n.The region of integration for x i is (0, ∞) when i > m and (−∞, ∞) when i ⩽ m.F is a homogeneous polynomial of x of degree L + 1.For integrals with momentum-space correspondences, L is the number of loops.For loop integrals, the polynomial F is related to the well-known Symanzik polynomials U and F through F = F + U x n+1 .But here we consider the more general parametric integrals which may not have momentum-space correspondences.This generalization is necessary, because some asymptotically expanded integrals may not have any momentumspace correspondence [76,77].By virtue of the homogeneity of the integrands, it can be shown that the parametric integrals satisfy the equations ) A parametric integral can be understood as a function of the indices λ i .Then we can define the following operators.R i I(λ 0 , . . ., λ i , . . ., λ n ) =(λ i + 1)I(λ 0 , . . ., λ i + 1, . . ., λ n ), It is understood that And we formally define operators ẑn+1 and xn+1 , such that ẑn+1 I = I, and We assume that xn+1 is always to the right of U (x) in F(x).By using these operators, we can write eq.(3.2) in the following form 2 : Let y be a kinematical variable, then it is easy to see that the parametric integrals satisfy the equation To get a nontrivial scale dependence, we insert a delta function into the parametric integral in eq.(3.1), and get Here the function E (n) is the one defined below eq.(3.1).Equation (3.2) also holds for eq.(3.5), since it is a consequence of the homogeneity of the integrand.In practical calculations, by a proper choice of the E (0) (x), we can eliminate one fold of integration by using the delta function.The resulting y-dependent integral is still of the form in eq.(3.1).Thus it can be reduced in the parametric representation and then calculated by using the differential-equation method.Compared with the original integral, the y-dependent integral is one less fold of integration.By successive applications of this method, we can calculate Feynman integrals recursively.
The method described in this subsection applies to both normal loop integrals and phase-space integrals.But we do not consider phase-space integrals hereafter.That is, we take m = 0 in eq.(3.1).

Rules for choosing E (0)
In principle, the function E (0) in eq.(3.5) can be chosen arbitrarily.Nevertheless, for a general choice of E (0) , it is not easy to express the right-hand sides of eqs.(3.2) in terms of regular parametric integrals.In practical calculations, we choose E (0) to be of the form Then we can eliminate the integration with respect to x i by using the delta function.That is, = dy dΠ (n) x j I (−n−1) dy For integrals that have momentum-space correspondences, this is equivalent to the method of combining two propagators with a Feynman parameter [78].The pair {x i , x j } can still be arbitrarily chosen.A good choice may greatly simplify the calculation.In this section, we provide a method to choose E (0) wisely such that the number of regions of the asymptotic expansion for the obtained y-dependent integral is minimized.Consequently, the boundary conditions of the differential equations are simplified.
A general F polynomial is of the structure where C F ,a are some x-independent constants.This polynomial may not depend on any physical scale.Thus it does not make sense to talk about asymptotic expansion for the corresponding parametric integrals.Nevertheless, we can still formally introduce the notion of "region" for this polynomial by using the idea of the convex hull described in ref. [77].Specifically, a region r is associated with a subset S r of {1, 2, • • • , A} and a n+2 dimensional vector k r , such that the number of elements of S r is not less than n + 1, and It is easy to see that Λ ai with a / ∈ S r is linearly independent of Λ ai with a ∈ S r .Since F is a homogeneous polynomial of degree L + 1, we have n+1 i=1 Λ ai = L + 1.Thus, if Λ bi = a∈Sr c ba Λ ai , we have a c ba = 1, and i Λ bi k r,i = k r,0 a c ba = k r,0 .Hence b ∈ S r .And it is easy to see that the cardinal number of S r should be smaller than A, because otherwise the corresponding parametric integral is scaleless.To see this, without loss of generality, we assume that k r,1 ̸ = 0. Then we rescale x i with i > 1 by . If S r = {1, 2, . . ., A}, the x 1 dependence of F can be factored out.Thus the integration with respect to x 1 is scaleless.
Presently we will show that those regions defined by eqs.(3.9) are intimately related to the regions of the y-dependent integrals.
We consider the integrals obtained by replacing x i with yx j .The corresponding F polynomial reads For simplicity, we formally denote y by x i for F ′ .Obviously, we have with It is easy to see that with According to the convex-hull algorithm described in ref. [77], a vector k ′ r with k ′ r,i > 0 gives exactly a region of asymptotic expansion in the limit of y → 0, because terms , by choosing the pair {i, j} such that the cardinal number of the set R ij is minimized, the number of regions of asymptotic expansion is minimized.
Obviously, after expanding the F polynomial asymptotically in a region r, only terms in S r survive.Thus, by choosing the pair {i , j} such that the cardinal number of S r (denoted by N r ) is minimized, the boundary integrals are simplified.
As a summary, we choose the pair {i, j} according to the following rules: (1) We choose the pair {i, j} such that the cardinal number of R ij is minimized, where R ij is defined in eq.(3.15).
(2) Among all the pairs satisfying the first rule, we choose the one such that max{N r |r ∈ R ij } is minimized, where N r is the cardinal number of S r .

Boundary integrals
By using the method described in the previous sections, we can construct differential equations for the parametric integrals.The boundaries of the solutions of the differential equations can further be expressed in terms of parametric integrals.Thus, this algorithm can be carried out recursively.The algorithm terminates when the F polynomial has exactly n + 1 monomials.In this case, the parametric integral can be expressed in terms of gamma functions.We have .
Here Λ ai and C F ,a are defined in eq.(3.8), and L is defined in the paragraph below eq.(3.1).
The derivation of eq.(3.16) is as follows.
We introduce a new set of variables The Jacobian is For the integration measure dΠ (n+1) , we choose n+1 .Then we have . (3.20)

Analytic continuation
While the analytic continuation is not a problem for the calculations in this paper, it needs to be considered in order to develop a general-purpose algorithm.For a F polynomial with both positive terms and negative terms, a Feynman parameter may cross a branch point in the region of integration.Generally speaking, it is not easy to determine the branch while a Feynman parameter crosses a branch point.A possible solution to this problem is as follows.We replace each negative coefficient of F, denoted by −C F ,a , by −yC F ,a , and construct differential equations with respect to y.The imaginary part of y should be i0 + due to the i0 + prescription of Feynman propagators.We determine the boundary conditions at y = 0 − .All the boundary integrals are with positive definite F polynomials and thus can further be evaluated by using the method described in previous subsections.

Examples
As an example of the application of the method described in this section, we consider the calculation of the following integral: The F polynomial for this topology reads ) .
The F polynomial for the integral family I 3 is The integral family I 3 does not have an evident momentum-space correspondence.It can further be calculated by using the method described in this section.The k r vectors for F 3 are . There are 7 pairs {x i , x j } with only one k r such that k r,i > k r,j , among which the pair {2, 6} has the minimal N r .Insertion of δ(y − x 2 x 6 ) leads to the integral , 0, 0, 0, 0, 1), with the F polynomial The integral I 4,0 can further be reduced to where the master integrals are , 0, 0, 0, 0, 0 .
The differential equations for these integrals are This differential-equation system can be converted into the canonical form [81] by using the package epsilon [82] which implements Lee's algorithm [83].The obtained differentialequation system is solved by using the standard differential-equation method.The boundary conditions can be determined by applying the method developed in this section recursively.
Here we do not go into more detail.
4 Soft theorem at three loops in QCD and N = 4 sYM In this section, we present the results for the three-loop soft factor in QCD up to O(ϵ 2 ).These results are necessary ingredients for QCD corrections or soft function calculation at N 4 LO.We also derive the corresponding soft factor in N = 4 sYM with full-color dependence, using the principle of leading transcendentality [84].Note that the principle of leading transcendentality has not been proved, but is known to work in many cases, including e.g.twist operator dimensions [50,84], soft functions or Wilson loops [36,85], form factors [86][87][88].We verify that the leading color contributions agree with a previous calculation [32] based on BDS ansatz [48], and determine a three-loop constant f analytically.

IR singularities of soft factor
Before presenting our results for soft factor at three loops, we first discuss its infrared singularities.IR singularities for scattering amplitudes have been understood to be factorized, as a result of soft-collinear factorization [89][90][91].Since the soft factor is simply the soft limit of scattering amplitude, we can extract the IR singularities of the soft factor by taking the soft limit in the IR singularities of scattering amplitude.To all orders in perturbation theory, the IR singularities of massless scattering amplitudes are governed by a multiplicative renormalization factor Z, which in general is a matrix in color space where in the equation above |M fin.,n ({p}, µ)⟩ is the IR renormalized finite amplitude, and |M n (ϵ, {p})⟩ is the UV renormalized amplitude.The IR renormalized factor is and where explicit data for the anomalous dimension will be provided in the appendix.The Mandelstam variable is defined as s ij = 2σ j p i p j + i0 with the sign factor The notation (i, j) refers to unordered pairs of distinct parton indices.In eq. ( 4.3), ∆ 3 refers to tripole contribution which is kinematics-independent and starts to contribute at the three-loop order [92]: Another contribution from the quadrupole term ∆ 4 involving four partons is also present at the three-loop order, for example in the three-loop four-parton scattering amplitude in N = 4 and QCD [93][94][95][96].However, because we only deal with the soft factor of two hard partons (only scattering amplitude involving three colored partons is required), it is not needed in this work.
The soft factor computed in this work refers to the soft gluon limit of a three-parton amplitude in QCD.We write the IR renormalization formula as where |M 3 (ϵ, p 1 , p 2 , p 3 )⟩ is the UV renormalized amplitudes with three massless QCD partons, for example γ * → q(p 1 )q(p 2 )g(p 3 ).In the soft gluon limit, the soft gluon factorization demands that the IR renormalized amplitude factorizes as where J(p 3 , µ) is the IR renormalized soft factor, and |M fin.,2 (p 1 , p 2 , µ)⟩ is the IR renormalized 2-parton amplitude, This leads to the relation The infrared singularities of the three-loop soft factor can then be read-off from Z s (ϵ, p 3 , µ).

Soft theorem to three loops in QCD
We are now ready to present our results for the soft factor to three loops.The results were verified to satisfy the (generalized) Casimir scaling principle, such that we are able to write them down in a unified form for both fundamental and adjoint representations.The corrections of B 12 in eq.(2.5) were calculated to two loops in [32,33], we list them here using the convention of eq.(2.12) for completeness.At one-loop order, the result can be expressed in terms of the following gamma functions, b For the two-loop corrections, we give the result to ϵ 4 and found full agreement with the ϵ expansion of all-order result in [33], b where several regular zeta values up to transcendental-weight 8 and one multiple zeta value are involved, The higher-order corrections of the Eikonal functions in eq.(2.7) can be readily expressed in terms of the above results, in the cases of one and two loops, and in the three-loop case, r where in gauge group SU(N c ), the quadratic color structures are evaluated to be the following explicit expressions, We emphasize that our results are for unrenormalized quantities.To perform the ultraviolet (UV) renormalization, we just need to renormalize the strong coupling constant, i.e., with , where β i is QCD beta function which will be collected in the appendix.After UV renormalization, the remaining poles stem from IR singularities.We checked that the IR poles in our explicit results agree with those predicted in section 4.1.

Soft theorem in N = 4 sYM and BDS ansatz at three loops
The soft theorem in N = 4 sYM can be easily extracted from QCD results assuming the principle of maximal transcendentality [50,84].At the one-loop order, the Eikonal function for the soft theorem in N = 4 sYM and QCD are identical, i.e., S 12 . (4.21) At two-loop order, by keeping only the maximal transcendentality part of eq.(4.12), the N = 4 sYM results up to ϵ4 reads, where only the leading color contributes.Similarly, at the three-loop order, we have where the three-loop N = 4 result receives contributions from both leading color and subleading color.The sub-leading color appears for the first time at three-loop order and comes from the fourth invariant tensor d abcd A d abcd A in (4.18) solely.While not rigorously proved, the validity of the principle of maximal transcendentality at subleading color was confirmed in the case of the four-loop Sudakov form factor [88,97]. 3 Interestingly, the soft limit of three-loop splitting amplitude in planar N = 4 sYMs can be reorganized into the following form [32] which is the soft limit of the well-known BDS ansatz [48],    previously known results obtained from BDS ansatz.In addition, we also analytically determine a three-loop constant f (3) 2 = 31ζ 2 3 + 1909ζ 6 /48 in BDS ansatz, which was only known numerically.Our new results are the full-color dependence of the three-loop soft factor, which can be used to check the three-loop form factor 1 → 3 decay once the relevant master integrals are known.Towards the full-color three-loop 1 → 3 form factor in QCD, we notice that the corresponding leading-color result became available quite recently [99].

B QCD beta function and anomalous dimensions
To predict the infrared singularities as shown in sec.4.1, we need the relevant anomalous dimensions and QCD beta function, which will be listed below for readers' convenience.
The QCD beta function is defined as and here we need to three-loop order [100] We perform a perturbative expansion for the anomalous dimension γ as follows, where a s is defined in eq.(2.8).The cusp anomalous dimension to three-loop order was first extracted from the three-loop non-singlet splitting functions [101], and is given as Finally, the quark and gluon anomalous dimensions of the three-loop order can be extracted from the three-loop quark and gluon form factors [102,103], where the explicit results can also be found in [104].

C Instructions of the supplementary material
We present our results for master integrals valid to transcendentality weight-8 and soft theorem to ϵ 2 in the supplementary material.MIsolutions.mcontains the solutions of all 49 three-loop master integrals from the six integral families as defined in section A. softTheorem.mcontains the results as shown in eq.(4.11) to eq. (4.15).

1 Introduction 1 2 2 . 1 5 3 13 4
Calculation of QCD soft theorem to three loops 3 Soft theorem from Soft-Collinear Effective Theory 3 2.2 Construction of loop integrand for Soft theorem Calculation Soft theorem at three loops in QCD and N = 4 sYM 16 4.1 IR singularities of soft factor 16 4.2 Soft theorem to three loops in QCD 18 4.3 Soft theorem in N = 4 sYM and BDS ansatz at three loops 21 0 are two light-like vectors.The form factor B 12 starts to contribute at one loop.The quadrupole invariant tensor d abcd A and d abcd F are defined by

2 Figure 1 .
Figure 1.Sample Feynman diagrams for three-loop soft theorem.From left to right, when contracting with tree diagrams with single gluon emission, the first diagram is zero, the second diagram contributes to d abcd R d abcd F , and the third diagram contributes to sub-leading color only. r

2
[74,75]that here the definition of ân+1 is slightly different from that in refs.[74,75].The reason is that, with the definition of ân+1 in refs.[74,75],eq.(3.4) is invalid if the first Symanzik polynomial U depends on y.And in this paper we do need to consider the situation where U depends on y.
n+1 I rather than xi n+1 xj n+1 I .In practical calculations, we always express xi n+1 in terms of ân+1 from the very beginning.