Double-Real-Virtual and Double-Virtual-Real Corrections to the Three-Loop Thrust Soft Function

We compute the ${\cal O}(\alpha_s^3)$ double-real-virtual (RRV) and double-virtual-real (VVR) soft contributions to the thrust/zero-jettiness event shape. The result clears up one of the most stubborn obstacles toward the complete ${\cal O}(\alpha_s^3)$ thrust soft function. The results presented here serve as the key input to realize the next-to-next-to-next-to-leading logarithmic prime (N${}^3$LL') % and even the next-to-next-to-next-to-next-to-leading logarithmic (N${}^4$LL) resummation of the thrust event shape. The obtained results also constitute the important ingredients of the $N$-jettiness-subtraction scheme at next-to-next-to-next-to-leading order (N${}^3$LO).


Introduction
The thrust distribution τ in e + e − annihilation process is among the fundamental tools to test the theory of strong interactions at high precision [1][2][3].It probes the global geometrical structure through an inclusive measurement of where t is known as the thrust axis that maximizes the last equation and ⃗ k i is the momentum of the particles produced in e + e − annihilation.Being theoretically clean and feasible to high perturbative order computations, the thrust distribution is particularly suitable for the precise determination of the strong coupling α s (M Z ) [4][5][6][7][8][9][10][11] and has been frequently measured at e + e − colliders with small experimental uncertainties.The precise extraction of the α s with unprecedented accuracy will continue to be the major scientific pillar at future e + e − facilities [12].
The thrust distribution τ probes the kinematic configurations between a perfectly isotropic multi-particle final state with τ close to the tail point 1/2, and the two collimated back-to-back energy flows (the di-jet limit) where τ → 0. When sufficiently away from the di-jet limit, the thrust distribution is well modeled by the fixed order calculation which has been known to O(α 3 s ) [4,13].Near the di-jet limit in which τ ≪ 1, resummation is required to secure the perturbative predictions.
When τ ≪ 1, the τ spectrum can be written in a factorized form [6,[14][15][16][17] within the soft collinear effective theory (SCET) [18][19][20][21][22], which gives dσ sing.dτ = σ 0 i=n,n,s where σ 0 contains all loop corrections to the matrix e + e − → q q and we assume that q is along the light-like direction n = (1, 0, 0, 1) while q along n = (1, 0, 0, −1).J is the inclusive jet function [23,24] and S is the soft function defined as the vacuum expectation of Wilson lines [6].The factorization is valid up to O(τ ) power suppressed contributions and with this approximation1 , the thrust measurement in Eq. (1.1) is simplified to The factorization theorem reproduces correctly the leading logarithmic singular contributions to the τ distribution such that where the logarithms L are given by and the C (n) k are some constants free of τ .The factorization theorem sets the foundation for precision prediction of the thrust distribution in the resummation region in which τ ≪ 1, and is the key ingredient to the precise determination of the α s .The logarithmic terms (i.e., C (n) k with k ⩽ 2n) can be calculated through the standard resummation techniques to all orders in α s for the jet and the soft functions, if the associated cusp Γ cusp and the non-cusp anomalous dimensions γ i are known, see Tab. 1.The cusp anomalous dimension is known analytically to 4loops [27,28] following earlier efforts [29][30][31][32][33][34][35][36][37][38].The thrust distribution is now known at the next-to-next-to-next-to-leading logarithmic (N 3 LL) accuracy [6,8,10].While the δ(τ ) terms (i.e., C (n) 2n+1 ) have to be determined through explicit fixed order calculations of the jet and the soft function order by order in α s and their accuracy also determines the logarithmic accuracy of resummation as shown in Tab. 1.

Anomalous Dimension
Table 1.The definitions for the logarithmic accuracy of the resummation calculation.
To push through the precision to the next level, the complete 3-loop calculation for the jet and the soft function have to be performed.The calculation of the N 3 LO quark jet function has been carried out in [39] years ago, while the soft function is only known to order O(α 2 s ) [40][41][42][43].Recently, as a first non-trivial step toward the three loops, the samehemisphere three-gluon real-real-real (RRR) contribution to the soft function is obtained analytically [44].
In this manuscript, we present the computation of the double-real-virtual (RRV) and the double-virtual-real (VVR) contributions to the thrust soft function.Since the 2-loop soft current is known in the literature for a while [45,46], it makes the double-virtual-real (VVR) correction to the thrust soft function easy to obtain, and we present the VVR results in the Appendix A.
The work is organized as follows: In Section 2, we introduce the convention we use in this paper.In Section 3, we review the method to evaluate the integral in Eq. (2.6).We demonstrate the methodology by reproducing the known double-real contribution to the O(α 2 s ) soft function in Section 3.3.We present the O(α 3 s ) double-real-virtual correction in Section 4 and we conclude in Section 5.

Theoretical set-ups
We focus on evaluating the RRV corrections, S RRV (τ ; ϵ), which is given by where we used the notation for the 2 body phase space integration.We have explicitly written out the loop integral d D l.The matrix element ω RRV is given by the interference between soft currents with 2 real emissions at the tree and loop levels, ω Representative Feynman diagrams for ω RRV are shown in fig. 1.In our calculation, the Feynman diagrams and amplitudes involving Wilson lines are generated by Qgraf [47].Additional eikonal propagators and vertexes are included as well as those for conventional QCD in the model.To generate the correct amplitudes, we discard the diagrams containing internal eikonal lines.Throughout this paper, a i0 + prescription is used for each propagator.That is, the denominator D i of a propagator should be understood as D i + i0 + .Thus, we write each propagator in M (L) † as − 1 −D i +i0 + .The i0 + prescription is also used for the eikonal propagator, which originated from an outgoing soft gluon emitted from an outgoing energetic parton.For brevity, the i0 + will be suppressed hereafter.We note that for the RRV contribution, the final result is independent of the pole position in the eikonal propagator [48,49] and therefore the result applies to the thrust in e + e − annihilation, as well as the jettiness distribution in Higgs/Drell-Yan [50,51] or DIS [52,53].
The phase space restrictions due to the thrust measure in Eq. (1.3) are given by where we introduced the notation We further introduce the dimensionless variables by re-scaling the momenta as τ i → τi τ , k i → ki τ and l → lτ , to factorize out the τ dependence from the rest of the integral.We get RRV ( l, k1 , k2 ; n, n) Θ(1; k1 , k2 ) . (2.6) In the rest of the work, we will suppress theˆnotation.And we will show in detail how the integration in Eq. (2.6) can be evaluated.
For future convenience, we will always multiply each Dirac δ-function and the Heaviside θ-function by a factor of 2π.We will associate each fold of loop integration with a factor of Throughout the work, we denote D = 4 − 2ϵ for the dimensional regularization in 4-dimension, while we use the notation d for the regularization in a more generic dimensional space-time such as d = 6 − 2ϵ.We relate each fold of integration with respect to a scalar, such as dτ , with a factor of 1 2π .That is, instead of calculating S (L) , we calculate RRV/RR (1, ϵ).

Linear reduction of integrals
In this part, we briefly review the methodology developed in [54][55][56] to perform the reduction of the phase-space integrals with kinematic cuts and the conventional Feynman loop integrals in a unified way.An alternative approach was developed in [57].In real calculations, the kinematic cuts are implemented by inserting the Dirac δ-function or the Heaviside θ-functions at the integrand level in the phase space integrals.One example of such kinematic cuts can be seen in Eq. (2.1) in the previous section.The unified reduction of the phase space and the Feynman integral roots in the well-known fact that either δ-function or the θ-function has its own integral representation similar to that of a propagator.By virtue of this similarity, integrals with δ-functions and θ-functions can be reduced and calculated by generating linear relations and solving differential equations in the parametric representation.

Parametric representation
Here we discuss the parametric representation of the Feynman integrals.To unify the treatment of the phase space integrals and the standard Feynman loop integrals, we note that the θ-function and the δ-function belong to the integral class whose definition is given by Obviously, we have for λ = 0, λ = −1 and λ = −2, respectively.First, we consider the scalar integrals with the form where L is the number of loops, l i are the loop momenta, D i are the denominators of propagators, and we denote the integration measure as We start with the all nonnegative λ i case.As was shown in Refs.[54,56], the scalar integral can be parametrized by with s g = det(η µν ) the determinant of the d-dimensional space-time metric, and is a homogeneous function of x of degree −n − 1 and is given by where the polynomial F is related to the Symanzik polynomials U and F through F(x) ≡ F (x) + U (x)x n+1 , with U (x) ≡ det A, and Here A, B, and C are linear in x and are determined from The parametrization of the scalar integral J(λ 1 , . . ., λ i , . . ., λ n ) in Eq. (3.4) can be generalized to the cases with negative values of λ i for i > m, in which the parametric integral I is defined through j=m+1,j̸ =i Γ(λ j + 1) Apparently A tensor integral can also be parametrized in terms of I(λ 0 , . . ., λ n ) recursively in a similar way.That is where with Here p µ i is an auxiliary vector introduced through the identity [54] and D 0 and xi are the index-shifting operators will be given in Eq. (3.10) and Eq.(3.11).

Linear relations
The parametric integral satisfies the linear relations ) We introduce the index-shifting operators R i , D i , and A i , with i = 0, 1, . . ., n, such that And we formally define operators D n+1 and R n+1 , such that D n+1 I = I, and ).The product of two operators X and Y are defined by (XY )I ≡ X(Y I).We further introduce the operators xi , ẑi and âi such that Obviously we have ân+1 = −(L + 1)A 0 − n i=1 (â i + 1).For i = 1, 2, . . ., n, we have the following commutation relations: With the operators xi , ẑi , and âi , it is easy to write the IBP identity in Eq. (3.9) in the following compact form One should keep in mind that these operator equations are valid only when they are applied to nontrivial parametric integrals.

Integral reduction
Similarly to the traditional integration-by-parts (IBP) technique [58,59], parametric integrals can be reduced by solving linear relations generated by Eq. (3.13).Based on Eq. (3.13), two methods to reduce the parametric integrals were proposed in Ref. [55].For conventional Feynman integrals, Method II therein is more efficient by virtue of the non-negativity of the indices.However, for integrals with cuts, we have to be able to handle the negative indices, because a delta function is of an index −1.So in this case Method I is more advantageous.
Here we give a brief introduction to this method.Let a ij and b ij be the solutions of where A and B are defined in the paragraph below Eq. (3.5).In general, solutions of these two equations could be linearly dependent.We denote the linearly dependent part of the solutions by c ij .By default, we assume that these solutions are excluded from a ij and b ij and therefore a ij and b ij are linearly independent.For brevity, we denote , where Q µ u are the linearly independent external momenta.We assume that the Gram determinant Q u • Q v is invertible.By introducing certain auxiliary parameters (or, equivalently, by introducing some auxiliary propagators), we can always make the matrices ∂B j ∂b i and ∂B ju ∂b i invertible.That is, there are matrices α and β such that Then we have the following identities where E is the number of external momenta, and where g uv is the inverse of the Gram matrix Comparing with Eqs.(3.13), Eqs.(3.17) are of lower degrees in x.And they are free of D 0 .Thus they will not shift the space-time dimension.In practical calculations, we will use Eqs.(3.17) to generate linear equations between the parametric integrals.The generated equations are then solved by using the package Kira [60][61][62][63].By using the operators Bi , Eq. (3.6) can be traded by where pi are vectors such that pi • Q j = 0.The B are free of D 0 , and thus will not shift the spacetime dimension.Due to the definition of β ij,k , the B commute with A's.
Thus by applying operators P µ i , tensor integrals are parametrized by integrals of the form f ( A)I(− d 2 , . ..),where f ( A) is a sum of chains of A ij .Chains of A ij can be traded by sums of chains of Āij by solving (3.20) The right-hand side of this equation is of degree n in A ij , while the left-hand side is of degree n − 1.Thus, by solving these identities, we can reduce the degrees of A ij recursively.

Calculation of master integrals
After carrying out the linear reduction, we express the amplitudes in terms of linear combinations of master integrals (MIs), which can be evaluated by using the differential-equation method [64][65][66].As was shown in ref. [55], differential equations can be constructed in the parametric representation.We let y be a Lorentz scalar, and assume that A ij are free of y.
Then we have As is mentioned in section 1, the τ dependence of soft functions can be factored out.Thus the differentiation of MIs with respect to τ is trivial.To get a nontrivial differentialequations system (DES), we need to introduce an auxiliary scale.This can be done by inserting into the MIs a trivial integral of the form dyδ(K − y), where K is a linear combination of Lorentz scalars of the loop momenta.Thus MI J i becomes where the integration measure [dl i ] is defined below Eq. (3.3).The region of integration for y depends on the choice of K.The choice of K depends on the integrals to be evaluated.For instance, for integrals free of theta functions (all the theta functions are reduced to delta functions), we insert ; etc..The obtained integrals J i (y) can be calculated by using the standard differentialequation technique.Different from some papers (see e.g.ref. [67]), which regularize the boundary conditions of the differential equations in 6 dimensions, we calculate J i (y) in 6 dimensions from the very beginning.That is d = D + 2 = 6 − 2ϵ.The dimensional recurrence [68] can easily be carried out in the parametric representation.To convert the DES to the canonical form [66], we need to rationalize it first.In the case K = k 1 • k 2 − y, we do the transformation y = 1 2 (1 − x) 2 .In the case K = k , the rationalization is unnecessary, and we do the transformation y = 1 x .The rationalized DES can be converted to the canonical form by using the package epsilon [69], which implements Lee's algorithm [70,71].After converting the differential-equation system to the canonical form, solving it is straightforward.In the case where K = k 1 •k 2 −y, the singularities of the DES are {0, 1, 2}.In the case where K = k , the singularities of the DES are {-1, 0, 1}.In either case, the solutions of the DES can be expressed in terms of the harmonic polylogarithms [72], of which some can be converted to normal polylogarithms by using the package HPL [73].
It remains to determine the boundary conditions.We chose the boundary to be at x = 0.The boundary conditions are determined by matching the asymptotic solutions of the DES to the asymptotic expansions [74] of the MIs.Most of the asymptotically expanded integrals can easily be calculated in d dimensions.For those that are not easy to calculate, we determine them through the regularization condition: all 6-dimensional integrals free of double propagators are free of infrared divergences.
Integrals J i can be obtained from J ′ i (y) by integrating out y.For integrals with K = k 1 • k 2 − y, the integration is straightforward.For integrals with K = k . The integration can be carried out by subtracting the singular part of J(y) (denoted by Sing{J(y)}) and integrating the singular part in d dimensions.That is, Sing{J ′ (y)}. (3.23) Here we use Ser to denote the series expansion of a function with respect to ϵ.

An example
To illustrate how to carry out the strategy described in the previous subsections in practice, we present a detailed calculation for the two-gluon-emission diagram shown in Fig. 2. The amplitude for this diagram reads where the colour Casimir constants 2N C , with N C the number of colors of a quark.Here we use a subscript a to distinguish it from the full two-gluon emission contribution Ŝ(2) RR,gluon .For this amplitude, no tensor reduction is needed.Generally, we use eqs.(3.19, 3.20) to parametrize tensor integrals.The resulted parametric integrals are reduced by solving linear relations generated by using Eq.(3.17).For example, we consider the following integral family. .
The operators appearing in Eqs.(3.19, 3.20) are It is easy to see that the tensor reduction is carried out by using eqs.(3.19, 3.20) is consistent with the traditional tensor reduction method.For example, for an integral with a numerator k And an integral with a numerator k 2 1 is parametrized by (3.28) The IBP identities for this integral family are given by After carrying out the IBP reduction we further express the resulted amplitude in terms of 6-dimensional integrals through dimensional recurrence: The reductions for all the other families of integrals are similar.We finally find Here the MIs J i are give by 2 For the readers' convenience, we have expressed the parametric integrals in terms of the normal momentum-space integrals.These MIs J i can be divided into two classes.The measure of one class includes a delta function δ(k . And the that of another class includes a delta function δ(k + 1 + k − 2 − 1). 3 Here we consider the first class.That is, integrals J 1−4 .
The first integral J 1 is free of theta functions.We insert to it a delta function δ(k 1 • k 2 − y), and replace y by y = 1 2 (1 − x) 2 .As in the convention used in the last subsection, we use J ′ i (y) to denote an integral resulting from inserting a delta function into J i .The differential equation for J ′ 1 is simple: which can be solved easily.The boundary condition can also be readily determined.
For the integrals J 2,3,4 , we insert a delta function δ(k − 1 +k − 2 −y).The resulting integrals 2 According to the i0 + prescription of the denominator of a Feynman propagator, an integral with a denominator, say, k − 1 + k − 2 , will be treated as independent of the integral with a denominator 3 For this calculation, all the integrals with delta functions δ(k can be identified with the integrals with delta functions δ(k , by virtue of the symmetries under permutations of indices in the parametric representation. J ′ 2 (y) and J ′ 4 (y) can further be reduced.We have To get a closed DES, we need to add one more integral The DES for the integrals which can easily be transformed to the canonical form by using epsilon [69].
To determine the boundary conditions, we change the variable to x = 1 y .The asymptotic solutions of the DES are of the form Among all the C 1,i,j , only C 1,10,−1 is independent, which can be determined by calculating the asymptotic expansions of J ′ 10 .And we note that all the C 2,i,j and C 3,i,j vanish.To see this, we rescale all the loop momenta by a factor of x.For example, Due to the constraints of the delta functions, it is easy to see that there is only one region, in which x, and k − 1 ∼ 1.Thus J ′ 10 (x) ∼ x ϵ−1 , so only C 1,i,j survive.Once the boundary conditions are determined, solving the DES is straightforward.
It remains to calculate J i by integrating J ′ i (y) with respect to y (or x).Obviously, the integration is divergent when y → ∞ (or x → 0).To carry out the integration, we need to subtract J ′ i (y) by their singular parts and integrate the singular parts in d dimensions.The singular parts are nothing but the asymptotic solutions we obtained before.
Finally, we solved all the MIs, J 1 (2π) For completeness, we present the complete double-real contributions (up to O(ϵ 2 )) in Appendix B. The results are in full agreement with the known calculations [40][41][42][43], which serves as a strong validation of our computational framework.

Results
Now we present the result for RRV.We start with the 2 gluon emission contribution.In the virtual loop, all parton propagators are included.We decompose the two-gluon-virtual contribution according to the color structures by4

Ŝ(3)
RRV,gg =(4π where N f is the number of the quark flavors, and C A and C F are color structures defined at the beginning of section 3. and Similarly, for the two-ghost-virtual contribution and the fermionic contribution, we organize the contribution according to the color decomposition as

Ŝ(3)
RRV,q q =2(4π) 6 cos(πϵ)α A is validated by the strongly ordered limit calculation.Also we observed the cancellation of the C F C 2 A ϵ −5 -pole among RRR (c.f.eqs.(35)(36)(37)(38) in ref. [44]), RRV and VVR (see Eq. (A.13)) as required by the maximal non-abelian exponentiation theorem.We note that in the direct calculation, the highest pole is obtained as a consequence of the non-trivial cancellation between higher spurious poles, therefore the correctness of the leading pole serves as a non-trivial check of our approach.Another strong check is on the C 2 F C A term, where the direct bruteforce computation completely agrees with the non-abelian exponentiation theorem which guarantees that the C 2 F C A term can be obtained by convolution of the O(α 2 s ) RV correction and the O(α s ) soft function in Eq. (A.18) and Eq.(A.17), respectively.Further checks have been performed for the class of contributions, e.g.vacuum polarization of the gluon propagator, in which we can work out the virtual loop integration first and then evaluate the double real emissions analytically.In those checks, again we find the complete agreement with the computational setups using tensor reduction and MIs in this manuscript.continuation of the VVR to these processes is straightforward.Furthermore, the results presented here are the indispensable components to fulfill the N 3 LO calculation using the N -jettiness subtraction scheme [75][76][77].
The calculation demonstrates the feasibility of the recently proposed computational framework for a unified treatment of the loop and the phase space integration with Heaviside functions [54][55][56].We note that beyond the thrust, the methodology used in this work can be used to other phase-space integrals that involve experimental cuts through the Heaviside functions, and we particularly expect its applications to future jet and substructure precision calculations [78].
Note: In the earlier versions of this manuscript, we did not present the result for the fermionic contribution, which is added in this version.We note that a similar calculation was recently reported in ref. [79].Perfect agreement is found between these two independent calculations.
where ψ is the digamma function and p F q are the generalized Hypergeometric functions.
The integration in Eq. (A.1) can be evaluated straightforwardly which gives the VVR correction to the thrust soft function with renormalized strong coupling constant α s,r

S
(3) where where the first term comes from the interference between the 2-loop and the tree-level soft current and the second contribution from the square of the 1-loop soft current.
Here we have performed the renormalization of the strong coupling constant α s , which is given by α s µ 2ϵ → α s,r (4π) −ϵ e γ E ϵ µ 2ϵ 1 − α s,r 2π

Figure 1 .
Figure 1.Representative Feynman diagrams for the RRV contribution at O(α 3 s ).Dotted lines represent the ghost propagators while the vertical dashed lines denote the cuts over the real emissions.

Figure 2 .
Figure 2. A diagram for the two-gluon emission.