Three-loop soft function for energetic electroweak boson production at hadron colliders

We calculate the three-loop soft function for the production of an electroweak boson (Higgs, $\gamma$, $W^\pm$, $Z$) with large transverse momentum at a hadron collider. It is the first time a soft function for a three-parton process is computed at next-to-next-to-next-to-leading order (N$^3$LO). As a technical novelty, we perform the calculation in terms of forward-scattering-type loop diagrams rather than evaluating phase space integrals. Our three-loop result contains color-tripole contributions and explicitly confirms predictions on the universal infrared structure of QCD scattering amplitudes with three massless parton legs. The soft function is a central ingredient in the factorized cross section for electroweak boson production near the kinematic endpoint (threshold), where the invariant mass of the recoiling hadronic radiation is small compared to its transverse momentum. Our result is required for predictions of the near-threshold cross sections at N$^3$LO and for the resummation of threshold logarithms at primed next-to-next-to-next-to-leading logarithmic (N$^3$LL$^\prime$) accuracy.


Introduction
The production of an electroweak (EW) boson (Higgs, γ, W ± or Z) with sizable transverse momentum is among the most fundamental scattering processes at hadron colliders like the LHC. The corresponding cross sections are experimentally and theoretically relatively clean observables and allow precision phenomenology. Their measurements provide excellent tests of the Standard Model, and probe the parton distributions inside the colliding hadrons at small distances. EW gauge boson production at large transverse momenta also represents an important background to Higgs measurements and new physics searches. On the other hand the Higgs transverse momentum spectrum plays e.g. an important role for the analysis of the Higgs couplings.
Moderate and large transverse momenta of EW bosons in hadron collisions are generated primarily by the recoil against hard QCD radiation (jets). Corresponding tree-level processes with one hard final state parton are depicted in fig. 1. Compared to the inclusive EW boson production the leading order (LO) contribution therefore contains an additional power of the strong coupling constant α s , and it is technically harder to achieve the same level of accuracy on the theory side. In the last couple of years fixed-order QCD predictions for EW boson production with nonzero transverse momentum have reached NNLO precision [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. Depending Figure 1. Tree-level Feynman diagrams for electroweak boson production with large transverse momentum at hadron colliders. For EW gauge boson (orange wiggly lines) production example diagrams of two different production channels are shown. on the kinematic region and the type of process the NNLO corrections can be substantial. For example the NNLO QCD corrections to Higgs+1 jet production at the LHC enhance the cross section by roughly 20%, while the uncertainties from scale variations still amount to about 10% [2]. In view of future precision measurements a further reduction of the theory error would be desirable.
Whereas full QCD N 3 LO results seem currently out of reach, threshold approximations and the corresponding resummation of threshold logarithms can certainly improve the current stateof-the-art fixed-order predictions, at least at large transverse momenta. 1 In the limit, where the transverse momentum of the EW boson becomes maximal (at fixed rapidity), the invariant mass of the recoiling hadronic radiation vanishes and vice versa. Close to this threshold all final state QCD radiation is either soft or anti-collinear to the direction of the EW boson and the production cross section can be factorized into a convolution of hard, jet, and soft functions [17,18]. Within the framework of soft-collinear effective theory (SCET) [19][20][21][22][23][24] the jet and soft functions are expressed as effective operator matrix elements, while the hard function is a matching coefficient. Each of the factorization functions obeys a corresponding renormalization group equation (RGE). Solving these and evolving the functions from the respective physical hard, jet, and soft scales to a common renormalization scale systematically resums large threshold logarithms (of the ratio between the transverse momentum and the hadronic invariant mass). In this way near-threshold predictions for cross sections with NNLL and partial N 3 LL resummation along with the corresponding NNLO threshold corrections were obtained [18,[25][26][27][28][29][30].
In order to improve the predictions for EW boson production near threshold to N 3 LO and (resummed) N 3 LL accuracy 2 the factorization functions are required at three-loop order. Given the recent progress in the calculation of N 3 LO four-particle scattering amplitudes in QCD [32,33], there is hope that also the three-loop virtual QCD corrections to EW boson + jet production, constituting the hard function, become available in the not-too-far future. The relevant quark and gluon jet functions were computed at three loops in refs. [34,35]. The present paper is dedicated to the calculation of the three-loop soft function.
The soft function is given by a phase space integral constrained by the threshold measurement over a squared matrix element with one outgoing and two incoming Wilson lines along widely separated lightlike directions. Although the relevant two-loop single [36], one-loop dou- 1 Often the threshold terms amount to the bulk of the corrections at a given order and the cross sections computed in the threshold limit represent indeed good approximations also at moderate transverse momenta. 2 For details and advantages of the primed counting, see e.g. ref. [31].
ble [37], and tree-level triple emission [38] soft currents are available, performing the necessary phase space integrations is a difficult task. We therefore follow a different approach. Generalizing the 'dispersive method' we used to calculate the two-parton soft function for heavy-to-light decays in ref. [39] we express our three-parton soft function in terms of a single Wilson line correlator corresponding to a forward scattering process. The phase space integrations then translate to loop integrations of Feynman diagrams with one internal and four external Wilson lines. This trick allows us to directly apply modern multi-loop technology and makes the three-loop calculation of the soft function relatively straight-forward to carry out. Our paper is organized as follows. In sec. 2 we briefly discuss the factorization formula for the near-threshold production of EW bosons at the LHC and define the soft function employing the color space formalism. In sec. 3 we devise the 'dispersive method' to compute the soft function in terms of loop diagrams. Section 4 describes the three-loop calculation in some detail and we present our results for the soft function in sec. 5, including a discussion of its divergence and color structure. We conclude in sec. 6.

Threshold factorization and soft function definition
In order to put the calculation performed in this paper into context we consider EW boson (Y = γ, W, Z, h) production at the LHC, i.e. the process pp → Y + X. The 'hadronic' invariant mass of the final state particles in X (including the proton remnants) is where P 1,2 are the proton momenta. For M X → 0 the particles in X form a collimated jet and any additional radiation must be soft. In this threshold limit the cross section factorizes, which allows the resummation of large logarithms ∼ log(M X /E cm ) with E cm being the center of mass energy of the colliding protons [17]. Defining the kinematic invariants at partonic level aŝ we can express the invariant hadronic mass as is the total invariant mass of the final state partons (excluding the proton remnants), x 1,2 are the fractions of the proton momenta carried by the incoming partons of the hard scattering process, m Y is the mass of the EW boson, and m ⊥ = p 2 T + m 2 Y is its transverse mass. The maximal transverse mass for fixed rapidity y of the EW boson is given by (m max 3) it is apparent that M X → 0 requires m X → 0 and x 1,2 → 1 simultaneously. We also see that this limit is approached for large transverse momenta of Y .
The corresponding factorized near-threshold cross section as derived in SCET [18] for the channel ab → cY with {a, b, c} = {q,q, g}, {q, g, q}, {q, g,q}, or {g, g, g} takes the form 3 The parton distribution functions (PDFs) f i (x, µ) are evaluated close to their endpoints at x = 1. The hard function H ab corresponds to a SCET Wilson coefficient obtained from matching SCET currents to full QCD matrix elements at the hard scale ∼ E cm ∼ p T ≈ p max T and contains the full QCD virtual corrections to the hard scattering process. It depends on the hard kinematics and is known to two-loop order for Y = γ, W, Z, h [40][41][42][43]. The jet function J c describes the collinear radiation along the jet direction n µ J , which is anti-collinear to the direction of Y , 4 and only depends on the type of parton initiating the jet and its virtuality ∼ M X . The respective three-loop results were obtained in refs. [34,35].
The effects of soft wide-angle radiation (with momenta ∼ M 2 X /E cm ) are encoded in the soft function S ab . It depends on the color representation of the three hard partons and the total momentum component inn µ J direction of the soft partons. We also indicated a dependence on n ij ≡ n i ·n j with i, j = 1, 2, J, where the n µ i denotes the lightlike proton and jet directions in Minkowski space (n 2 i = 0), respectively. This dependence is however such that it exactly cancels the dependence on the jet energy E J = [(ŝ +t)(ŝ +û)n 12 /(2n 1J n 2Jŝ )] 1/2 in eq. (2.4) as can be made manifest by rescaling ω [18]. The soft function is the same for any Y and was calculated at two loops in ref. [44]. In the present paper we compute it at three loops, i.e. O(α 3 s ).
For later reference we also quote the Laplace transform of eq. (2.4) d 2σ where the convolutions turned into simple products off i ,J c , ands ab , which denote the PDFs, jet, and soft functions in Laplace space [18], respectively, and In section 5.1, we derive the three-loop anomalous dimension of the soft function in Laplace space via the renormalization group (RG) invariance of eq. (2.5) and check the consistency with our explicit calculation. In SCET the PDFs, the jet function, and the soft function are defined as operator matrix elements. For illustration we quote here the relevant (unrenormalized) expressions for the quark and anti-quark PDFs of a proton (p n i ) with momentum P − n µ i /2 employing the SCET 3 The factorization was performed here on the partonic level for small mX . The factorized partonic cross section is then convoluted with the (threshold) PDFs to obtain eq. (2.4). 4 We define n 2 i = 0 andni ·ni = 2 for i = 1, 2, J.
label formalism [23,45]: The composite quark field operator χ n i is gauge-invariant w.r.t. collinear gauge transformations and includes a collinear Wilson line in its definition, see e.g. ref. [45] for details. The label momentum operator P n i returns the total large light-cone (minus) momentum component of the n i -collinear fields it acts on. A similar expression holds for the gluon PDF [23,45]. The jet function and PDFs have trivial color structure, i.e. they do not have open color indices. In contrast, the hard and soft functions can be regarded as color tensors with one index for each parton in the hard scattering amplitude and one for each parton in the complex conjugate amplitude. Corresponding color indices of soft and hard functions are pairwise contracted in the factorized cross section. (In eqs. (2.4) and (2.5) color indices have been suppressed for the sake of compactness.) A convenient way to deal with the colors structure of scattering processes with more than two partons is the color-space formalism [46,47]. Applying it to our factorized cross section the soft function represents a square matrix in color space acting to the left and right on color vectors corresponding to the hard scattering amplitude and its complex conjugate. In general we write for the hard function and the corresponding hard amplitude 5 where the a i (and b i ) are color indices of the parton legs and a 1 a 2 . . . k are (not necessarily orthonormal) vectors spanning the color space. In the special case of three hard partons the color space is only one-dimensional and we define 6 for the qg → q channel, if a 2 a J a 1 for the gg → g channel, where t a bc and if bac are the SU (N c ) generators in the fundamental and the adjoint representation, respectively. 7 We thus have (for each production channel) where S a 1 a 2 ... represents the soft function as a tensor with color indices and S represents the soft function as an operator acting on color space vectors, which is proportional to the unit operator (S = S 1) in our one-dimensional color space. Using this notation the soft function operator is defined by [18] The hard amplitude C equals the QCD scattering amplitude with hard parton legs, where in the course of the SCET matching procedure, the IR divergences are subtracted. 6 Here and in the following we often neglect theqg →q channel, because it is trivially related to the qg → q channel by charge conjugation. In QCD therefore all physical results for the two channels are the same. 7 The totally-symmetric tensor d abc = tr[{t a , t b }t c ]/TF is ruled out by charge conjugation symmetry of QCD as another possible color structure for the gg → g amplitude, see e.g. ref. [48].
where T (T) indicates (anti-)time ordering,p µ is the soft momentum operator, and are the soft Wilson lines for an incoming and outgoing parton i, respectively. The symbol P (P)denotes (anti-)path ordering of the color charge operators T c i and the associated (ultra)soft SCET gluon field operators A c µ (x). The action of the T c i on color space vectors is defined by if parton i is a outgoing quark or incoming antiquark, if parton i is an incoming quark or outgoing antiquark, where C R is the quadratic Casimir of the color representation R of parton i (R = F for (anti)quarks and R = A for gluons) and

Dispersive method
The soft function defined in eq. (2.12) represents an integral of a squared matrix element of one outgoing and two incoming Wilson lines. This can be made explicit by inserting a complete set of (soft final) states (X), Rather than evaluating phase space integrals we would like to compute the soft function in terms of forward-scattering-type loop diagrams. A similar 'dispersive' method was used for the calculation of the soft function for heavy-to-light decays near the kinematic endpoint as detailed in ref. [39]. Following this approach, we will rewrite the δ-function in eq. (2.12) as a (momentum-space Wilson line) propagator-type expression using where iδ represents an infinitesimally small positive imaginary part. Our goal is to express the soft function in terms of a single time-ordered (forward-scattering) matrix element, which can be evaluated using ordinary (Wilson line) momentum space Feynman rules. As usual in momentum space perturbation theory the time-ordering is effectively implemented via the causal 'i0' prescriptions for the propagators in the corresponding Feynman graphs. Using eq. (3.2) means that we will introduce another independent imaginary infinitesimal in the computation and special care has to be taken in order to avoid unphysical imaginary parts due to interference of iδ and i0 in loop graphs giving rise to physical thresholds. 8 Concretely, we implement the δ-function in eq. (2.12) as follows At first sight taking the real part in eq. (3.3) may seem redundant. It will however become important due to our treatment of iδ and i0 in the following derivation. We will illustrate the necessity of taking the real part in our dispersive approach for the soft function explicitly at two loops at the end of this section. Note that the operator matrix element in the last line of eq. (3.3) alone is ambiguous. We therefore conveniently define Here we inserted +i0 in the denominator in order to obtain a well-defined analytical expression for real ω and to make σ(ω) directly calculable via Feynman diagrams as we will demonstrate below. With eq. (3.4) we can write where the i0 → 0 limit inside σ(ω) is understood to be taken before the discontinuity. Note that we could just as well define σ(ω) with i0 → −i0, which would lead to Disc ω σ(ω) → [Disc ω σ(ω)] * . Taking the real part in eq. (3.5) removes this ambiguity and is therefore necessary to ensure the consistency with eq. (3.3).
Let us now consider the evaluation of σ(ω). First, we notice that J,out are already (anti)time-ordered by default as a consequence of the (anti)path-ordering in eq. (2.14), and because the fields at positive times 8 By 'physical threshold' (not to be confused with the large transverse momentum 'threshold' of the EW boson) we refer to a branch cut along physical values of some kinematic invariant the soft function depends on besides ω. The corresponding amplitude for the heavy-to-light soft function discussed in ref. [39] is free of such thresholds. The subtleties related to the difference of iδ and i0 are therefore absent in that case and eq. (3.2) can directly be applied with iδ i0. from Y ( †) J,out are already to the left (right) of the fields at negative times in the T (T) product. Using we can thus write For an evaluation of σ(ω) via forward scattering Feynman diagrams all field operators in eq. (3.8) must be time-ordered. In the present form they are already time-ordered except for the operators in the T product. We can bring them into time order by the procedure sketched in the following, see also ref. [39] for a similar argument. We first recombine σ(ω) with the n 1 -and n 2 -collinear matrix elements, i.e. the PDFs. We can now undo the BPS field redefinition [22] of the collinear fields, which made the decoupling of soft and collinear sectors at leading order in the SCET expansion manifest. While the following arguments hold in general, let us consider the qq → g channel for concreteness. The resulting soft-collinear matrix element then has the schematic form where we suppressed the Dirac structure as well as the delta functions fixing the label momentum of the collinear (χ) fields from eqs. (2.7) and (2.8) for brevity. As r > 0 the fields in eq. (3.9) are in fact time-ordered. We can make this manifest by inserting the time-ordering operator T and redo the BPS field redefinitions χ n 1 → Y 1,in χ n 1 and χ n 2 → Y 2,in χ n 2 . 9 Note, however, that the consistent interpretation of eq. (3.9) as forward matrix element now requires that p n 1 p n 2 | = out p n 1 p n 2 |. This entails that the field redefinition of the interpolating (anti)quark fields at time +∞ generating part of the outgoing proton-proton state gives rise to an additional factor Y 1∞ Y 2∞ [49] with to the right of the Y † i,in . We emphasize that the color charge operator T i in Y i,out is still that of an incoming parton i = 1, 2 w.r.t. eq. (2.16).
We can now factorize out the PDFs again, which are matrix elements of local operator products in SCET and as such unaffected by the time-ordering, and finally end up with (3.11) 9 To properly keep track of the color indices, see e.g. the detailed derivation of the factorization theorem based on the BPS field redefinition in ref. [18]. This expression for σ(ω) admits a straightforward evaluation in terms of forward-scatteringtype loop diagrams using QCD Wilson-line Feynman rules in momentum space.
One of the relevant Feynman graphs at two-loop order is shown in fig. 2. Applying the Feynman rules and dimensional regularization (d = 4 − 2 ) we have The left-right and up-down mirror graphs yield the exact same result. The color factor C S depends on the channel and is given in eq. (5.3). To obtain the corresponding contribution to the soft function according to eq. (3.5) we first take the discontinuity in ω (after ω + i0 → ω) using The resulting expression has an imaginary part from the branch cut in the (−n 12 −i0) −2 factor for (physical) n 12 > 0 (related to a physical threshold in the full QCD one-gluon emission amplitude for the process ab → cY ). In a calculation based on eq. (3.1) as performed e.g. in ref. [44] such imaginary parts also appear in the soft amplitude, but cancel in the product with the complex conjugated amplitude. Equivalently, using two-loop diagrams with a unitarity cut (through a single gluon line) the imaginary part cancels between the diagram corresponding to fig. 2 and its left-right mirror diagram. In our approach the imaginary part is removed by explicitly by taking the real part in eq. (3.5). At two loops the diagram in fig. 2 (together with its mirror graphs) is the only one with an imaginary part after taking the discontinuity. At three loops we found quite a number of such diagrams.

Calculation
In this section, we present details of our three-loop calculation of the soft function based on eqs. (3.5) and (3.11). The calculation is performed in general covariant gauge with gauge parameter ξ, where ξ = 0 corresponds to Feynman gauge. Ultraviolet (UV) and infrared (IR) divergences are regulated using dimensional regularization (d = 4 − 2 ). We generate the relevant three-loop Feynman diagrams with qgraf [50]. Some examples of non-trivial three-loop graphs are shown in fig. 3. Based on dimensional counting and investigating the behavior of the corresponding loop integrals under rescaling we can directly identify already at this stage classes of diagrams that vanish as scaleless integrals, as illustrated in fig. 4. We are then left with 1290 non-trivial diagrams to be computed. These diagrams are further processed by a private Mathematica code [51] which first assigns the momentum space Feynman rules and performs simplifications of the Dirac, Lorentz and color structure. After that the diagrams are expressed as linear combinations of scalar Feynman integrals. Exploiting the n 1 ↔ n 2 interchange symmetry these can be mapped onto twenty integral families corresponding to different sets of fifteen linearly independent linear (Wilson line) and quadratic propagators. The mapping of scalar integrals to specific families requires partial-fraction decomposition of linear propagators followed by suitable shifts of the loop momenta. In order to automate the extensive partial fractioning we employed the algorithm described in ref. [52].
Next, we use the public program FIRE5 [53] to perform the integration-by-parts (IBP) reduction of the integrals in each of the twenty integral families to a set of master integrals (MIs). In three of the families we find an additional partial-fraction identity among the MIs returned by FIRE5, reducing their number by one, respectively. We then identify redundant MIs across the families, i.e. subsets of MIs with the same integrand up to loop momentum shifts despite being members of different families. Finally, the sum of all three-loop diagrams contributing to σ(ω) can be expressed as a linear combination (with real d-dependent coefficients) of 73 linearly independent MIs belonging to eleven different integral families. In this expression the gauge parameter ξ manifestly cancels out, which represents a strong check of our setup. The integrals in each family can be written as where the denominators D 1,··· ,6 correspond to quadratic propagators, and D 7,8,9 , D 10,11,12 , and D 13,14,15 correspond to propagators of Wilson lines along the n 1 , n 2 and n J directions, respectively. To give an example, one of our integral families is defined by where D i → D i − i0 in eq. (4.2) is understood. Our definitions of propagator denominators D i for the remaining ten integral families are given in app. C. For all eleven families, the result of the generic integral in eq. (4.2) is constrained by the scaling properties of its integrand w.r.t. eq. (4.1) to be of the form 10 where The dependence on the external kinematics is thus completely fixed and factored out. What is left to be computed is the dimensionless, in general complex, function I( a, b, c, e, ). In fact, because of eq. (3.5), we can safely set ω = −1 and only evaluate the real parts of the relevant MIs as an expansion in .
The ω dependence of the MIs can finally be restored according to eq. (4.4).
For the analytic computation of the MIs we follow the approach of ref. [34], which was inspired by refs. [54][55][56]. See also refs. [39,57] for recent applications and some more details of this method. The basic strategy is to express each MI in terms of known integrals (with less propagators) and integrals that are quasi-finite in d = 4 − 2 or higher, in practice d = 6 − 2 or d = 8 − 2 , dimensions using dimensional recurrence relations [58][59][60] and IBP reduction. The integrands of quasi-finite integrals in the Feynman parameter representation are by definition free of endpoint singularities. For each MI we determine a suitable set of related quasi-finite integrals using the program Reduze2 [61]. We then expand their integrands and perform the integrations over the Feynman parameters order-by-order in using the Maple package HyperInt [62]. In order to deal with threshold singularities (related to the 'physical thresholds' discussed in sec. 3) HyperInt automatically performs a contour deformation by adding an infinitesimal imaginary part δ k iε with δ k = ±1 to the kth Feynman parameter. This procedure introduces some ambiguity (dependence on δ k ) in the imaginary part of the result, which is however irrelevant for our calculation, because only the unambiguous real part of the integral eventually contributes to our soft function. We discuss the evaluation of a three-loop integral with non-zero imaginary part as an example in more detail in app. D. To determine the required dimensional recurrence relations we use LiteRed [63,64].

Results
Evaluating all relevant three-loop Feynman diagrams as described in the previous section and summing up their contributions yields σ(ω) in eq. (3.11) at O(α 3 s ). According to eq. (3.5) and using eq. (3.13) we obtain the bare N 3 LO soft function 11 where the coefficients K X are given in app. B. The subscript η = qq, qg, gg indicates the partonic (production) channel. All dependence of the soft function on the kinematic variables n 12 , n 1J and n 2J is encoded inω which is invariant under the rescaling in eq. (4.1). This follows from the (re)scaling invariance of the Wilson lines in eq. (2.12). Following ref. [44] we have included a factor of √ 2 in the definition ofω for convenience.
The color factor C Sη depends on the partonic channel and can be expressed in terms of the quadratic Casimir invariants C F = (N 2 c −1)/(2N c ) and C A = N c of SU (N c ) with T F = 1/2 [44]: The color factor C η 3 originates from tripole color structures that first appear at three loops and is computed below. The number of light (massless) quark flavors is denoted by n f . For practical reasons our code performing the color algebra for the three-loop Feynman diagrams returns the color factors expressed in terms of N c . Due to non-Abelian exponentiation [66,67] and Casimir scaling of the color dipole contributions we are nevertheless able to uniquely reconstruct the color factors in eq. (5.1) in terms of Casimir invariants as explained in sec. 5.2.
We define α s ≡ α s (µ) to be the renormalized QCD coupling constant. We stress that Z α µ 2 α s (= α bare s ) and therefore the bare soft function is independent of the renormalization scale µ. For renormalized quantities we use the MS scheme throughout this work. The relevant terms of the strong coupling renormalization factor Z α are Following ref. [18] we conveniently perform the renormalization of the soft function in Laplace space. This is, because the Laplace transformation according tõ

Anomalous dimension
In Laplace space, the bare and renormalized soft function are related by the local renormalization factor Z sη ,s bare η (κ) = Z sη (µ)s η (L, µ) . with the soft anomalous dimension Γ Sη . The Laplace space RGEs for the PDFs, the hard, and the jet function take the same form. RG invariance of the near-threshold cross section in eq. (2.5) implies ({a, b, c} = {q,q, g}, {q, g, q}, or {q, g,q}, or {g, g, g}, η = ab) where 12 are the (x → 1 threshold) PDF and jet function anomalous dimensions in Laplace space, respectively. The variables Q 2 and τ i are given in eq. (2.6) and we conveniently define Γ q cusp = C F γ cusp and Γ g cusp = C A γ cusp for the universal cusp anomalous dimension. The anomalous dimension of the hard function is denoted by Γ Hη . For scattering amplitudes with three massless partons the corresponding operator in color space reads [69,70] where T i is the color charge operator for the i-th parton leg of the scattering amplitude defined in eq. (2.16) and s ij ≡ 2σ ij p i · p j + i0, where the sign factor σ ij = +1, if the momenta of i-th and j-th parton p i and p j are both incoming or outgoing, and σ ij = −1 otherwise. The sums are over all combinations of distinct parton indices. The non-cusp anomalous dimension γ i is associated with each external massless (anti)quark (γ i = γ q ) or gluon (γ i = γ g ). The last term in eq. (5.14) corresponds to a color tripole contribution and is non-zero starting from three loops. The color structure T ijkl is given by [69] T The three-loop contribution to the associated coefficient was computed in ref. [71] and reads a 1 a 2 a J T iijk a 1 a 2 a J a 1 a 2 a J a 1 a 2 a J (5.19) 12 In this paper we use the same convention for γ J i as refs. [18,68]. In order to switch to the convention of refs. [34,39] one has to multiply our γ J i by −1/2.
we obtain from eq. (5.11) where γ Sη is the non-cusp piece of the soft anomalous dimension, and we have for the different parton channels. In order to obtain the explicit expressions for the color factor C η 3 in terms of the quartic Casimir invariants R denotes the fully symmetric rank-four tensor of the SU (N c ) representation R, and N A = N 2 c −1, N F = N c , we employed the color code [72]. All anomalous dimensions on the right-hand sides of eq. (5.21) are known at least to three-loop order, see app. A. 13 Just like the cusp also the non-cusp anomalous dimension γ Sη obeys Casimir scaling up to three loops, i.e.
Knowing Γ Sη we can determine the soft function renormalization factor Z sη via Γ Sη = −d ln Z sη /d ln µ. Up to three loops we find [69,76] 14 25) and we have expanded the anomalous dimensions as (5.26) 13 Starting from four loops the cusp and non-cusp anomalous dimensions violate Casimir scaling [73][74][75]. 14 for brevity we suppress the subscript η of the soft function in following.
As expected, Z s absorbs all divergences of our explicit result for the bare Laplace space soft functions(κ). This represents a strong check of our calculation and at the same time confirms the universal infrared structure of QCD scattering amplitudes with three massless parton legs, as predicted by eq. (5.14), at three loops.

Renormalized results
As mentioned above our results for the bare soft function in the different parton channels are initially expressed in terms of N c . In order to rewrite the color factors in terms of (quadratic and quartic) Casimir invariants we proceed as follows. Non-Abelian exponentiation [66,67] restricts the three-loop contribution to our soft function for all parton channels to be a universal linear combination of the color factors First, we discuss how to determine the coefficients of the color factors for n f = 0. For the gg → g channel we have according to eqs. (5.3) and (5.21) C S = N c /2 and C 3 = −9N c . So, all color factors except for C 3 are proportional to N 3 c . We can therefore directly identify the color tripole contribution ∝ C 3 in our soft function result. Once this is fixed we take our result for the qq → g channel and subtract the tripole term with C gg 3 → C qq 3 . The remaining terms in S bare qq can now be uniquely expressed in terms of the dipole-type color factors in eq. (5.27) using the replacement rules Concerning the color factors involving n f , we first map the C S C F n f T F contribution to the n f T F N 0 c term in our result for S bare gg . After that we can rewrite the remaining n f dependent terms in S bare qq by replacing In this way S bare qq is completely expressed in terms of Casimir invariants and matches eq. (5.1). We successfully checked the universality of this expression against our explicit results in terms of N c for S bare qg and S bare gg by adjusting C S and C 3 according to eqs. (5.3) and (5.21). Solving the RGE in eq. (5.10) with the anomalous dimension in eq. (5.20), the renormalized soft function in Laplace space can be written as Here we have expanded the anomalous dimensions in eq. (5.21) as Our calculation determines the non-logarithmic terms in eq. (5.30) to be The expression for c S 3 is new and represents the main result of this work. For completeness we also give here the renormalized soft function in momentum space obtained from inverting eq. (5.6). The coefficients in the α s expansion take the form We have calculated the universal three-loop soft function contributing to factorized cross sections for EW boson production in the (threshold) limit, where the transverse momentum of the boson is close to its kinematically allowed maximum for a given (not too large) rapidity.
To the best of our knowledge this represents the first result of a soft function for a threeparton process at three loops. We have derived a novel expression for the soft function in terms of a forward-scattering-type matrix element of Wilson line operators. This allowed us to avoid any phase space integrations and to straight-forwardly take advantage of well-established multi-loop technology in our calculation. We present our results for the renormalized soft function both in Laplace space and momentum space in sec. 5.2. The three-loop non-logarithmic terms in eq. (5.32) represent the genuinely new information at this order. Our explicit three-loop calculation validates the relation between the anomalous dimensions of the PDFs and the corresponding hard, jet, and soft functions inferred from RG consistency of the near-threshold cross section. In this way we also confirm the nonzero (Casimir scaling violating) color tripole contribution to the three-loop soft and corresponding hard anomalous dimensions found in ref. [71].
The threshold approximation of the EW boson production cross section is strictly valid when the invariant mass of the hadronic radiation (including the proton remnants) is small compared to the transverse momentum of the boson. Phenomenologically, however, it usually performs well even far away from this limit, see e.g. refs. [25,26]. Together with the threeloop jet function [34,35] and the yet-unknown three-loop hard function our soft function result will allow to determine the N 3 LO corrections to the transverse-momentum spectra of γ, W ± , Z and Higgs bosons in the (far-)tail region. Furthermore, it is necessary to carry out the resummation of threshold logarithms to N 3 LL (or N 4 LL) accuracy. Once the QCD three-loop virtual corrections to the qq → Y g (or crossing-symmetric) and gg → Hg scattering amplitudes, i.e. the hard functions, become available, this may constitute significant progress in the precision phenomenology for Higgs, photon and EW gauge boson production in association with a jet at the LHC. Acknowledgments Z.L.L. is grateful to Thomas Becher, Matthias Neubert and Ding Yu Shao for valuable discussions. We are indebted to Robin Brüser for providing the codes to perform the diagram generation, topology mapping, and partial fractioning and for his help with the color code. This research has been supported by the Cluster of Excellence PRISMA + (project ID 39083149), funded by the German Research Foundation (DFG). The research of Z.L.L. is supported by the U.S. Department of Energy under Contract No. DE-AC52-06NA25396, the LANL/LDRD program and within the framework of the TMD Topical Collaboration. All graphs were drawn using JaxoDraw [77].

A Three-loop anomalous dimensions
For completeness we list here the explicit expressions for all anomalous dimensions in eq. (5.21) to three-loop order. The convention for the loop expansions is analogous to eqs. (5.26) and (5.31). The coefficients of the cusp anomalous dimension are [78,79] The coefficients of the soft non-cusp anomalous dimension γ S are [18] The hard anomalous dimensions γ q and γ g can be determined [69] from the divergent part of the on-shell quark and gluon form factors in QCD [80,81] and read up to three loops At three loops the anomalous dimensions of the quark and gluon jet functions were derived in refs. [18,68] based on RG consistency arguments and confirmed by the jet function calculations in refs. [34,35]. Their coefficients are and respectively. The anomalous dimensions describing the evolution of the quark and gluon PDFs near x → 1 can be determined from the QCD splitting functions [79,82] and were extracted in refs. [83,84] up to three-loop order:

B Bare data
Here we present our expressions for the coefficients of the different color structures in the bare soft function, eq. (5.1). We show the results as an expansion in = (4 − d)/2 to the order required for the calculation of the renormalized three-loop soft function:

C Definition of integral families
In section 4, we defined one integral family for illustration. All eleven families share the same propagator denominators D 1,··· ,6 . Here we present the definition of the other propagator denominators D 7,··· ,15 for the remaining ten families.