Next-to-leading power endpoint factorization and resummation for off-diagonal"gluon"thrust

The lack of convergence of the convolution integrals appearing in next-to-leading-power (NLP) factorization theorems prevents the applications of existing methods to resum power-suppressed large logarithmic corrections in collider physics. We consider thrust distribution in the two-jet region for the flavour-nonsinglet off-diagonal contribution, where a gluon-initiated jet recoils against a quark-antiquark pair, which is power-suppressed. With the help of operatorial endpoint factorization conditions, we obtain a factorization formula, where the individual terms are free from endpoint divergences in convolutions and can be expressed in terms of renormalized hard, soft and collinear functions in four dimensions. This allows us to perform the first resummation of the endpoint-divergent SCET$_{\rm I}$ observables at the leading logarithmic accuracy using exclusively renormalization-group methods. The presented approach relies on universal properties of the soft and collinear limits and may serve as a paradigm for the systematic NLP resummation for other $1\to 2$ and $2\to 1$ collider physics processes.


Introduction
Hadronic event shape variables in the two-jet region have played a key role for the development of resummation techniques in QCD [1,2]. 1 Later, the utility of soft-collinear effective theory (SCET) in extending the resummation order beyond what was achieved with diagrammatic methods has first been demonstrated for the thrust variable T [7]. In the present work we develop the factorization of thrust for the particular contribution to the two-jet region, when a gluon recoils (at leading order) against a quark-antiquark pair ("gluon thrust"). The motivation for this derives from the fact that this process is of next-to-leading power in 1 − T 1 in the two-jet region, and thus not covered by the existing factorization theorem. It represents the hadronic e + e − annihilation analogue of the off-diagonal qg parton channels in deep-inelastic scattering (DIS) at large Bjorken-x, or in Drell-Yan (DY) production near threshold. These and the power corrections to the corresponding diagonal channels have recently received much attention [8][9][10][11][12][13][14][15][16][17][18][19][20] in an effort to extend the classical leading-power (LP) factorization theorems for these observables to next-to-leading power (NLP) in the kinematic regime, where large logarithms spoil the weak-coupling expansion. We focus on gluon thrust here, since it presents the same difficulties as the forementioned off-diagonal processes, but does not require the factorization of parton distribution functions.
The thrust of a hadronic event in e + e − annihilation with center-of-mass energy Q is defined as where i runs over all hadrons (partons) in the final state. The plane orthogonal to the thrust axis n divides space into a left and right hemisphere. The total invariant mass (squared) of the hadrons in these hemispheres will be denoted by M 2 L and M 2 R , respectively. As the particles cluster into a pair of back-to-back jets, and perturbation theory in the strong coupling α s (Q) breaks down due to large logarithms ln τ at every order. This has been extensively studied for the leading quark-antiquark two-jet process. All-order resummations are essential to provide reliable calculations of thrust and jet masses in the two-jet regions. In this work we consider instead the phase-space region where the direction of gluon jet is defined to be the "collinear direction" and the direction of the recoiling qq-jet the "anti-collinear" one. This process starts at order α s only (see Fig. 1). More importantly, as τ → 0, it does not have the leading-power [ln τ /τ ] + γ * (Q) q(p a ) Figure 1: Leading O(α s ) QCD diagrams contributing to gluon thrust.
behaviour from soft-gluon emission. Instead the leading term is α s ln τ , and the entire process is of next-to-leading power in τ , 2 similar to the off-diagonal DIS and DY processes related by crossing. "Gluon thrust" and its unconventional all-order logarithmic structure was first studied in [14] and an expression for the leading double logarithms was written down, which was subsequently derived from d-dimensional consistency relations [18]. The result of [14] exhibits an unconventional "quark" Sudakov form factor, related to an enhancement when either the quark or anti-quark in the qq-jet becomes soft. In this case the octet colour charge of the gluon jet is balanced by anti-collinear particles in a colour-triplet or anti-triplet. The colour mismatch of the back-to-back energetic particles causes a new type of double logarithms proportional to the difference C A − C F of the colour charges at the leading double-logarithmic order. Despite these insights, and sketches of the form of a factorization theorem within SCET [11,14,18], it is not yet known how to proceed beyond the double-logarithmic order. The issue is an endpoint-divergence in the convolution integrals that link the hard, (anti-) collinear and soft functions in the factorization theorem. In the present work we build on refactorization ideas developed for DIS [18], and Higgs decay to two photons through bottom-quark loops [21,22] to derive a factorization theorem suitable for systematic resummation. To this end we combine standard SCET factorization with endpoint factorization.
The outline of this paper is as follows: In Sec. 2 we provide a heuristic discussion of the factorization formula, which motivates the proposed endpoint subtraction. In Sec. 3 we derive the "bare" factorization theorem using standard SCET methods. For endpoint factorization to work, the hard, collinear and soft functions must satisfy certain asymptotic conditions, which are derived in Sec. 4, and subsequently employed to rearrange the factorization formula such that the convolution integrals are finite. The renormalization-group equations required for the various functions are given in Sec. 5 and solved with leading-logarithmic (LL) accuracy. Finally, in Sec. 6, we obtain the resummed gluon thrust distribution and display the numerical effect of resummation.
Throughout the paper, we employ dimensional regularization with d = 4 − 2 spacetime dimensions, and the MS subtraction scheme whenever renormalization is required.

Heuristic discussion
At leading order (LO) a virtual photon couples to a quark-antiquark pair producing two back-to-back jets in the center-of-mass frame. In this work, we consider instead the different two-jet situation, where one of the jets is initiated by a gluon. This process is not only suppressed by one power of α s (see Fig. 1) relative to the LO process, but also by one power of τ .
We define the hemisphere of the gluon jet to be in the "collinear" direction, the other jet is then in the "anti-collinear" direction. There are two possibilities to generate the gluon jet at O(α s ): I Quark and anti-quark have both large anti-collinear momentum and form a single jet, which recoils against the gluon.
II The quark or anti-quark is anti-collinear and balances the gluon momentum, while the other of the two is soft.
Both configurations are suppressed by a factor of τ . Their further evolution is determined by the standard leading-power collinear splittings and soft emissions. The above separation introduces an ambiguity from the precise meaning of "soft" and "anti-collinear". We can start with situation I, and decrease the large anti-collinear momentum component of the quark or anti-quark. At some point, it becomes soft (and anti-collinear) and should be counted as part of contribution II. In perturbative computations in powers of α s , this ambiguity is not a problem as long as the jet definition is infrared-safe. If, however, one wants to perform resummation, the effective field theory approach requires the separation of soft and collinear modes to split the large logarithms into single-scale quantities. This conflict between physical and mathematical mode separation leads to a so-called endpoint divergence in the convolution integrals appearing in the factorization theorems. Let us analyze these two possibilities from the SCET point of view by looking first at the Feynman diagrams shown in Fig. 1. We begin with I (later called "B-type" contribution, following the SCET notation for the corresponding hard vertex). The intermediate propagator in the diagram on the left carries momentum q a = p g + p a . Since p g is assumed to be collinear and p a is anti-collinear, q a is a hard momentum with virtuality q 2 a ∼ Q 2 . The q a -internal line is integrated out, and the primary hard SCET vertex produces a qqg state directly, see the first diagram in Fig. 2. Since there are two partons (qq) in the anti-collinear direction, only their total momentum is fixed and the amplitude depends on the fraction of anti-collinear momentum carried by each parton. When one of these fractions tends to zero, the parton becomes effectively soft and one moves to possibility II (called "A-type" later). The q a -internal line (or q b , if the antiquark becomes soft) is no longer hard. Hence the primary vertex is the standard γ * → qq process. The energetic quark or anti-quark subsequently transfers its entire momentum to the gluon, while becoming soft. This process is represented by a power-suppressed SCET Lagrangian term L  The fact that in some intermediate region the same physical process can be represented by the limits of two different mathematical expressions is the key to the idea of endpoint factorization. To anticipate what will be discussed in technical terms further below, we state the factorization theorem for the two-hemisphere invariant mass distribution of gluon thrust in schematic form in Laplace space: Here C are hard matching coefficients, J are jet functions and S are soft functions.
The precise definitions will be given in subsequent sections. We dropped all function arguments except for the convolution variables, which have divergent integrals. The soft-collinear ω, ω convolution integrals are logarithmically divergent for ω, ω → ∞, while the hard-anti-collinear integrals over r, r diverge logarithmically when r, r → 0, 1. However, both integrals are well-defined in dimensional regularization as long as the limit to four space-time dimensions is not taken.
In the intermediate region, the soft quark in the soft function S NLP that appears in possibility II has large soft momentum ω and can equally well be regarded as part of the anti-collinear function J qq c in the second line with small anti-collinear momentum fraction r. Removing the quark from S NLP leaves the single-particle soft function S (g) , hence S NLP → S (g) and J

Bare factorization theorem
The distribution of an event shape e in e + e − annihilation into hadrons through a photon with virtuality Q 2 = q 2 is with the leptonic tensor and the electromagnetic current J µ =ψγ µ ψ, for simplicity of a single quark flavour with electric charge e ψ = +1. For a given final state X,ê(X) returns the value of the event shape variable. For small τ , and the thrust distribution is related to the two-hemisphere invariant-mass distribution by Eq. (5) is the starting point of the perturbative computation in terms of partonic final states X. The LP factorization theorem for the resummed event shape in the two-jet limit has been known for a long time [1,2]. The resummation of NLP doublelogarithmic corrections to the quark-antiquark back-to-back jet configuration, which is present already at LP, has been derived recently [9]. Defining the left hemisphere to be the collinear one and employing the SCET framework, the hemisphere invariant-mass distribution for M 2 R , M 2 L Q 2 is expressed at LP as in terms of the convolution of the jet (anti-) collinear functions with the hemispheresoft function, multiplied by a universal hard function [7,[23][24][25]. 3 The distribution is normalized to the leading-order total cross section with the number of colours N c = 3. Factorization theorems at NLP are more complicated due to multi-local convolutions, which further exhibit endpoint divergences, when expressed in terms of renormalized hard, collinear and soft functions [15].
Before turning to endpoint factorization, we provide in this section the SCET derivation of the factorization formula for off-diagonal gluon thrust, which has already been outlined in [14]. We introduce two light-like vectors n µ ± , satisfying n + n − = 2 and pointing into the directions of the back-to-back jets. The vector n µ − defines the direction of the collinear modes of the gluon jet. (Anti-) collinear modes therefore have large components n + p c ∼ Q (n − pc ∼ Q). The transverse momenta in the jet are of order The momenta of soft particles scale uniformly, p s ∼ τ Q. In the adopted reference frame, q µ ⊥ = 0. We further note that since gluon thrust vanishes at LP, there are no kinematic NLP corrections, and (7), (8) continue to hold. After integrating out the hard modes, the electromagnetic current matches tō in SCET. Here A c⊥ν denotes the gauge-invariant gluon field operator dressed with Wilson lines, which in light-cone gauge is related to the gluon field by A c⊥ν = g s A c⊥ν . The omitted terms are higher than NLP corrections and purely gluonic "flavour-singlet" operators. For simplicity, we neglect the mixing of qq into gg operators in the second line. Singlet terms cause well-understood technical complications but do not touch the essence of the factorization discussed in this work. The originally local current is now represented by light-ray operators. The first line corresponds to the LP hard vertex, which emits a back-to-back quark-antiquark pair. The second line of (11) represents the O(λ) suppressed "B-type" SCET operator, which directly produces the configuration relevant to gluon thrust, a collinear gluon and an anticollinear quark-antiquark pair. The first line can nevertheless contribute to gluon thrust through a time-ordered product with the O(λ) suppressed SCET interaction [26,27] which converts a collinear quark into a collinear gluon and a soft quark. Here and we further define the scalar quantity x ∓ ≡ n ± ·x/2. The gluon jet now recoils against a quark-antiquark pair in a highly asymmetric configuration, in which the antiquark carries almost all the momentum of the jet. 4 The Dirac structures in the power-suppressed B-type operator are given by We note that boost invariance together with the projection property / n − χ c = 0, / n + χc = 0 implies that only the transverse components of the vector current contribute.
Eq. (11) is the starting point for deriving the factorization formula and the main steps are fairly standard. After the collinear field redefinition with the soft Wilson line [28] for anti-collinear fields), the collinear, anti-collinear and soft fields are decoupled at LP. The final state |X = |X c |Xc |X s consists of a collection of corresponding modes, and hence the matrix element factorizes into suitably defined collinear, anti-collinear and soft functions. Performing these steps on the two terms in (11), more precisely the time-ordered product with (12) for the first, gives rise to an "A-type" and a "B-type" term in the factorization formula, which we discuss separately next. Note that due to the different field/mode content in the final state (soft quark or not), there is no interference between the two terms when squaring the amplitude.

A-type term (soft quark)
This term stems from the square of the matrix element summed and integrated over all possible final states. As mentioned above, there are two A-type terms, because either the quark or the anti-quark in the current can transfer its momentum to the collinear gluon. Due to the different final state, the two terms do not interfere. In (15) we consider explicitly only the first contribution. We comment on the second at the end. To separate the above matrix element into collinear, anti-collinear and soft factors, we perform the field redefinitions [28] Inspection of (11), (15) shows that the following hard, (anti-) collinear and soft functions appear.
Hard function This term involves the hard function C A0 (t,t) of the LP A-type current. We define the momentum-space coefficient It depends only on the product Q 2 = n + p c n − pc of the large components of the collinear and anti-collinear momenta. At LO, we find C A0 (Q 2 ) = 1 + O(α s ). The two-loop result can be found in [29].
Anti-collinear function The L ξq Lagrangian insertion acts in the collinear sector, hence the anti-collinear sector is unaffected. After squaring the amplitude and integrating over the anti-collinear final state, we obtain the LP anti-quark jet function 1 2π At LO, we have J . This object is well understood and computed up to the third order in α s [30]. Here and below we use small Latin letters for colour and Greek letters for Dirac spinor indices.
Collinear function This function picks up the collinear fields from the L ξq insertion and represents a new non-local jet function that appears at NLP. Defining the non-local operator where y − = n + · y/2, its jet function defined as in the LP case (19) by 1 2π The factor 1/g 2 s has been introduced for convenience, such that the leading term is O(α 0 s ). The Dirac matrix / n + /2 α α on the left-hand side ensures that only a single Dirac structure can appear on the right-hand side, and avoids a discussion of evanescent structures that would otherwise be present. In the derivation of the A-type term, this factor will arise from the decomposition of the soft function, see (24) below. 5 At LO in α s , the collinear final state X c consists of a single gluon, 6 and only the unhatted jet function is non-vanishing: Soft function Collecting the soft Wilson lines from the soft-decoupling field redefinition, as well as the soft quark field from the Lagrangian insertion (12) leads to the following definition of the NLP soft quark function: where x − is defined below (13), and contains the measurement function on the soft final state and the final-state sum pertaining to the two-hemisphere invariant mass distribution. With this definition, the left hemisphere contains the collinear gluon jet, and the right hemisphere is on the anticollinear side. The dots in (24) denote terms proportional to / n − , which do not contribute due to the projection property of collinear fields in other parts of the process. The prefactor is chosen for convenience, such that the soft function starts at O(α s ). The soft function for thrust is obtained from There is an analogous definition with obvious modifications for the NLP soft anti-quark function relevant to the second A-type contribution. The colour decomposition is chosen such that at leading O(α s ), S NLP (l + , l − , ω, ω ) = 0, while The two-hemisphere NLP soft-quark function is introduced here for the first time, but we take note of the calculation of the NLP soft function that appears for soft gluon emission in the Drell-Yan process [10], which is already known to O(α 2 s ) [10,31].
Factorization formula With these definitions, we obtain the soft-quark contribution to gluon thrust in the two-jet region in the form where is a d-dimensional factor defined such that f (0) = 1 in four dimensions.
To turn the convolution in l + , l − into a product, it is convenient to consider the Laplace transform of the two-hemisphere invariant-mass distribution. In general, the Laplace transform of collinear functions with variable p 2 are defined with ∞ 0 dp 2 e −sp 2 /Q , while for soft functions we use The total A-type term is twice the above (28), (31), since there is an equal contribution from the soft anti-quark phase-space region.

B-type term
The B-type term is obtained under the assumption that the intermediate propagator in Fig. 1 is hard, which leads to the B1 operator in the matching equation (11). We then need the square of the matrix elements summed and integrated over all possible final states. Inspection of (11), (32) shows that the following hard, (anti-) collinear and soft functions appear.
Hard function There are two relevant operators, i = 1, 2, which differ by their Dirac structures (14). We introduce the Fourier transforms of the corresponding hard functions. The B-type operators contain two fields (quark and anti-quark) in the same collinear direction. The momentum-space coefficient depends on the product Q 2 = n + p c n − pc of the large components of the collinear and total anticollinear momentum, and r, which denotes the fraction of total anti-collinear momentum carried by the quark. Accordingly,r = 1 − r is the momentum fraction of the anti-quark. At the leading order in α s , we find 7 From the transformation of (11) under charge conjugation, we find that charge conjugation invariance implies that C B1 2 (Q 2 , r) = −C B1 1 (Q 2 ,r) to all orders in the strong coupling. We notice that the matching coefficient C B1 2 ) is singular in the limit r → 0 (r → 1), in which the quark (anti-quark) becomes soft. From this observation, we can already anticipate a relation to the A-type contributions.
Further insight can be obtained by a helicity consideration. To this end, we separate the electromagnetic current in (11) into a left-and right-handed piece by inserting 1 = P L + P R . In four space-time dimensions we have 8 where projects on the left (−1) / right (+1) helicity of the collinear gluon in the B1 operator (11). Since helicity is conserved and amplitudes with gluons of different helicity do not interfer (see (41) below), the virtual photon and gluon always have the same helicity. Focusing on the negative helicity case for definiteness, we conclude from (36) -(39) that the operator with Dirac structure Γ µν 1 produces a left-handed out-going quark, while for Γ µν 2 it is right-handed. It follows that the two operators cannot interfer in the non-singlet channel, where theχc(t 1 n − )χc(t 2 n − ) fields in the operator in (11) cannot be contracted among themselves. This strictly holds only in four space-time dimensions. In other words, a non-singlet interference term must be "evanescent" and vanish as → 0.
The limit r → 0 in the B1 matching coefficient implies that the momentum of the out-going quark in the left diagram of Fig. 1 goes to zero, and the intermediate quark 7 We recall that the definition of A µ c⊥ absorbs the factor g s from the tree-level diagram. 8 Convention: 0123 = +1 propagator goes on-shell. The helicity argument above implies that only C B1 1 (Q 2 , r) corresponding to Γ µν 1 can be singular in this limit. For definiteness, consider again the case of the negative helicity gluon and virtual photon. The nearly on-shell intermediate propagator implies that the B1 operator factorizes into γ * → q +q and q → q + g. As r → 0 the outgoing quark in the latter process has zero momentum, and the incoming quark and outgoing gluon have the same large collinear momentum. Since g has negative helicity, the incoming quark must be left-handed, since a helicity change by 3/2 is not possible. By helicity conservation of the strong interaction, the final soft quark must then also be left-handed. As shown above, this singles out Γ µν 1 . The situation is reversed in the soft anti-quark limit r → 1 and implies right-handedness, which singles out Γ µν 2 . We conclude that, to all orders in perturbation theory, only C B1 1 (Q 2 , r) (C B1 2 (Q 2 , r)) can be singular in the endpoint r → 0 (r → 1).

Collinear function
The collinear function is simply the standard LP gluon jet function, defined as 1 2π At LO, we have J (g) c (p 2 ) = δ + (p 2 ). This object is well understood and known up to the third order in α s [32].
Anti-collinear function The anti-collinear fields consist of the quark-anti-quark pair, which recoils against the gluon and is therefore in a colour-octet state. The corresponding jet function can be defined as for the gluon, but in terms of the composite field At leading order, r can be interpreted as the momentum fraction of the out-going quark. The Dirac structures refer to (14). Except for the colour and Dirac matrix, the operator coincides with the one that enters the definition of light-cone distribution amplitudes of mesons. However, instead of the vacuum-to-meson matrix element, here we need the inclusive jet function which is a NLP object. It is easy to see that the four possible values of i, i = 1, 2 give rise to only two independent jet functions, which we decompose as above into a diagonal and off-diagonal term in ii . At lowest non-vanishing order O(α s ), the Xc = qq final state results in The hatted function, which describes the interference of the two B1 operators, is given as this order by J (p 2 , r, r ), and vanishes in four dimensions, as expected from the helicity discussion above. In principle, it should be possible to remove this evanescent anti-collinear function by a finite, non-minimal counterterm related to its mixing into the "physical" jet function. This relies crucially on a renormalized factorization theorem, i.e. being able to take → 0 at the level of the individual collinear, soft etc. functions. Since for the time being, we continue to work with d-dimensional quantities, we keep J qq(8) c (p 2 , r, r ) in the following. We will see later that it is irrelevant for the discussion of endpoint divergences.
Soft function Performing the soft decoupling field redefinition on (32) results in the field productχc The Wilson lines can be combined into the adjoint Wilson lines Y n ∓ using After squaring the matrix element and summing over the soft final state, we obtain the LP two-hemisphere soft function in the adjoint representation: Factorization formula Inserting these matrix element definitions, we obtain the factorization formula for the B-type term In Laplace space: Only the first term in the curly brackets proportional to δ ii is non-vanishing at leading O(α s ). 9

Tree-level evaluation
Inserting tree-level results for all functions given above, the two terms in equations (28) and (48) reduce to the following expressions. The A-type soft-quark term is given by where the pole in arises from the logarithmic divergence of the integral ∞ l + dω/ω for large values of the soft function variable. The total A-type contribution is twice the above after adding the soft-antiquark term. Similarly, the B-type term is given by In this case, the singularity arises from a logarithmic divergence in the integral over the momentum fraction r as r → 0 and r → 1. Adding both terms and taking the limit → 0, we find The poles cancel in the sum of both the A and B-type contributions in equations (50) and (51) (once we account for the soft-antiquark contribution to the A-type term) as it should be, as the off-diagonal gluon-thrust is infrared safe. After converting the twohemisphere invariant mass distribution to the thrust distribution according to (8), we reproduce the coefficient −α s C F /π of the ln τ term in [33]. The single logarithm arises from the dimensionally regulated convolution integrals, which are logarithmically divergent in d = 4. For resummation we would like to define renormalized hard, soft and (anti) collinear functions and take the limit → 0 before performing the convolution integrals. Proceeding in this way, however, the convolution integrals are ill-defined. In the following sections we explain how this problem can be solved.
To prepare this discussion we make the following observation. The integrand of the two-dimensional convolution in ω, ω in (50) in the limit of large ω, ω (where the divergence arises) is given by In the B-type term (51), the divergence of the convolution arises from small momentum fraction. The integrand of the r, r integrals in the limit of small r reads which is identical to the previous expression provided we identify r = ω/Q, r = ω /Q. A similar agreement holds for the identical soft anti-quark contribution to the A-type term and the r → 1 limit of the B-type term (51). From the heuristic discussion of Sec. 2 it is evident that this agreement is not a coincidence, but holds to all orders in the expansion in α s .

Endpoint factorization
Eqs. (28), (48) provide factorized expressions for the two-hemisphere mass distribution in the two-jet limit in the gluon-thrust region from which the NLP logarithms can be computed to any logarithmic accuracy order by order in α s by evaluating the convolutions of the d-dimensional hard, (anti-)collinear and soft function. The limit → 0 must be taken at the end, because the convolutions in the A-type term are divergent when ω and ω both tend to infinity. In the B-type term this happens when r and r both approach 0 or 1. This form is not quite sufficient to sum the logarithms to all orders in α s . To apply the standard SCET renormalization-group method to the hard, (anti-)collinear and soft factors, one must renormalize them and take the limit → 0 before the convolution. In this section we employ the coincidence of the integrands in certain asymptotic limits, which should hold on physics grounds and has been demonstrated above at lowest order, to factorize and rearrange the endpoint contributions such that the convolution integrals are finite. Similar methods to derive endpoint factorization at the level of factorization formulas with functions defined by matrix elements of operators have been derived before only for exclusive B decays to P-wave charmonia [34], which involves SCET and nonrelativistic QCD, and Higgs decay to two photons [22], which is a SCET II process. These are amplitude-factorization problems. The present case of gluon thrust is somewhat different as it is a SCET I problem, and involves inclusive collinear and soft functions defined at cross-section level, similar to [18]. Nevertheless, the main mechanism that achieves endpoint factorization is similar to the above-mentioned cases. Interestingly, it partly involves functions, which are the non-abelian generalization of some that appear in H → γγ.

Soft-collinear limit of the B1 matching coefficients
An important ingredient in this discussion is the factorization property of the coefficient function of qqg SCET B1 operators, when the light-cone momentum fraction carried by the quark or anti-quark becomes small. The endpoint divergence in the B-type contribution originates from the singular behaviour of the C B1 i matching coefficients in this limit, which in turn arises because the intermediate quark or anti-quark goes on-shell, see Fig. 1.
By its original definition, the matching coefficient is a single-scale, hard function. For r → 0 (soft quark) and r → 1 (soft anti-quark), it becomes a two-scale object, which can itself be factorized according to [18] where the two scales Q 2 and rQ 2 are now separated into the hard matching coefficient of the leading-power A0 SCET current, and a new coefficient D B1 (rQ 2 ), which depends only on the endpoint-scale √ rQ. We recall that C B1 1 (Q 2 , r) is regular in the soft antiquark limit r → 1. An analogous factorization holds for the second B1 operator in the limit r → 1: r → 0q D B1c s q p 2 Figure 3: Graphical representation of the factorization (55) of the B1 matching coefficient in the soft-quark limit.
Since the two cases are related by r ↔r (plus a global minus sign), the following is phrased for the soft-quark limit only. A graphical representation of (55) is shown in Fig. 3.
For r → 0 the outgoing quark represented byχc field in the second line of (11) becomes soft-anticollinear. The soft-collinear coefficient D B1 (p 2 ), which appears in (55), (56), can be defined by matching the time-ordered product [18] with p 2 = (p c + p sc ) 2 = rQ 2 and the leading soft-quark Lagrangian [26,27] withq s replaced byq sc . This relation becomes evident by comparing the second diagram in Fig. 2 to the right-hand side in Fig. 3. From its definition in terms of a time-ordered product it is clear that the matching coefficient D B1 (p 2 ) is a universal function that will appear in different processes involving soft quark emission.
The leading double-logarithmic resummation of D B1 (p 2 ) to all orders in α s has been derived in [18]. It exhibits the all-order (C A − C F ) n colour coefficient, which seems to be characteristic for soft quark emission [14,18,35]. Interestingly, the same coefficient enters the ggH amplitude with a bottom-quark loop and was recently computed at the two-loop level [36], together with its one-loop evolution kernel, which was obtained from renormalization-group invariance of the full ggH amplitude.
The D B1 (p 2 ) coefficient as well as its evolution equation can also be obtained from the corresponding B1 operator coefficients and anomalous dimension by taking the limit r → 0 and provide the relevant one-loop expressions. Since the anti-collinear part of the B1 operator has the same field content as the operator relevant to light-cone distribution amplitudes, the extraction of the anomalous dimension is similar to the derivation of the asymptotic kernel for the QED light-meson light-cone distribution amplitude [37], see also [22]. We refer to App. A for this derivation, which results in with in agreement with [36]. When the gluon field is replaced by a photon, the corresponding abelian version of D B1 (p 2 ) is identical to the jet function, which appears at LP in the factorization of the radiative semi-leptonic B decay B → γ ν [38] and in the NLP endpoint factorization of the photonic Higgs decay H → γγ amplitude through a bottom-quark loop [22]. In the abelian case, the two-loop evolution has been inferred from renormalization-group consistency of the B → γ ν observable [39]. The direct computation of the one-loop evolution kernel of this jet function can be found in [40].

Endpoint factorization consistency conditions
As discussed in Sec. 2 and checked explicitly above at O(α s ), we expect the integrand of the A-and B-type terms of the factorization formula to have identical asymptotic limits to all orders, which is a prerequisite for endpoint factorization. More precisely, the limit of anti-collinear momentum component n − · pc = rQ, r Q → 0 in the B-type term, matches the limit n − · k = ω, ω → ∞ of the corresponding soft momentum component in the A-type term. 10 We start from the expressions (31), (49) in Laplace space, in which the integrand of the ω, ω and r, r integrals involves no further integrals and the consistency relations are algebraic. The coincidence of asymptotic limits implies that the following must hold: Figure 4: Factorization of the radiative collinear function J c (p 2 , ω, ω ), which appears in the A-type term, in the limit ω, ω → ∞.
The left-hand side of this equation represents the leading asymptotics of the B-type term (49) for r, r → 0, already accounting for the fact that this arises only from the square of the amplitude of the first B1 operator, since the coefficient of the second is non-singular as r → 0. The right-hand side is one half of the total A-type term, which we can identify with the soft quark contribution. An analogous equation with r, r →r,r and C B1 1 → C B1 2 applies to the limit r, r → 1. The right-hand side of (62) contains two terms in the curly brackets, which correspond to different colour structure in the jet and soft function, (21) and (24), respectively. We shall now show that endpoint factorization consistency implies that the hatted term in (62) is subleading as ω, ω → ∞, and can be dropped.
We recall from (21) that J c (p 2 , ω, ω ) is an inclusive radiative jet function of the nonlocal, time-ordered product operator defined in (20). Graphically, we may represent it by the diagram to the left of the arrow in Fig. 4, where the dashed external lines to the left and right of the cut denotes the soft momentum of the composite field A ⊥c χ c in the time-ordered product, which is carried away by the soft quark. We now claim that the leading asymptotic behaviour proportional to 1/(ωω ) for large ω, ω requires that the diagram is "one-collinear-particle reducible", as displayed to the right of the arrow in Fig. 4. Then colour conservation implies that the collinear line must be a colour-octet, that is a gluon, and hence only the unhatted colour term (21), which defines J c (p 2 , ω, ω ) contributes to the asymptotic behaviour, as was to be shown.
Suppose the diagram was not one-collinear-particle reducible, such as the one-loop example displayed in Fig. 5. Endpoint factorization consistency requires that the large ω, ω must match a corresponding B-type contribution, obtained by making the dashed lines representing external soft quark anti-collinear. This turns internal collinear propagators into hard propagators, and the entire diagram into a contribution to the matching coefficient of a hard operator with external collinear and anti-collinear fields. The key observation is that if the diagram was not one-collinear-particle reducible, the corresponding hard operator would contain more than one external collinear field, in contradiction with the B1 operators available at NLP, which contain only a single collinear gluon field, see (11). Consistency therefore requires that the asymptotic behaviour takes the diagrammatic form shown in Fig. 4. This proves not only that the hatted jet function is irrelevant for endpoint factorization, but further that J c (p 2 , ω, ω ) factorizes into the s s c c c c c c Figure 5: Example of a one-collinear-particle reducible diagram, which does not contribute to the large-ω ( ) behaviour of J c (p 2 , ω, ω ).
where the function D B1 (p 2 ) is the same as the one that appears in the factorization of the hard B1 operator coefficient (55). Eq. (63) is the first endpoint factorization consistency condition. Inserting the Laplace-transformed version of this relation together with (55) for the B1 hard function into (62) immediately implies a similar consistency relation for the NLP soft quark function, The same identity holds with r, r →r,r . Relations (I) and (II) formalize the heuristic picture developed in Sec. 2. For large ω, the soft quark field in S NLP becomes anticollinear,q s Y n − →χc, hence it moves from S NLP to the anti-collinear function J qq(8) c . Removingq s Y n − from S NLP leaves the leading-power soft function S (g) , while adding it asχc to J . Overall, this results in , which is relation (II). At the same time, the quark fields in the A-type collinear function (20) become highly off-shell, which removes them from J c (p 2 , ω, ω ), leaving only the collinear gluon, and consequently C A0 turns into C B1 .
c , which is relation (I).

Endpoint factorization formula
We are now in the position to derive the endpoint-finite factorization formula. For its concise formulation, we employ the double-bracket notation introduced in [21] to denote the asymptotic behaviours of the various functions. The precise definitions are as follows: In functions of ω, ω , rescale ω → κω, ω → κω and take κ → ∞. Then where δ(ω − ω ) counts as κ −1 . The right-hand side of the previous equation equals the right-hand side of the consistency relation (I). Similarly, in functions of r, r , rescale r → rκ, r → r κ and take κ → 0, or the corresponding rescaling is applied tor,r . Which of the two is meant, will be indicated by the subscript 0 or 1 on the double bracket. Then With this definition the right-hand sides of these equations are given by (55) and (56), respectively. Finally, we have due to the symmetry r ( ) ↔r ( ) . In the d-dimensional expressions there are -dependent powers of κ, which turn into logarithms in the renormalized functions. Thus counts as an infinitesimal variable for the purpose of endpoint κ power-counting.
To implement the rearrangement of endpoint-singular terms we start with the scaleless integral which vanishes in d dimensions. We split this integral in two terms I 1,2 , which can be done in two ways, illustrated in Fig. (6): (1) ω and ω smaller than an endpoint factorization parameter Λ (integral I 1 ), and the complement region I 2 .
(2) ω or ω smaller Λ and the complement region. Since the mixed regions where one variable is larger and the other smaller than Λ are not endpoint-singular, both ways could be employed for the following considerations, although the expressions take a different form. We will use the second version. In this case the complement region is ω, ω > Λ and the double-bracket asymptotic behaviour can be used for functions of ω, ω in the A-type term. 12 In both versions I 1 + I 2 = 0. The endpoint rearrangement now consists of subtracting I 1 from the B-type term and I 2 from the A-type term. The subtracted expressions are now separately endpoint-finite, but depend on Λ. However, as long as no approximations are made, the Λ dependence cancels exactly between the two terms.  (2) [right] of the split of (70) into I 1 + I 2 according to the correspondingly indicated regions in the ω − ω plane as described in the text below (70). In the overlap region in green the asymptotic behaviour of the A-and B-type term must agree.
We begin with the A-type contribution (31), which accounts for soft-quark emission, to the factorization theorem and subtract from it the complement region I 2 , ω, ω > Λ, of the integral (70), resulting in For ω, ω → ∞, we use the endpoint factorization condition (I) from (63) to rewrite the jet function in the first term of the integrand. In this limit it is also justified to replace the soft function by its asymptotic form (65). Thus the first and second term in the integrand cancel. The third term is subleading in this limit, and the convolution integrals over ω and ω are now endpoint-finite, as the logarithmically divergent part of the integrand in the large ω and ω region has been removed. To make this explicit, we use (63) in (71), and obtain If we further choose Λ 1/s R , 1/s L , we may simplify the previous equation to where the equality signs hold up to corrections of O(1/(s L Λ), 1/(s R Λ)). It is then understood that in evaluating (73), terms suppressed by powers of Λ are dropped. This is simply (31) with integration region ω, ω > Λ removed, see Fig. 6. The remaining part I 1 of the integral (70) can now be combined with i = i = 1 part of the B-type term (49). Similarly, we proceed for the A-type soft-antiquark contribution and the i = i = 2 part of the B-type term after using the symmetry under the exchange r ↔r. The term i = i is endpoint-finite so it remains unaffected.
Beginning with the i = i = 1 term, we have We substitute ω = rQ and ω = r Q in the subtraction term, and employ (55) to express the D B1 functions in terms of the asymptotic limits of the B1 matching coefficients.
Next, we make use of the endpoint factorization condition (II) to express the asymptotic soft-quark function in terms of the B-type jet function J qq(8) c (s R , r, r ) in the asymptotic limit to obtain where now the integrals over r, r are evaluated from 0 to ∞. The convolution integrals are now convergent, but they must be evaluated after the integrands are combined. Provided Λ Q, this expression can be simplified to up to corrections of O(Λ/Q). The first integral in the curly brackets is simply (49) with integration region r, r < Λ/Q removed, see Fig. 6, which is endpoint-finite. The second integral is a finite left-over from the subtraction term. For completeness, we provide the term i = i = 2. In the bare factorization formula (49) we change integration variables tor,r to map the singular point r = r = 1 tō r =r = 0. The subtraction term is then obtained after mapping ω =rQ and ω =r Q and following the same steps as before we find which coincides with (75) up to the expected substitutions. Similar simplifications hold when Λ Q is assumed. The total endpoint-subtracted B-type is the sum of (75), (77) and the mixed term i = i in (55).
By construction, the dependence on Λ cancels in the sum of the A-and B-type contributions. In App. B we provide the corresponding version of the endpoint-rearranged factorization formulas, when (70) is split according to the first version mentioned below this equation, that is, for ω and ω smaller than Λ.

Renormalization-group equations
In this section we collect the renormalization-group functions, which are known for some of the hard, soft and (anti-) collinear functions of the two terms in the factorization formula, and infer the others from renormalization-group consistency, namely the condition that the renormalized integrands of the A-and B-type term must be separately independent of the scale µ. The endpoint-factorization conditions then provide further relations between the anomalous dimensions in the respective large-ω and small-r limits. We then obtain the solutions to leading-logarithmic (LL) accuracy, which includes running-coupling effects.

RGEs for the A-type functions
In the A-type term, the hard function and the anti-collinear function are LP objects. Their evolution is therefore well-known. The evolution of the two NLP objects, the collinear function and the soft function, is currently unknown. We cannot reconstruct both RGEs from just consistency arguments since we miss two anomalous dimensions.
Instead, we will use the endpoint factorization condition (I) from (63) to derive the evolution of the collinear function J c (p 2 , ω, ω ) for large ω, ω from the known RGEs for the LP gluon jet function and the D B1 coefficient. Through consistency, we can then derive the asymptotic evolution of the soft function. This allows us to resum the A-type functions in the endpoint region, which is sufficient for LL accuracy.
Hard function The evolution of the LP hard matching coefficient C A0 is well known [41]. It obeys the RGE where −Q 2 in the logarithm should be read as −Q 2 − iε. The cusp anomalous dimension γ cusp and the hard non-cusp anomalous dimension γ A0 are given by The one-loop cusp anomalous dimension sums the leading logarithms. Non-cusp terms become relevant only at NLL accuracy. 13 Dropping them, we obtain the LL solution where µ h ∼ Q denotes the hard initial scale, such that C A0 (Q 2 , µ 2 h ) free from large logarithms. The functions S (ν, µ) and A γ i (ν, µ) are defined as [42] S (ν, µ) = − αs(µ) The QCD beta-function is Anti-collinear function The anti-collinear function appears already at LP. In Laplace space, it obeys the local RGE [41] where the non-cusp term is given by At LL accuracy, the RGE is solved by where µc denotes the anti-collinear initial scale.
Collinear function The collinear function J c (p 2 , ω, ω ) is a new NLP object and its anomalous dimension is currently unknown. In the asymptotic region of large ω, ω , it factorizes into the LP gluon jet function and two D B1 matching coefficients as shown in (63). We therefore derive its asymptotic evolution from the LP gluon jet function, which renormalizes in Laplace space as [43] with the gluon non-cusp term given by and the evolution of the D B1 coefficient. Making the cusp term explicit, (60), (61) or (161), (162) imply the non-local RGE where the non-cusp anomalous dimensionγ D (ω, ω) is given by the second line of (162). The first endpoint factorization consistency condition (63) now implies the following asymptotic RGE for the collinear function J c , There is no mixing of the unhatted collinear function into the hatted collinear function, which is irrelevant asymptotically. At LL accuracy, we drop the non-cusp terms (except for those proportional to β 0 ) and obtain the solution We introduce two initial scales, µ c and µ cΛ , connected to the two different cusp logarithms in (91), allowing for the possibility to choose different initial scales for the gluon jet function and D B1 coefficient, which will be used in Sec. 6.1.

Soft function
The soft quark function S NLP (l + , l − , ω, ω ) appears first at NLP and its complete RGE is currently unknown even at lowest order. We can obtain its asymptotic evolution for large ω, ω from RGE consistency. For fixed values of ω, ω and r, r , respectively, the integrands of the A-and B-type term must be individually RGE-invariant. Imposing this requirement on the A-type term, we obtain the asymptotic RGE for the soft function, Consistency requires that there is no mixing of the unhatted soft function into the hatted soft function. At LL accuracy, the asymptotic RGE is solved by where again we allow for two separate initial soft scales, µ s and µ sΛ . The lowest-order initial condition is proportional to α s , which we choose to evaluate at µ s . Together with the prefactor α s (µ)/α s (µ s ), which arises from the β 0 -term in γ J (g) , this results in an overall factor of α s (µ). The point of keeping the non-cusp β 0 -term in γ J (g) is that it renders the result independent of the initial scale -choosing to evaluate α s at µ sΛ instead, we would obtain the same resummed soft function.

RGEs for the B-type functions
In the B-type contribution, the collinear function and the soft function are known LP objects while the hard functions C B1 i and the anti-collinear function J qq(8) c (p 2 , r, r ) are NLP objects. The evolution equation of the hard matching coefficients is of the form discussed in [44,45] for B1 SCET operators, hence we can derive the anomalous dimensions for the new anti-collinear function from consistency, i.e. the scale-independence of the entire B-type term for given r, r . However, the leading logarithms to gluon thrust arise only from the endpoint region r, r → 0, 1, and therefore we will only need the equations for the asymptotic hard and anti-collinear functions.
Collinear function The anomalous dimension of the collinear function is well-known since it appears in LP factorization theorems. The RGE equation and anomalous dimension have already been given in (88) and (89), respectively. The LL solution to the RGE reads where µ c denotes the collinear initial scale.
Soft function The LP two-hemisphere soft function for gluons, which is defined with adjoint Wilson lines, obeys the evolution equation in Laplace space. The non-cusp part of the anomalous dimension is given by [46] At LL accuracy, the solution to the RGE reads Hard function Following the calculations of [44,45], the specific anomalous dimension matrix for the two coefficients C B1 i , i = 1, 2 needed here can be inferred from App. A.2. For the flavour non-singlet case, which is assumed here, ∆ F can be set to zero there, and the two coefficients evolve independently, with the same anomalous dimension (see (146)) where γ(r, s) = α s 2π The non-cusp part of this RGE contains a hidden large logarithm of r (r), which must be extracted and combined with the cusp part when one wants to sum the hard function for small r orr, see App. A.2. This results in the evolution equation (60) for the D B1 coefficient. The evolution equation for the asymptotic C B1 i coefficients then follows from (55), (56). Since the leading logarithms to gluon thrust arise only from these asymptotic regions, we provide the RGEs for the asymptotic coefficients in the following.
The asymptotic B1 coefficients factorize into a product of A0 coefficient and the D B1 coefficient. Combining (78) and (90), we find The LL solution reads where as in the case of similar A-type RGE solutions, we allowed for two different initial scales, µ h and µ hΛ , connected with the different cusp logarithms.
For each of the two endpoint limits r → 0, 1, we only need one of the two B1 coefficients since the other one is regular. The RGE equation and solution for C B1 2 (Q 2 , r, µ 2 ) is obtained from the above by replacing r →r, and similarly forr including the integration measure dr.
Anti-collinear function The evolution of the hard, collinear, and soft functions is known, hence we can derive the evolution of the anti-collinear function through consistency. We demand that the B-type contribution before integration over r, r is RGEinvariant on its own, which allows us to derive the RGE and anomalous dimension of the anti-collinear function in the form where the non-cusp part of the B1 anomalous dimensionγ B1 is given bŷ As for the hard function, the non-cusp part of the full anomalous dimension contains hidden leading logarithms coming from the endpoint region. We therefore use asymptotic anomalous dimensions for the LL resummation. RGE-consistency yields in the two endpoint regions At LL accuracy, the asymptotic RGE is solved by where µc and µc Λ denote the anti-collinear initial scales. As for the A-type soft function, we may evaluate the factor of α s in the fixed-order initial condition either at µc or µc Λ , and µc has been chosen above. This choice does not affect the overall resummed anticollinear function since the ratio of couplings in front always guarantees a global factor of α s (µ). Once again the RGE and solution for the anti-collinear function for r → 1, J qq (8) c,LL (s, r, r , µ 2 ) 1 , are obtained from the above by replacing r →r, and similarly forr including the integration measure dr as well as for the corresponding primed quantities.

Resummation
With the individual functions resummed to LL accuracy, we proceed to the resummation of the leading logarithms of gluon thrust at NLP. All functions are renormalized, so from now on we set d = 4.
We begin by rearranging the endpoint-subtracted expressions (72) and (75) for the A-and B-type term, respectively. In (72), we rewrite where σ(ω, ω ) is an auxiliary function, which equals unity whenever ω ( ) s R , ω ( ) s L 1 and which goes to zero for ω ( ) → 0. We find In this expression, the first line in curly brackets is the subtraction term, which contains the endpoint divergence, if Λ → ∞, and therefore results in a logarithmically enhanced contribution from the ω, ω convolution. The next term in square brackets resembles a "plus-distribution" type term -this combination has no support for large ω, ω and does not produce an extra logarithm. Finally, the last term in curly brackets involving the hatted soft function is regular by itself and does not require subtraction. After this rewriting, only the first line in curly brackets needs to be kept at LL accuracy. The introduction of the auxiliary function σ is necessary to eliminate a spurious singularity for small ω, ω , which arises if one wants to integrate the subtraction and plus-distributionlike term separately, and which cancels between the two. Similarly, we write (75) in the form The first integral in curly brackets represents the subtraction term, which contributes to the leading logarithms and cancels the Λ-dependence of the corresponding A-type term. The second integral with support in [0, 1] is of the plus-distribution type and can be safely integrated to r, r = 0. This term can be dropped at LL accuracy. An auxiliary function is not required for the B-type term. We remind the reader that there is an identical soft anti-quark A-type term and the i = i = 2 term, which is identical to the i = i = 1 term due to the symmetry in r ( ) ↔r ( ) . The i = i interference terms can be dropped at this point for LL accuracy.
The previous two equations are general and provide a suitable starting point for the resummation of NLP logarithms beyond the leading ones. Since some of the nonasymptotic anomalous dimensions in the plus-distribution-like terms are not yet available, the following focuses on the leading logarithms, including a new set of next-toleading logarithms associated with the running coupling in the leading terms.

Counting of logarithms and proper choice of scales
Standard Sudakov exponentiation puts an observable into the form where λ = α s ln τ counts as O(1), while τ 1. The function g 0 is needed to sum the leading logarithms, g 1 is NLL etc. F collects the initial condition and has an expansion without large logs. The O(α s ) correction to the leading approximation of F is already a next-to-next-to-leading logarithmic (NNLL) effect. The large-logarithm counting primarily refers to the exponent. Figure 7: Graphical illustration of the resummation of the initial condition of the B1 operator for r 1 (left) and the LL-accurate approximate procedure employed here.
The integrands of (109), (110) are of this form when the solutions to the RGEs from Sec. 5 are used. However, one now has to take integrals of the form and the assumption underlying standard Sudakov resummation that ω = O(τ ) and r = O(1) is not justified, since either Λ τ Q or Λ Q, or both. The fact that the integrals are logarithmic for large ω and small r implies two modifications of the standard resummation scheme.
First, the logarithmic integrals combine to an extra logarithm promoting their logcounting by one order. Thus, inserting the LL-resummed RGE functions in the subtraction terms in (109), (110) will result in leading logarithms for the observables. On the other hand, inserting them into the plus-distribution type terms will result only in next-to-leading logarithms, since these terms are not sensitive to the large ω and small r region by construction and therefore lack one power of logarithm. This counting extends to higher logarithmic order in an obvious manner and justified the previous assertion that to LL accuracy, it suffices to focus on the subtraction terms in (109), (110).
Second, the fact that ω = O(τ ) and r = O(1) cannot be assumed can be seen as the need for resumming the initial condition of the RGE evolution. Consider for example the hard function C B1 1 (Q 2 , r, µ 2 ). Standard RGE summation implies that in the initial condition, one chooses the scale µ h ∼ Q, at which there are no large logarithms and evolves to µ Q. However, the initial condition for µ h ∼ Q may still contain an infinite series of ln r terms. While these normally count as O(1), they lead to a breakdown of the perturbative expansion in integrals such as (112), which receive contributions from r 1. These logarithms are summed by evolving with asymptotic evolution kernel for the D B1 coefficient from µ 2 ∼ rQ 2 to µ 2 h . At LL accuracy, combining this evolution with the standard evolution from µ h to µ, amounts to the expression (103) with split initial scales -µ h ∼ Q in C A0 , and the running scale µ 2 hΛ ∼ rQ 2 for the part that arises from the D B1 coefficient, whose natural scale is rQ 2 . The procedure is illustrated in Fig. 7. A similar reasoning applies to the choice of initial conditions for the anti-collinear function J qq(8) c (s R , r, r ) and the functions J c (s L , ω, ω , Q 2 , µ 2 ), S NLP (s R , s L , ω, ω , µ 2 ), which appear in the A-type term. Hence, we choose the initial scales in the hard, soft and (anti-) collinear functions to be of order of To ensure the cancellation of the rearrangement scale Λ between the A-and B-type term, the integrands of the subtractions terms must match, which requires which is satisfied by (113) with the identification rQ ↔ ω.
We now insert the resummed functions from Sec. 5, evolved from their respective initial scales to a common scale µ, into the integrands of (109), (110), and simplify the expressions employing the identities and the special case µ 3 = µ 1 . We also abbreviate A ≡ A γcusp . We find for the A-type integrand, and for the B-type, we find As it should be with all functions evaluated consistently to LL accuracy the above expressions are manifestly independent of the arbitrary renormalization scale µ. Note that if we would not have included the terms proportional to β 0 in the non-cusp parts of the anomalous dimensions, the only change would be that the overall prefactor of α s in the square bracket in (117) would be evaluated at scale µ s (or µ sΛ ) instead of µ c . Similarly, in (118) the prefactor would be evaluated at µc (or µc Λ ) instead of µ c . The scale at which the prefactor is evaluated formally contributes at the same order in the logarithmic counting as the LL terms in the argument of the exponential factors, when including the running coupling. Therefore, it is necessary to take into account the terms proportional to β 0 in the non-cusp parts of the anomalous dimensions to obtain a consistent LL result. In addition, we note that this also ensures that the final result is insensitive to a possible redefinition J /α s and C B1 1 → C B1 1 × √ α s , that corresponds to a reshuffling of the overall factor of α s from the anti-collinear to the hard function, and similar ambiguities in the definitions of the functions in the A-type term.
We note the similarity of the two expressions. In fact, enforcing the relations (114), both expressions are exactly equal up to an overall factor of Q with the identification ω ( ) = r ( ) Q. We thus sum the subtraction terms from (109), (110), and change variables from r to ω/Q in the B-type contribution (110). We also choose the auxiliary function to be of the form σ(ω, ω) = θ(ω − σ), where σ is of order 1/s R , or 1/s L . The two integrals to be added then combine as and we obtain where µ 2 hΛ ∼ ωQ and µ 2 sΛ ∼ ω/s R according to (113). We also added an overall factor of 2 to account for the identical soft anti-quark and i = i = 2 contribution. This is our main result for the LL resummed two-hemisphere invariant mass distribution in the NLP gluon thrust region. The Laplace-space gluon thrust distribution itself is obtained as which follows from (8) and the definition of the Laplace transformation.

Double-logarithmic limit
The double-logarithmic limit corresponds to evaluating the renormalization group functions at leading order in an expansion in α s : Setting scales as prescribed by (113), the LL result (120) simplifies to The scale of α s remains undetermined in the double-logarithmic approximation.
To obtain the thrust distribution we set s R = s L = s and choose σ = 1/(se γ E ). For the inverse Laplace-transformation from s to τ = 1 − T we use that the Laplace transform of τ a is Γ(1 + a)/(sQ) 1+a . Taylor-expanding this relation in a gives where the arrow denotes inverse Laplace-transformation. At the double-logarithmic level, the inverse Laplace-transformation is therefore simply accomplished by multiplying with sQ and substituting se γ E → 1/(Qτ ) resulting in This agrees with the previous result [14,18], which, however, was obtained through ddimensional factorization and consistency relations. 14

Leading logarithms and running coupling effects
The present approach based on renormalized functions satisfying standard RGEs allow us to go beyond the double-logarithmic limit. As a first application, we consider the LL approximation, which includes one-loop running-coupling effects. The RGE functions read and the strong coupling α s at a scale ν is given by in terms of the coupling at a reference scale µ. In the following α s without argument refers to α s (µ). These expressions are now to be used in (120), which upon implementing (113) (with "∼" replaced by "=" and s i by s i e γ E ) reads 15 This is our main result for the resummed two-hemisphere invariant mass distribution in the gluon-thrust region. The integral over ω can be performed only numerically. An approximate analytic expression can be obtained by expanding the one-loop running coupling to linear order in β 0 in the exponent and prefactor, resulting in where Once again, the thrust distribution in Laplace-space is obtained by identifying s R = s L ≡ s, in which case we also choose σ = 1/(se γ E ).
To display the effect of resummation and of running coupling effects, we show in the upper panel of Fig. 8 the gluon-thrust distribution 1 σ 0 dσ dτ in the leading O(α s ) approximation (blue-dotted line) and in the resummed double-logarithmic approximation (126) (light-blue band). We adopt Q = m Z = 91.1876 GeV corresponding to final states from hadronic Z-boson decay, and α s (m Z ) = 0.1179. The coupling is evolved with one-loop accuracy and n l = 5 according to (128). In the leading-order approximation we choose α s at the collinear scale √ τ Q. Since the scale of α s is not fixed in the double-logarithmic approximation, we vary it from the soft scale τ Q to the hard scale Q, which produces the shown band and a considerable uncertainty. For comparison, we also show (light-green band) the numerical inverse Laplace transform of the double-logarithmic approximation (124) in Laplace space. To this end, we compute where c > 0 is chosen to enclose the singularities of the Laplace-space distribution, excluding the one from the Landau pole of the running coupling, which corresponds to the "minimal prescription" of [47]. We also set σ = 1/(se γ E ) in (124). We observe a significant difference between the two bands, which must be attributed to formally nextto-leading logarithmic effects in the relation between Laplace and momentum space.
Finally, the upper plot shows (red line) the full leading-logarithmic result (129) (with the ω integral evaluated numerically for σ = 1/(se γ E )) after taking the inverse Laplace transformation (132) numerically. The LL result turns out to be very close to the DL one when the inverse Laplace transformation is computed in the same way. The green band ends at τ ≈ 0.03, since the DL result at the soft scale becomes sensitive to the strong coupling regime at smaller values τ . We also recall that perturbative resummation requires Qτ Λ QCD ≈ 0.5 GeV and τ 1. We can gain some insight on the importance of intrinsic next-to-leading logarithms by varying in the LL result the various matching scales around the values adopted in (129). For this purpose, we vary the three pairs of scales (µ h , µ hΛ ), (µ c , µc), (µ s , µ sΛ ) by a factor of 1/2 and 2 around their default scales. We then take the minimum and maximum values of the 15 possible combinations of (1, 1 2 , 1 2 ) etc. to compute the scale variation, excluding the cases where different scales are varied in different directions. For simplicity, we show this result for the normalized Laplace-space distribution Qs σ 0 dσ ds in the lower panel of Fig. 8 as the light-red band around the red curve (LL) that represents (129). For comparison the tree-level (LO) and linear-β 0 truncation (130) of the LL expression are displayed, but the difference between the latter approximation and the full result is hardly visible. The sizeable scale variation seen in the figure emphasizes the need for NLL resummation. The renormalized and endpoint-rearranged factorization formula derived in this work provides the starting point for this systematic improvement. However, the actual implementation requires the calculation of the presently unknown non-cusp parts of the O(α s ) anomalous dimension of J c (p 2 , ω, ω ) and S NLP (l + , l − , ω, ω ) for general ω ( ) , which require two-loop calculations beyond the scope of this work, as well as the O(α 2 s ) cusp parts of the anomalous dimensions of the NLP objects.

Conclusion
The lack of convergence of the convolution integrals appearing in subleading-power factorization theorems is currently the biggest obstacle to resummation of the powersuppressed large logarithmic corrections in collider physics. This applies in particular to the classic 1 → 2 and 2 → 1 processes, such as event shapes in e + e − annihilation, deep-inelastic scattering for x → 1, and the Drell-Yan threshold. In earlier work [18] on off-diagonal deep-inelastic scattering for x → 1, we showed that the hard function of the subleading power B-type operator factorizes in the limit when the collinear momentum fraction carried by one of the quark fields tends to zero. To avoid dealing with the non-perturbative parton distributions, in this work we considered the analogue of the off-diagonal channel for e + e − annihilation, "gluon thrust". Earlier results on the double-logarithmic resummation of this next-to-leading power kinematic configuration employed d-dimensional factorization and consistency relations [14,18], which do not readily generalize beyond the leading logarithms.
In this paper, we derived a novel endpoint factorization relation for the NLP thrust distribution in the two-jet region for the flavour-nonsinglet off-diagonal contribution, where a gluon-initiated jet recoils against a quark-antiquark pair, which involves subleading-power soft and jet functions defined at the cross-section level. The abovementioned factorization of the subleading power B-type operator plays again a crucial role for the consistency of the rearrangement that renders the expressions free from endpoint divergences. The framework developed in the present paper allows for the first time to remove systematically endpoint divergences in the convolution integrals of the SCET I factorization theorems and hence opens the path to systematic next-to-leadingpower resummation for collider observables involving soft quark emission. Employing the subtraction of endpoint divergences with the help of operatorial endpoint factorization conditions, we were able to reshuffle the factorization theorem such that the individual terms are free from endpoint divergences in convolutions and can be expressed in terms of renormalized hard, soft and collinear functions in four dimensions. At this point, standard renormalization-group techniques can be used to obtain the resummed integrands. The basic structure of the rearrangements turns out to be surprisingly similar to the one for Higgs decay through light-quark loops into two photons [21,22], which however is a SCET II process, where factorization applies to the amplitude. For the 1 → 2 and 2 → 1 processes considered here, there are no rapidity divergences releated to momentum modes of the same virtuality, and dimensional regularization suffices to make the convolutions well-defined [10,15], yet their divergence as → 0 requires rearrangement and refactorization as discussed here.
For the NLP thrust distribution in the power-suppressed two-jet region, where a gluon-initiated jet recoils against a quark-antiquark pair, we derived the necessary anomalous dimension of the NLP jet and soft functions using renormalization-group consistency and endpoint factorization relations. We also explicitly computed the one-loop anomalous dimension for the hard matching coefficients. These ingredients allowed us to perform the first resummation of the endpoint-divergent SCET I observables at the LL accuracy using exclusively renormalization-group methods. After resummation of the integrands, a judicious choice of initial conditions is necessary for the endpoint-subtraction terms, requiring a scale setting depending on the value of the convolution variable. We verified that our results simplified to double logarithmic accuracy agrees with the earlier results [14,18] and evaluated the numerical impact of the new LL corrections. Our main technical result for gluon thrust is (129), which provides a concise expression for the two-hemisphere invariant mass distribution in Laplace space.
The presented method relies on universal properties of the soft and collinear limits and can be applied to other SCET I problems. An immediate target for the further exploration of NLP resummation is the extension of the present work for gluon thrust and related soft-quark emission processes to NLL accuracy, which requires the computation of renormalization kernels for the NLP soft and jet functions at order O(α s ), i.e., at the two-loop level (counting phase-space integrals as loops). For the diagonal processes, the leading NLP logarithms have been considered for the thrust distribution and the DY threshold in [9,10,12,13], and they turn out to be almost trivial in comparison with the ones for the off-diagonal processes, as there are no endpoint divergences. This is no longer expected beyond the LL accuracy [15]. The present framework supplies a method to address the endpoint divergences that appear in the soft gluon limit at NLP, which is relevant to the diagonal processes. A C B1 and D B1 In this appendix we provide explicit results for the hard matching coefficient C B1 i (Q 2 , r) associated to the SCET operators defined in the second line of (11) and (14). In momentum space, where r andr = 1 − r denote the momentum fractions of the anti-collinear quark and anti-quark, respectively, and We then obtain the one-loop correction and one-loop evolution equation of the refactorization coefficient D B1 (p 2 ) by taking the r → 0 limit of the corresponding results for C B1 1 (Q 2 , r) and (55), Here endpoint-regular contributions diverge at most logarithmically in the limit r → 0.

A.1 One-loop results
Allowing for the one-loop correction to (34), we write We recall that charge conjugation implies C B1 2 (Q 2 , r) = −C B1 1 (Q 2 ,r), hence it suffices to discuss the case i = 1. A standard matching calculation of the γ * → qqg amplitude results in The limit r → 0 is of particular interest, and results in which shows that the B1 coefficient depends on two scales, Q 2 and rQ 2 . We can identify with the one-loop contribution to the leading-power hard matching coefficient coefficient Comparing to (55), we find For completeness, we check that the r → 1 limit of ∆ B1 (Q 2 , r) does not contain 1/r terms: A.2 Derivation of the evolution equation of the asymptotic refactorization coefficient D B1 Our goal is to show that the RGE (60) for the matching coefficient D B1 (rQ 2 ) can be derived from the standard evolution equation of the hard-matching coefficient C B1 i (Q 2 , r) for generic value of r. We note that this finding is distinct from the case of B1 operators containing a quark and a gluon in the same collinear direction, where the soft-gluon limit needs to be treated separately from the endpoint-regular piece [49].
The renormalization of SCET operators of the B1 type has been studied in detail in [44,45], except for the fermion-number zero case, to which the operators J B1,µ i (r) in (133) belong. The anomalous dimension matrix is determined by computing the MS renormalization factors Z µν ij (r, s) = δ ij g µν ⊥ δ(r − s) + δZ µν ij (r, s) defined via J B1,µ i,ren (r) = ds Z µν ij (r, s)J B1,ν j,bare (s) .
Following the lines of [44,45], we find at one-loop order δZ µν ij (r, s) = g µν ⊥ δZ ij (r, s) with where γ(r, s) = α s 2π The momentum fractions have support in the interval 0 ≤ r, s ≤ 1. Furthermore, T F = 1/2, and ∆ F = 1 (0) in the flavour-singlet (non-singlet) projection of the anticollinearχcχc fields in J B1,µ i . The expression in square brackets defines a symmetric distribution, which applies to integration over test functions of r or s, with definition The definition applies to any integration range, which is left unspecified here. In (146) the integrations are limited to 0 ≤ r, s ≤ 1, but we will also consider the case 0 ≤ r, s ≤ ∞ below.
The anomalous dimension can then be obtained as giving the RGEs Here we already used that the current renormalization for the specific operators J B1,µ i (r) considered here is diagonal in Lorentz indices, such that they can be dropped in the Wilson coefficients and anomalous dimension matrix.
We checked that convolving the tree-level coefficients C B1 i (Q 2 , r) with δZ ij (r, s) agrees with the divergent part of the explicit one-loop result for C B1 j (Q 2 , s) as given in the previous subsection for the flavour non-singlet projection, including both the contributions that are singular as well as those that are regular with respect to the dependence on s. This confirms the observation that the standard evolution equation for the B1 current can be used in this context, despite of the soft (anti-)quark singularity for r → 0 (r → 0). In the following we are interested in extracting from the evolution equation the part that describes the renormalization of the endpoint-singular piece captured by D B1 . The all-order identity C B1 2 (Q 2 , r) = −C B1 1 (Q 2 ,r) is consistent with the relations δZ 11 (r, s) = δZ 22 (r, s) = δZ 22 (r,s) and δZ 12 (r, s) = δZ 21 (r, s) = δZ 21 (r,s). Furthermore it can be checked that δZ 12 (r, s) does not contribute to the renormalization of D B1 (rQ 2 ), due to the independence on s (see below). Therefore it is sufficient to consider Z(r, s) ≡ δ(r − s) + δZ 11 (r, s) in the following.
The matching coefficient C B1 1 (Q 2 , r) for the problem at hand features a power-like endpoint singularity ∝ 1/r for r → 0, that may be modulated by logarithmic corrections, and is captured by the D B1 contribution in (55). In order to extract the renormalization factor for D B1 , we are interested in the convolution with Γ 11 (r, s) (or equivalently Z(r, s) at order α s ). The result of the convolution can be split into an endpoint divergent piece ∝ 1/s, and a regular part. The former describes the renormalization of D B1 , and the latter yields a mixing of D B1 into the endpoint-finite part of C B1 1 (Q 2 , r). Here we are interested only in the endpoint-singular contribution. In order to isolate it, we consider test functions that represent the most general possible dependence of the Wilson coefficient on the momentum fraction, and investigate the integral dr f (r) r Z(r, s) = F (s) s + endpoint-regular .
Here f (r) is an arbitrary test function defined for r ≥ 0, that depends logarithmically on r, and represents the D B1 coefficient. For example, f (r) could be given by some power of ln(r) or a more general polylogarithmic dependence. The right-hand side states that the result can be written as an endpoint-singular piece, with some function F (s) that also depends logarithmically on s, and an endpoint-regular contribution. Below we show that this is indeed the case, and that F (s) can be obtained via with a renormalization factor Z D B1 D B1 (r, s), and Z A0A0 the renormalization constant for the leading-power coefficient C A0 (Q 2 ), being given at order α s by The strategy for obtaining Z D B1 D B1 is to apply the method of regions [50] to (151), and expanding the integrand such that a homogeneous scaling in the limit r, s → 0 is obtained. As a first step, we need to check whether the integral over r exists, despite the endpoint-singular factor 1/r. This is trivially the case for the local contribution to Z(r, s), and we therefore only need to check the piece containing γ(r, s), dr f (r) r γ(r, s) = α s 2π The integrals in the second line exist due to the explicit factor r in the second line of (146), that cancels the endpoint-singularity 1/r. The first line has no endpoint singularity due to the lower integration boundary.
As the next step, we need to expand the integrand assuming a homogeneous scaling for r, s → 0. The first integral in the second line of (154) has already a homogeneous scaling, and only the first term in the round bracket contributes to F (s)/s, while the second term yields an endpoint-regular contribution s 0 drf (r)/s. In addition, the integral in the "singlet" term proportional to ∆ F is independent of s and can therefore not contribute to F (s) in (151).
Let us now turn to the first line of (154). Expanding in the region of small r, s allows us to replacer/s → 1. Next, we see that the integral 1 s drf (r)/r contributes only to the endpoint-finite part on the right-hand side of (151), and can therefore be dropped in the present discussion. However, the first term in the round bracket in the first line of (154) yields contributions that scale as 1/s. To isolate them, we replace the upper integration boundary 1 by a Heaviside function θ(1−r) within the integrand. The method of regions instructs us to homogenize also the Heaviside function, corresponding to extending the upper integration boundary to +∞. At this point, we use that the test function f (r) depends logarithmically on r. It can therefore be extended over the interval 0 ≤ r < ∞, in accordance with the method-of-region expansion. However, this procedure is too naive, since it implies that the integral diverges for r → ∞. This problem is linked to the observation that the respective integral yields contributions that are logarithmically enhanced as ln(s)/s for s → 0. They can in turn be isolated by the rewriting After this rewriting we can extend the integration range of the remaining integral to s ≤ r < ∞, yielding the desired homogeneous scaling. 16 Note that the relation above could equivalently be written at the level of plusdistributions, but we prefer to keep the explicit form here in order to be clear about the treatment of the integration domain in the various steps.
Altogether, this implies the contributions that scale as 1/s up to logarithmic corrections have been isolated, and are given by By comparing with (151) we can read off the renormalization coefficient of D B1 , whereγ (r, s) = α s 2π and the plus distributions are defined with respect to the interval 0 ≤ r, s < ∞ here. The anomalous dimension is, at order α s , given by which gives We can express the result equivalently in terms of ω andω using s = ω/Q, r =ω/Q. This yields d d ln µ D B1 (ω) = ∞ 0 dω γ D (ω, ω)D B1 (ω) , with γ D (ω, ω) = α s (C F − C A ) π δ(ω −ω) ln µ 2 −ωQ − iε This result is consistent with Eq. (3.2) in [36] obtained by a different method, when rewritingω → xω and ωQ → p 2 in terms of the variables x and p 2 used there, and accounting for the different convention of the order of arguments of γ D in (161).

B Alternative version of the endpoint-finite factorization formula
In this appendix consider the alternative split of the integration regions mentioned in Sec. 4.3, i.e. we present the factorization formula for the version (1), corresponding to ω and ω smaller than an endpoint factorization parameter Λ (integral I 1 ), and the complement region I 2 , see Fig. 6. Subtracting the complement region I 2 , the A-type term takes the from If we further assume that Λ 1/s R , 1/s L , then we can keep only the leading terms in 1/(s R Λ) and 1/(s L Λ) and the previous equation simplifies to − J c (s L , ω, ω ) S NLP (s R , s L , ω, ω ) + J c (s L , ω, ω ) S NLP (s R , s L , ω, ω ) .
The remaining part I 1 of the scaleless integral (70) must be combined withe the Btype term. To make notation more concise, we introduce r Λ ≡ Λ/Q. The B-type term for i = i = 1 is