\section{Introduction}
\label{sec:introduction}

Hadronic event shape variables in the two-jet region have played a key 
role for the development of resummation techniques in 
QCD \cite{Catani:1991kz,Catani:1992ua}.\footnote{See \cite{Dokshitzer:1978dr,Parisi:1979se,Collins:1981uk,Kodaira:1981np} for other early work discussing the break-down of the expansion in the 
strong coupling for back-to-back jets and semi-inclusive particle production.} 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$~\cite{Becher:2008cf}. 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\ll 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 \cite{Bonocore:2016awd,Moult:2018jjd,Beneke:2018gvs,Moult:2019mog,Bahjat-Abbas:2019fqa,Beneke:2019mua,Moult:2019uhz,Beneke:2019oqx,Moult:2019vou,Ajjath:2020ulr,Beneke:2020ibj,Ajjath:2020sjk,vanBeekveld:2021mxn}
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 
\begin{align}
T=\mbox{max}_{\vec n}\,\frac{\sum_{i}
\left|\vec{p_{i}}\cdot\vec{n}\right|}
{\sum_{i}\left|\vec{p_{i}}\right|},
\end{align}
where $i$ runs over all hadrons (partons) in the final state. 
The plane orthogonal to the thrust axis $\vec{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_L^2$ and $M_R^2$, respectively. As 
\begin{align}
\label{eq:tau}
\tau=1-T \to 0,
\end{align}
the particles cluster into a pair of back-to-back jets, 
and perturbation theory in the strong coupling $\alpha_s(Q)$ 
breaks down due to large logarithms $\ln\tau$ 
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
\begin{equation}
e^+ e^-\to\gamma^*\to [g]_{c} + [q\bar q]_{\bar c}\,,
\end{equation}
where the direction of gluon jet is defined to be the 
``collinear direction'' and the direction of the recoiling 
$q\bar q$-jet the ``anti-collinear'' one. This process 
starts at order $\alpha_s$ only (see Fig.~\ref{fig:QCD}). 
More importantly, as $\tau\to 0$, it does not have the 
leading-power $[\ln\tau/\tau]_+$ behaviour from
soft-gluon emission. Instead the leading term is 
$\alpha_s\ln \tau$, and the entire process is of 
next-to-leading power in $\tau$,\footnote{For clarification, 
we note that we do not consider here non-perturbative 
power corrections in the strong 
interaction scale $\Lambda_{\rm QCD}$, which are  
expected to be  $\mathcal{O}(\Lambda_{\rm QCD}/(\tau Q))$, that is, 
we assume the strong hierarchy $\Lambda_{\rm QCD}/Q\ll \tau\ll1$.} 
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 \cite{Moult:2019uhz} and an 
expression for the leading double logarithms was written down, 
which was subsequently derived from $d$-dimensional consistency 
relations \cite{Beneke:2020ibj}. The result 
of \cite{Moult:2019uhz} exhibits an unconventional ``quark'' 
Sudakov form factor, related to an enhancement when 
either the quark or anti-quark in the 
$q\bar q$-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 
\cite{Moult:2019uhz,Moult:2019mog,Beneke:2020ibj}, 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 \cite{Beneke:2020ibj}, and Higgs decay to two photons 
through bottom-quark loops \cite{Liu:2019oav,Liu:2020wbn} 
to derive a factorization theorem suitable for systematic 
resummation. To this 
end we combine standard SCET factorization with endpoint 
factorization.

\begin{figure}
\begin{center}
\includegraphics[width=0.22\textwidth]{figs/QCD1}\qquad\qquad
\includegraphics[width=0.22\textwidth]{figs/QCD2}
\caption{\label{fig:QCD} Leading $\mathcal{O}(\alpha_s)$ QCD 
diagrams contributing to gluon thrust.}
\end{center}
\end{figure}

The outline of this paper is as follows: In Sec.~\ref{sec:heuristic} 
we provide a heuristic discussion of the factorization formula, 
which motivates the proposed endpoint subtraction. In 
Sec.~\ref{sec:bare} 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.~\ref{sec:refactorization}, 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.~\ref{sec:consistency} and 
solved with leading-logarithmic (LL) accuracy. Finally, in 
Sec.~\ref{sec:resummation}, 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\epsilon$ space-time dimensions, and the 
$\overline{\rm MS}$ subtraction 
scheme whenever renormalization is required.


\section{Heuristic discussion}
\label{sec:heuristic}

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 $\alpha_s$ (see Fig.~\ref{fig:QCD}) relative to the LO 
process, but also by one power of $\tau$.

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 $\mathcal{O}(\alpha_s)$:
\begin{enumerate}
\item[I] Quark and anti-quark have both large  anti-collinear 
momentum and 
form a single jet, which recoils against the gluon.
\item[II] The quark or anti-quark is anti-collinear and balances
the gluon momentum, while the other of the two is soft.
\end{enumerate}
Both configurations are suppressed by a factor of $\tau$. 
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 $\alpha_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.~\ref{fig:QCD}. 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_a^2\sim Q^2$. The $q_a$-internal 
line is integrated out, and the primary hard SCET vertex produces 
a $q\bar{q} g$ state directly, see the first diagram in 
Fig.~\ref{fig:SCET}. Since there are two partons ($q\bar q$) 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 anti-quark becomes soft) 
is no longer hard. Hence the primary vertex is the standard 
$\gamma^*\to q\bar {q}$ 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 $\mathcal{L}^{(1)}_{\xi q}$ describing 
soft (anti-) quark emission, see the second and third diagram in  
Fig.~\ref{fig:SCET}.  

\begin{figure}
\begin{center}
\includegraphics[width=0.28\textwidth]{figs/B1.pdf}\hspace{0.15cm}
\includegraphics[width=0.28\textwidth]{figs/A11.pdf}\hspace{0.15cm}
\includegraphics[width=0.28\textwidth]{figs/A12.pdf}
\caption{\label{fig:SCET} SCET representation of the 
``off-diagonal'' gluon-thrust amplitude in the two-jet region. 
The diagrams show the NLP structure 
of the process with all LP interactions removed.}
\end{center}
\end{figure}

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:
\begin{align}
\frac{1}{\sigma_{0}}\frac{\widetilde{d\sigma}}{ds_R ds_L}
&=
\int_{0}^{\Lambda} d\omega d\omega'\,
\left|C^{A0}\right |^{2}\times
\mathcal{J}^{(\bar{q})}_{\bar{c}}\times 
\mathcal{J}_{c}\left(\omega,\omega'\right)
\otimes S_{\rm NLP} \left(\omega,\omega'\right)
\nonumber \\&+
\int_{\Lambda/Q}^{1-\Lambda/Q} drdr'\,
C^{B1}(r) C^{B1}(r')^{*}\otimes
\mathcal{J}^{q\bar q}_{\bar{c}}\left(r,r'\right)\times 
\mathcal{J}_{c}^{(g)}\times S^{(g)}\,.
\label{eq:ffschematic}
\end{align}
Here $C$ are hard matching coefficients, $\mathcal{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 $\omega$, $\omega'$ 
convolution integrals are logarithmically divergent for 
$\omega,\omega'\to \infty$, 
while the hard-anti-collinear integrals over $r$, $r'$ diverge 
logarithmically when $r,r'\to 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_{\rm NLP}$ that appears in possibility II has {\em large} soft 
momentum $\omega$ and can equally well be regarded as part of 
the anti-collinear function $\mathcal{J}^{q\bar q}_{\bar{c}}$ 
in the second line with {\em small} 
anti-collinear momentum fraction $r$. Removing the  
quark from $S_{\rm NLP}$ leaves the single-particle soft 
function $S^{(g)}$, hence $S_{\rm NLP}\to S^{(g)}$ and 
$\mathcal{J}^{(\bar{q})}_{\bar{c}} \to 
\mathcal{J}^{q\bar q}_{\bar{c}}$. At the same time, the hard 
process changes from A0-type to B1-type. 
Thus, the {\em integrands} of the two terms in 
\eqref{eq:ffschematic} should become identical in the 
respective limits. A rearrangement at the integrand level 
can then be performed in the singular limits, which 
introduces the factorization parameter $\Lambda$, such that 
both integrals are separately finite. At this point 
one can remove the dimensional regulator and use standard 
renormalization-group techniques to sum the 
logarithms in the hard, jet and soft functions. These 
involve only LP interactions. NLP interactions occur 
only once (at amplitude level), leading to the 
two terms in \eqref{eq:ffschematic}, and therefore provide 
only one additional logarithmic, but finite integration from the 
convolutions.

 
\section{Bare factorization theorem}
\label{sec:bare}

The distribution of an event shape $e$ in $e^+ e^-$ annihilation 
into hadrons through a photon with virtuality $Q^2=q^2$ is
\begin{align}
\frac{d\sigma}{de} =
\sum_{X}\int \!d\mbox{PS}_{X}\left(2\pi\right)^{d}\delta^{(d)}\left(q-p_{X}\right)L_{\mu\nu}\left(q\right)
\delta(e-\hat{e}(X))\left\langle 0|J^{\dagger\nu}(0)|X\right\rangle \left\langle X|J^{\mu}(0)|0\right\rangle ,
\label{eq:dsigma}
\end{align}
with the leptonic tensor
\begin{align}
L_{\mu\nu}\left(q\right)=\frac{8\pi^{2}\alpha_{\rm em}^{2}}
{3Q^{4}}\left(-g_{\mu\nu}+\frac{q_{\mu}q_{\nu}}{Q^{2}}\right)\,
\label{eq:leptontensor}
\end{align}
and the electromagnetic current $J^\mu=\bar\psi\gamma^\mu\psi$, 
for simplicity of a single quark flavour with electric charge 
$e_\psi=+1$.
For a given final state $X$, $\hat{e}(X)$ returns the value of the 
event shape variable. For small $\tau$, 
\begin{equation}
\tau = \frac{M_R^2+M_L^2}{Q^2}+\mathcal{O}(\tau^2),
\label{eq:tauMRML}
\end{equation}
and the thrust distribution is related to the two-hemisphere 
invariant-mass distribution by
\begin{equation}
\frac{d\sigma}{d\tau} =
\int dM_R^2dM_L^2 \,\delta\!\left(\tau-\frac{M_R^2+M_L^2}{Q^2}
\right)\frac{d\sigma}{dM_R^2 dM_L^2}\,.
\label{eq:taudist}
\end{equation}

Eq.~(\ref{eq:dsigma}) 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 \cite{Catani:1991kz,Catani:1992ua}. The resummation of NLP 
double-logarithmic corrections to the quark-antiquark 
back-to-back jet configuration, which is present already at LP,
 has been derived recently 
\cite{Moult:2018jjd}. Defining the left hemisphere to be the 
collinear one and employing the SCET framework, the hemisphere 
invariant-mass distribution for $M_R^2,M_L^2\ll Q^2$ is expressed 
at LP as 
\begin{equation}
\frac{1}{\sigma_0}\frac{d\sigma}{dM_R^2 dM_L^2} = 
|C^{\rm A0}(Q^2)|^2 \int dl_+ dl_-\,
\mathcal{J}_c^{(q)}(M_L^2-Q l_-)\mathcal{J}_{\bar c}^{(\bar q)}(M_R^2-Q l_+)\,
S_{\rm LP}(l_+,l_-)
\end{equation}
in terms of the convolution of the jet (anti-) collinear functions 
with the hemisphere-soft function, multiplied by a 
universal hard function\cite{Fleming:2007qr,Schwartz:2007ib,Bauer:2008dt,Becher:2008cf}.\footnote{By charge conjugation symmetry, the 
quark and anti-quark jet functions coincide, $\mathcal{J}_{\bar c}^{(\bar q)}(p^2) = \mathcal{J}_c^{(q)}(p^2)$.} The distribution is normalized to 
the leading-order total cross section
\begin{equation}
\sigma_0=\frac{4\pi N_c\alpha_{\rm em}^{2}}
{3Q^2}\,,
\end{equation} 
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 \cite{Beneke:2019oqx}.

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 \cite{Moult:2019uhz}. We introduce two light-like vectors 
$n_\pm^\mu$, satisfying $\npn_-=2$ and pointing into the 
directions of the back-to-back jets. The vector $n_-^\mu$ defines 
the direction of the collinear modes of the gluon jet. (Anti-) 
collinear modes therefore have large components $n_+ p_c\sim Q$ 
($n_- p_{\bar c}\sim Q$). 
The transverse momenta in the jet are of order 
$p_\perp \sim \lambda Q\sim\sqrt{\tau} Q \ll Q$, where 
$\lambda\sim\sqrt{\tau}$ denotes the SCET power-counting 
parameter. The momenta of soft particles scale uniformly, 
$p_s\sim \tau Q$. In the adopted reference frame, $q_\perp^\mu=0.$
We further note that since gluon thrust vanishes at LP, there 
are no kinematic NLP corrections, and \eqref{eq:tauMRML}, 
\eqref{eq:taudist} continue to hold.

After integrating out the hard modes, the electromagnetic current 
matches to 
\begin{eqnarray}
\bar\psi\gamma_\perp^\mu\psi(0) &=& \int dtd\bar{t}\,
\widetilde{C}^{\rm A0}(t,\bar{t})\,
\bar{\chi}_c(tn_+)\gamma_\perp^\mu\chi_{\bar c}(\bar tn_-) 
+ (c\leftrightarrow \bar{c})
\nonumber\\
&&+\,\sum_{i=1,2}\int dt d\bar{t}_1d\bar{t}_2\,
\widetilde{C}^{\rm B1}_i(t,\bar{t}_1,\bar{t}_2)\,
\bar{\chi}_{\bar c}(\bar{t}_1n_-)\Gamma_i^{\mu\nu}\mathcal{A}_{c\perp\nu}
(tn_+)\chi_{\bar c}(\bar{t}_2n_-)
+\ldots \qquad
\label{eq:hardmatching}
\end{eqnarray}
in SCET. Here $\mathcal{A}_{c\perp\nu}$ denotes the gauge-invariant 
gluon field operator dressed with Wilson lines, which in light-cone gauge is related to the gluon field by $\mathcal{A}_{c\perp\nu}=
g_s A_{c\perp\nu}$. The omitted terms are higher than NLP 
corrections and purely gluonic ``flavour-singlet'' operators. 
For simplicity, we neglect the mixing 
of $q\bar{q}$ 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 \eqref{eq:hardmatching} 
represents the $\mathcal{O}(\lambda)$ suppressed 
``B-type'' SCET operator, which directly produces the 
configuration relevant to gluon thrust, a collinear gluon and an 
anti-collinear quark-antiquark pair. The first line can 
nevertheless contribute to gluon thrust through a time-ordered 
product with the  $\mathcal{O}(\lambda)$ suppressed SCET 
interaction  \cite{Beneke:2002ph,Beneke:2002ni}
\begin{align}
\mathcal{L}_{\xi q}(x)=\bar{q}_s(x_{-})\slashed{\mathcal{A}}_{c\perp}(x)\chi_{c}(x)+\rm{h.c.}\,,
\label{eq:Lxiq}
\end{align} 
which converts a collinear quark into a collinear gluon and 
a soft quark. Here
\begin{equation}
x_\mp^\mu = (n_\pm\cdot x)\frac{n_\mp^\mu}{2}\,,
\label{eq:xminusdef}
\end{equation}
and we further define the scalar quantity 
$x_\mp \equiv n_\pm\cdot 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.\footnote{The hermitian conjugate in \eqref{eq:Lxiq} provides 
a second contribution, in which the antiquark turns into the 
collinear gluon and a soft antiquark.} The Dirac structures in the 
power-suppressed B-type operator are given by
\begin{equation}
\Gamma_1^{\mu\nu}=\frac{\slashed n_{-}}{2}\gamma_{\perp}^{\nu}\gamma_{\perp}^{\mu}\,,
\qquad\quad
\Gamma_2^{\mu\nu}=\frac{\slashed n_{-}}{2}\gamma_{\perp}^{\mu}\gamma_{\perp}^{\nu}\,.
\label{eq:BtypeDirac}
\end{equation}
We note that boost invariance together with the projection 
property $\slashed n_{-} \chi_c=0$, 
$\slashed n_+ \chi_{\bar c}=0$ implies that only the transverse 
components of the vector current contribute.

Eq.~\eqref{eq:hardmatching} 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 \cite{Bauer:2001yt} $Y_{n_-}(x)=\mathcal{P} \exp\left[i g_s \int_0^\infty ds n_-A_s(x+sn_-) \right]$ 
($Y_{n_+}(x)$ for anti-collinear fields), the collinear, 
anti-collinear and soft fields are decoupled at LP. The final 
state  $|X\rangle = |X_c\rangle|X_{\bar c}\rangle 
|X_s\rangle$ 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 
\eqref{eq:hardmatching}, more precisely the time-ordered 
product with \eqref{eq:Lxiq} 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.

\subsection{A-type term (soft quark)}

This term stems from the square of the matrix element
\begin{equation}
\langle X_c|\langle X_{\bar c}|\langle X_s| \int d^d x \,T\,[\bar{\chi}_c(tn_+)\gamma_\perp^\mu\chi_{\bar c}
(\bar tn_-), i\mathcal{L}_{\xi q}(x)]\,|0\rangle\,,
\label{eq:atypetproduct}
\end{equation}
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 \eqref{eq:atypetproduct} 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 \cite{Bauer:2001yt}
\begin{eqnarray}
\chi_c(x) \to Y_{n_-}(x_-)\chi_c(x),&\qquad& 
\chi_{\bar c}(x) \to Y_{n_+}(x_+)\chi_{\bar c}(x), 
\\[0.1cm]
&&\hspace*{-4cm}
\mathcal{A}_{c\perp}^\mu(x)\to Y_{n_-}(x_-)\mathcal{A}_{c\perp}^\mu(x)
Y^\dagger_{n_-}(x_-)\,.
\end{eqnarray}
Inspection of \eqref{eq:hardmatching}, 
\eqref{eq:atypetproduct} shows that the following hard, (anti-) 
collinear and soft functions appear.

\paragraph{Hard function} This term involves the hard function 
$\widetilde{C}^{\rm A0}(t,\bar{t})$ of 
the LP A-type current. We define the momentum-space 
coefficient
\begin{equation}
C^{\rm A0}(n_+p_c, n_-p_{\bar{c}}) = 
\int dtd\bar{t}\,e^{it n_+ p_c+i \bar{t}n_- p_{\bar c}}\,
\widetilde{C}^{\rm A0}(t,\bar{t})\,.
\label{eq:CA0momentumspace}
\end{equation}
It depends only on the product $Q^2 = n_+p_c n_-p_{\bar{c}}$  of 
the large components of the collinear and anti-collinear momenta. 
At LO, we find $C^{\rm A0}(Q^2)=1+\mathcal{O}(\alpha_s)$. 
The two-loop result can be found in \cite{Becher:2007ty}.

\paragraph{Anti-collinear function}

The $\mathcal{L}_{\xi 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 
\begin{eqnarray}
&&\frac{1}{2\pi}\sum_{X_{\bar{c}}}\int \!d\mbox{PS}_{X_{\bar{c}}}
\langle 0|\bar{\chi}_{\bar{c}}(x)_{b\beta}|X_{\bar{c}}\rangle 
\langle X_{\bar{c}}|\chi_{\bar{c}}(0)_{a\alpha}|0\rangle 
\nonumber\\
&& \equiv\,\delta_{ab}\int\!\frac{d^d p}{(2\pi)^d}\,n_{-}p\,e^{-ipx}
\mathcal{J}^{(\bar{q})}_{\bar{c}}(p^2)\left(\frac{\slashed n_{+}}{2}\right)_{\alpha\beta}.
\label{eq:LPjet}
\end{eqnarray}
At LO, we have $\mathcal{J}^{(\bar{q})}_{\bar{c}}(p^2)=\delta^{+}(p^2)
\equiv \theta(p^0)\delta(p^2)$. This object is well understood and computed up to the third order in $\alpha_s$ \cite{Bruser:2018rad}. Here and below we use small Latin 
letters for colour and Greek letters for Dirac spinor indices. 

\paragraph{Collinear function} This function picks up the collinear fields from the 
$\mathcal{L}_{\xi q}$ insertion and represents a new non-local 
jet function that appears at NLP. Defining the non-local operator
\begin{equation}
\mathcal{O}_{a\alpha;b\beta}(\omega,x) = 
\int d^dy \,e^{i y_-\omega}\,T\big\{\bar{\chi}_{c,b\beta}(x),
\left[\,\not\hspace*{-1.5mm}\mathcal{A}_{\perp c}\chi_c\right]_{a\alpha}(x+y)\big\}\,,
\label{eq:nonlocaljetfield}
\end{equation}
where $y_-=n_+\cdot y/2$, its jet function defined as in the LP 
case \eqref{eq:LPjet} by 
\begin{eqnarray}
&&\frac{1}{2\pi}\sum_{X_{c}}\int d\mbox{PS}_{X_{c}}\,
\frac{1}{g_{s}^{2}}\left\langle 0\right|
\mathcal{O}^\dagger_{b'\beta';a'\alpha'}(\omega^\prime,x)
\left|X_{c}\right\rangle 
\,\left[\frac{\slashed n_{+}}{2}\right]_{\alpha' \alpha} \!\!
\left\langle X_{c}\right|\mathcal{O}_{a\alpha;b\beta}(\omega,0)\left|0\right\rangle 	
\nonumber\\&& \hspace*{0.5cm} = \,(d-2)\,
\left[\frac{\slashed n_{-}}{2}\right]_{\beta'\beta}\,
\int\frac{d^dp}{(2\pi)^d}\,e^{-ipx}\,\Big\{
\left[t^{A}\right]_{ab}\left[t^{A}\right]_{b'a'}
\,\mathcal{J}_{c}(p^2,\omega,\omega^\prime)
\nonumber\\
&& \hspace*{1cm}+
\left[t^{A}\right]_{aa'}\left[t^{A}\right]_{b'b}
\,\widehat{\mathcal{J}}_{c}(p^2,\omega,\omega^\prime)
\Big\}\,.
\label{eq:defjetNLP}
\end{eqnarray}
The factor $1/g_s^2$ has been introduced for convenience, such 
that the leading term is $\mathcal{O}(\alpha_s^0)$. The Dirac matrix 
$\left[\slashed n_{+}/2\right]_{\alpha' \alpha}$ 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 \eqref{eq:NLPsoftfndef} below.\footnote{This 
would remain true even if we 
had allowed for flavour-singlet contractions of the quark fields. 
These do not exist for the A-type term, since the contraction of 
different modes, that is the collinear quark field with the 
soft quark field from the 
soft function, is not allowed in SCET.} At 
LO in $\alpha_s$, the collinear final state $X_c$ 
consists of a single gluon,\footnote{The $q\bar{q} g$ final state 
at this order corresponds to a disconnected diagram, which does 
not contribute to the scattering amplitude.} and only the 
unhatted jet function is non-vanishing:
\begin{align}
\mathcal{J}_{c}(p^2,\omega,\omega')&= 
 \frac{1}{\omega\omega'}\,\delta^+(p^{2})\,,\\
\widehat{\mathcal{J}}_{c}(p^2,\omega,\omega')&=0\,.
\label{eq:AtypeJctree}
\end{align}

\paragraph{Soft function} 

Collecting the soft Wilson lines from 
the soft-decoupling field redefinition, as well as the 
soft quark field from the Lagrangian insertion \eqref{eq:Lxiq} 
leads to the following definition of the NLP soft quark function:
\begin{eqnarray}
&&g_s^2
\int\frac{dx_-}{2\pi}\frac{dx_-^\prime}{2\pi}
\,e^{-i (x_-\omega-x_-^\prime\omega^\prime)}	
\left\langle 0\right|\overline{T}\left\{
\left[Y_{n_+}^{\dagger}(0)Y_{n_-}(0)\right]_{cb'}
\left[Y^\dagger_{n_-}q_s\right]_{\alpha'a'}(x_-^\prime)
\right\}
\nonumber \\ 
&&\hspace*{1cm}\times\,
\mathcal{P}_s(l_+,l_-)\,
T\left\{\left[\bar{q}_sY_{n_-}\right]_{\alpha a}(x_-)\left[Y_{n_-}^{\dagger}\left(0\right)Y_{n_+}
\left(0\right)\right]_{bc}\right\}
\left|0\right\rangle 	\nonumber\\
&& = \, \left(\frac{\slashed n_{+}}{2}\right)_{\alpha'\alpha}\!
\left\{
\delta_{a'a}\delta_{bb'} \,
S_{\text{NLP}}(l_+,l_-,\omega,\omega^\prime) +
\delta_{ba}\delta_{a'b'} \,
\widehat{S}_{\text{NLP}}(l_+,l_-,\omega,\omega^\prime)
\right\} +\ldots\,,\qquad\quad
\label{eq:NLPsoftfndef}
\end{eqnarray}
where $x_-$ is defined below \eqref{eq:xminusdef}, and 
\begin{equation}
\mathcal{P}_s(l_+,l_-) \equiv \sum_{X_s}\int d\mbox{PS}_{X_s}\,
\delta(l_+-n_+ p_{X_s}^R)\delta(l_--n_- p_{X_s}^L)\,|X_s\rangle
\langle X_s|\,,
\label{eq:softmeasurement}
\end{equation}
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 anti-collinear side. 
The dots in \eqref{eq:NLPsoftfndef} denote terms proportional to 
$\slashed n_{-}$, which do not contribute due to the projection 
property of collinear fields in other parts of the process.
The pre\-factor is chosen for convenience, such that 
the soft function starts at $\mathcal{O}(\alpha_s)$. 
The soft function for thrust is obtained from
\begin{equation}
S_{\text{NLP}}^{T}(k,\omega,\omega') = \int_0^\infty 
dl_+d l_-\,\delta(k-l_+ -l_-) \, 
S_{\text{NLP}}(l_+,l_-,\omega,\omega')\,.
\end{equation}
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 $\mathcal{O}(\alpha_s)$, 
$\widehat{S}_{\text{NLP}}(l_+,l_-,\omega,\omega^\prime)=0$, 
while
\begin{eqnarray}
S_{\text{NLP}}(l_+,l_-,\omega,\omega^\prime) &=&
\frac{\alpha_s}{4\pi}\frac{\theta(\omega)\delta(\omega-\omega')}
{e^{-\epsilon\gamma_E}\Gamma(1-\epsilon)}\,
\bigg[\delta(l_-)\theta(\omega-l_+)\theta(l_+)\,
\omega\left(\frac{l_+\omega}{\mu^2}\right)^{\!-\epsilon}
\nonumber\\
&& -\,\frac{\omega^2}{1-\epsilon}\,\delta(l_+)\delta(l_--\omega)
\left(\frac{\omega^2}{\mu^2}\right)^{\!-\epsilon}\,\bigg]\,.
\label{eq:NLPsofttree}
\end{eqnarray}
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 \cite{Beneke:2018gvs}, which is already 
known to $\mathcal{O}(\alpha_s^2)$ 
\cite{Beneke:2018gvs,Broggio:2021fnr}.

\paragraph{Factorization formula}
With these definitions, we obtain the soft-quark contribution 
to gluon thrust in the two-jet region in the form
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{d\sigma}{dM_R^2 dM_L^2}|_{\rm A-type} &=& 
\frac{2 C_F}{Q}\,f(\epsilon)\,
|C^{\rm A0}(Q^2)|^2 \int_0^\infty dl_+ dl_-\,
\int d\omega d\omega'\,
\mathcal{J}_{\bar c}^{(\bar q)}(M_R^2-Q l_+)
\nonumber\\
&&\hspace*{-2cm}\times\,
\bigg\{\,\mathcal{J}_{c}(M_{L}^{2}-Ql_-,\omega,\omega')
\,S_{\rm NLP}(l_+,l_-,\omega,\omega')
\nonumber\\
&&\hspace*{-1.2cm}
+\,\widehat{\mathcal{J}}_{c}(M_{L}^{2}-Ql_-,\omega,\omega')\,
\widehat{S}_{\rm NLP}(l_+,l_-,\omega,\omega')\,\bigg\}\,,
\label{eq:Atypefactformula}
\end{eqnarray}
where 
\begin{eqnarray}
f(\epsilon) = \bigg( \frac{Q^{2}}{4\pi} \bigg)^{-\epsilon} \frac{ (1-\epsilon)^2 \Gamma(1-\epsilon)}{ \Gamma(2- 2\epsilon)}
\end{eqnarray}
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 
\begin{equation}
\frac{\widetilde{d\sigma}}{ds_R ds_L} = 
\int_0^\infty dM_R^2 dM_L^2\,\,e^{-s_R M_R^2/Q}\,
e^{-s_L M_L^2/Q}\,\frac{d\sigma}{dM_R^2 dM_L^2}
\label{eq:laplacetrafo}
\end{equation}
of the two-hemisphere invariant-mass distribution. In general, 
the Laplace transform of collinear functions with variable 
$p^2$ are defined with $\int_0^\infty dp^2\,e^{-s p^2/Q}$, while 
for soft functions we use $\int_0^\infty dl_+ dl_-\,
e^{-s_R l_+-s_L l_-}$. Then
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{\rm A-type} &=& 
\frac{2 C_F}{Q} \,f(\epsilon)\,|C^{\rm A0}(Q^2)|^2 \,
\widetilde{\mathcal{J}}_{\bar c}^{(\bar q)}(s_R)
\,\int d\omega d\omega'\,
\nonumber\\
&&\hspace*{-3cm}\times\,
\bigg\{\,\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
\,\widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')
+\;\,\widetilde{\!\!\widehat{\mathcal{J}}}_{\!c}(s_{L},
\omega,\omega')
\,\,\widetilde{\!\widehat{S}}_{\rm NLP}(s_R,s_L,\omega,\omega')\,\bigg\}\,. \quad
\label{eq:AtypeLaplace}
\end{eqnarray}
The total A-type term is twice the above \eqref{eq:Atypefactformula}, 
\eqref{eq:AtypeLaplace}, since there is 
an equal contribution from the soft anti-quark phase-space 
region. 

\subsection{B-type term}

The B-type term is obtained under the assumption that the 
intermediate propagator in Fig.~\ref{fig:QCD} is hard, which 
leads to the B1 operator in the matching equation 
\eqref{eq:hardmatching}. We then need the square of   
the matrix elements
\begin{equation}
\langle X_c|\langle X_{\bar c}|\langle X_s|\,\bar{\chi}_{\bar c}(\bar{t}_1n_-)\Gamma_i^{\mu\nu}
\mathcal{A}_{c\perp\nu}(tn_+)\chi_{\bar c}
(\bar{t}_2n_-) \, |0\rangle\,,
\label{eq:btypeME}
\end{equation}
summed and integrated over all possible final states.
Inspection of \eqref{eq:hardmatching}, 
\eqref{eq:btypeME} shows that the following hard, (anti-) 
collinear and soft functions appear.

\paragraph{Hard function} 
There are two relevant operators, $i=1,2$, 
which differ by their Dirac structures \eqref{eq:BtypeDirac}. We 
introduce the Fourier transforms
\begin{equation}
C^{\rm B1}_i(n_+p_c, r n_-p_{\bar{c}},\bar{r} n_-p_{\bar{c}}) = 
n_- p_{\bar c} \int dtd\bar{t}_1d\bar{t}_2\,e^{it n_+ p_c
+i (\bar{t}_1 r + \bar{t}_2 \bar{r})n_- p_{\bar c}}\,
\widetilde{C}^{\rm B1}_i(t,\bar{t}_1,\bar{t}_2)
\label{eq:CB1momentumspace}
\end{equation}
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_-p_{\bar{c}}$  of 
the large components of the collinear and total 
anti-collinear momentum, and $r$, which denotes the fraction 
of total anti-collinear momentum carried by the quark. 
Accordingly, $\bar{r}=1-r$ is the momentum fraction of the 
anti-quark. At the leading order in $\alpha_s$, we find\footnote{We 
recall that the definition of $\mathcal{A}_{c\perp}^\mu$ absorbs 
the factor $g_s$ from the tree-level diagram.}
\begin{align}
C_1^{\rm B1}(Q^2,r) &=\frac{1}{r}\,,
\label{eq:CB11tree}\\
C_2^{\rm B1}(Q^2,r) &=-\frac{1}{\bar{r}}\,.
\end{align}
From the transformation of \eqref{eq:hardmatching} under charge 
conjugation, we find that charge conjugation invariance implies that 
$C_2^{\rm B1}(Q^2,r) = -C_1^{\rm B1}(Q^2,\bar{r})$ to 
all orders in the strong coupling. 
We notice that the matching coefficient $C_1^{\rm B1}$ ($C_2^{\rm B1}$) 
is singular in the limit $r\to0$ ($r\to1$), 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 
\eqref{eq:hardmatching} into a left- and right-handed piece 
by inserting $1=P_L+P_R$. In four space-time 
dimensions we have\footnote{Convention: $\epsilon^{0123} = +1$}
\begin{eqnarray}
\label{eq:b1helicity1}
\Gamma^{\mu\nu}_1 P_L &=& g^{\mu\nu}_{L} \slashed n_{-} P_L\,,
\\
\Gamma^{\mu\nu}_2 P_L &=& g^{\mu\nu}_{R} \slashed n_{-} P_L\,,
\\
\Gamma^{\mu\nu}_1 P_R &=& g^{\mu\nu}_{R} \slashed n_{-} P_R\,,
\\
\Gamma^{\mu\nu}_2 P_R &=& g^{\mu\nu}_{L} \slashed n_{-} P_R\,,
\label{eq:b1helicity2}
\end{eqnarray}
where 
\begin{equation}
g_{L/R}^{\mu\nu} = \frac{1}{2}\left(g_\perp^{\mu\nu} 
\pm\frac{i}{2}\epsilon^{\rho\sigma\mu\nu} n_{-\rho} n_{+\sigma}\right)
\label{eq:gLR}
\end{equation}
projects on the left $(-1)$ / right $(+1)$ 
helicity of the collinear gluon in the 
B1 operator \eqref{eq:hardmatching}. 
Since helicity is conserved and amplitudes with 
gluons of different helicity do not interfer (see 
\eqref{eq:LPgluonjet} below), the virtual photon and gluon 
always have the same helicity. Focusing on the negative  
helicity case for definiteness, we conclude from 
\eqref{eq:b1helicity1}~--~\eqref{eq:b1helicity2} that 
the operator with Dirac structure $\Gamma^{\mu\nu}_1$ 
produces a left-handed out-going quark, while for 
$\Gamma^{\mu\nu}_2$ it is right-handed. It follows that 
the two operators cannot interfer in the non-singlet channel, 
where the $\bar{\chi}_{\bar c}(\bar{t}_1n_-)
\chi_{\bar c}(\bar{t}_2n_-)$ fields in the operator in 
\eqref{eq:hardmatching} 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 $\epsilon\to 0$. 

The limit $r\to 0$ in the B1 matching coefficient implies 
that the momentum of the out-going quark in the left 
diagram of Fig.~\ref{fig:QCD} goes to zero, and the 
intermediate quark propagator goes on-shell. The helicity 
argument above implies that only $C_1^{\rm B1}(Q^2,r)$ 
corresponding to $\Gamma^{\mu\nu}_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 $\gamma^* \to 
q+\bar{q}$ and $q\to q+g$. As $r\to 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  $\Gamma^{\mu\nu}_1$. The situation is reversed 
in the soft anti-quark limit 
$r\to 1$ and implies right-handedness, which 
singles out  $\Gamma^{\mu\nu}_2$. We conclude 
that, to all orders in perturbation theory, only 
$C_1^{\rm B1}(Q^2,r)$ ($C_2^{\rm B1}(Q^2,r)$) can be  
singular in the endpoint $r\to 0$ ($r\to 1$). 
  
\paragraph{Collinear function}

The collinear function is simply the standard LP gluon jet 
function, defined as  
\begin{eqnarray}
&&\frac{1}{2\pi}\frac{1}{g_{s}^2}
\sum_{X_c}\int \!d\mbox{PS}_{X_c}\,
\langle 0|\mathcal{A}_{c\perp\mu}^B(x)|X_{c}\rangle 
\langle X_c|\mathcal{A}_{c\perp\nu}^C(0)|0\rangle 
\nonumber\\
&& \equiv\,\delta^{BC}\,(-g_{\mu\nu}^\perp)
\int\!\frac{d^{d}p}{(2\pi)^d}\,e^{-ipx}
\mathcal{J}^{(g)}_{c}(p^2).
\label{eq:LPgluonjet}
\end{eqnarray}
At LO, we have $\mathcal{J}^{(g)}_{c}(p^2)=
\delta^+(p^2)$. This object is well understood and known up 
to the third order in $\alpha_s$ \cite{Banerjee:2018ozf}. 

\paragraph{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
\begin{equation}
\mathcal{Q}_i^{A\mu\nu}(x,r) = 
\frac{1}{2\pi}\int_0^\infty \!\!d\bar{t}\,
e^{-i r \bar{t} n_-\cdot p_{\bar{c}}}\,
\bar{\chi}_{\bar{c}}(x+\bar{t} n_-) t^A\,\Gamma_i^{\mu\nu}
\chi_{\bar{c}}(x)\,.
\end{equation}
At leading order, $r$ can be interpreted as the momentum fraction 
of the out-going quark. The Dirac structures refer to 
\eqref{eq:BtypeDirac}. 
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 
\begin{eqnarray}
&&\frac{g_s^2}{2\pi}
\sum_{X_{\bar{c}}}\int \!d\mbox{PS}_{X_{\bar{c}}}\,
\langle 0|\mathcal{Q}_{i'\mu\nu}^{\dagger B}(x,r')
|X_{\bar{c}}\rangle 
\langle X_{\bar{c}}|\mathcal{Q}_{i}^{A\mu\nu}(0,r)|0\rangle 
\nonumber\\
&&\hspace*{0cm}=\,
\delta^{AB}\,(d-2)^2\int\!\frac{d^{d}p}{(2\pi)^d}\,e^{-ipx}
\,\bigg\{
\delta_{i i'} \mathcal{J}_{\bar{c}}^{q\bar{q}(8)}(p^2,r,r')
+(1-\delta_{i i'})\,
\mathcal{\widehat{J}}_{\bar{c}}^{q\bar{q}(8)}(p^2,r,r^\prime)
\bigg\}\,,\quad
\label{eq:NLPqqbarjet}
\end{eqnarray}
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 $\mathcal{O}(\alpha_s)$, 
the $X_{\bar{c}} = q\bar{q}$ final state results in 
\begin{align}
\mathcal{J}_{\bar{c}}^{q\bar{q}(8)}(p^2,r,r')=\frac{\alpha_{s}}{4\pi}
\,\frac{1}{\Gamma\left(1-\epsilon\right)}\,\delta(r-r')\,r\bar{r}\left(\frac{\mu^{2}e^{\gamma_E}}{p^{2}r \bar{r}}\right)^{\!\epsilon},
\end{align}
where $\bar r =1-r$. 

The hatted function, which describes the 
interference of the two B1 operators, is given as this order by 
$\mathcal{\widehat{J}}_{\bar{c}}^{q\bar{q}(8)}(p^2,r,r') = 
(-\epsilon)/(1-\epsilon)\,\mathcal{J}_{\bar{c}}^{q\bar{q}(8)}(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 $\epsilon\to 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 
$\mathcal{\widehat{J}}_{\bar{c}}^{q\bar{q}(8)}(p^2,r,r^\prime)$ 
in the following. We will see later that it is irrelevant for 
the discussion of endpoint divergences.

\paragraph{Soft function} 
Performing the soft decoupling field redefinition on 
\eqref{eq:btypeME} results in the field product
\begin{equation}
\bar{\chi}_{\bar{c}a} [Y_{n_+}^\dagger Y_{n_-}]_{ac} 
\big[t^B\big]_{cd}\, 
[Y_{n_-}^\dagger Y_{n_+}]_{db} \chi_{\bar{c}b}\,
\mathcal{A}^B_{c\perp\nu}\,.
\end{equation}
The Wilson lines can be combined into the adjoint Wilson 
lines $\mathcal{Y}_{n_\mp}$ using 
\begin{equation}
\big[Y_{n_\mp} t^A Y_{n_\mp}^\dagger\big]_{ab}  = 
\big[Y_{n_\mp}^\dagger t^A Y_{n_\mp}\big]_{ab}  = 
\mathcal{Y}_{n_\mp}^{AB} \big[t^B\big]_{ab}\,.
\end{equation}
After squaring the matrix element and summing over the soft 
final state, we obtain the LP two-hemisphere soft 
function in the adjoint representation:
\begin{equation}
S^{(g)}(l_+,l_-)=\frac{1}{N_c^2-1}\langle 0|\overline{T}\left\{
\mathcal{Y}_{n_+}^{BD}(0)\mathcal{Y}_{n_-}^{DA}(0)\right\}
\,\mathcal{P}_s(l_+,l_-)\,
T\left\{\mathcal{Y}_{n_-}^{AC}(0)\mathcal{Y}_{n_+}^{CB}(0)
\right\}|0\rangle \,.
\end{equation}
The LO result is $S^{(g)}(l_+,l_-)=\delta(l_+)\delta(l_-)$.

\paragraph{Factorization formula}

Inserting these matrix element definitions, we 
obtain the factorization formula for the B-type term 
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{d\sigma}{dM_R^2 dM_L^2}|_{\rm B-type} &=& 
\frac{2 C_F}{Q^2} \,f(\epsilon)\int_0^\infty dl_+ dl_-\,
\sum_{i,i'=1,2} \int dr dr'\,
C^{\rm B1*}_{i'}(Q^2,r')C^{\rm B1}_{i}(Q^2,r) 
\nonumber\\
&&\hspace*{-2cm}\times\,
\bigg\{\,
\delta_{i i'}\,\mathcal{J}_{\bar{c}}^{q\bar{q}(8)}(M_R^2-Ql_+,r,r^\prime)
+\,(1-\delta_{i i'})\,\mathcal{\widehat{J}}_{\bar{c}}^{q\bar{q}(8)}(M_R^2-Ql_+,r,r^\prime)
\bigg\}
\nonumber\\[0.2cm]
&&\hspace*{-2cm}\times\,
\mathcal{J}_{c}^{(g)}(M_L^2-Q l_-)\,S^{(g)}(l_+,l_-)\,. 
\label{eq:Btypefactformula}
\end{eqnarray}
In Laplace space:
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{\rm B-type} 
&=& 
\frac{2 C_F}{Q^2}  \,f(\epsilon)\sum_{i,i'=1,2} \int dr dr'\,
C^{\rm B1*}_{i'}(Q^2,r')C^{\rm B1}_{i}(Q^2,r) 
\nonumber\\
&&\hspace*{-2.5cm}\times\,
\bigg\{\,
\delta_{i i'}\,\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)
+\,(1-\delta_{i i'})\;\;\widetilde{\!\!\mathcal{\widehat{J}}}{}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)
\bigg\}\,
\widetilde{\mathcal{J}}_{c}^{(g)}(s_L)\,\widetilde{S}^{(g)}(s_R,s_L)\,. \quad
\label{eq:BtypeLaplace}
\end{eqnarray}
Only the first term in the curly brackets proportional to 
$\delta_{i i'}$ is non-vanishing at leading 
$\mathcal{O}(\alpha_s)$.\footnote{Since the {\em sum} of 
the A-type and B-type term are finite in the limit $\epsilon\to 0$, 
the common prefactor $f(\epsilon)$ could be set to its 
four-dimensional value 1 everywhere at this point.} 

\subsection{Tree-level evaluation}
\label{sec:treelevelevaulation}

Inserting tree-level results for all functions given above, the
two terms in equations \eqref{eq:Atypefactformula} and
\eqref{eq:Btypefactformula} reduce to the following expressions. The
A-type soft-quark term is given by 
\begin{eqnarray} 
\label{eq:AtypesqTree}
\frac{1}{\sigma_0}\frac{d^2\sigma}{dM^2_R\,dM^2_L}|_{
\scriptsize \begin{array}{l}
$\rm A--type$,\\[-0.1cm] \rm tree \end{array}} 
&=&\frac{2C_F}{Q}\, f(\epsilon) \int dl_+\,dl_-
\, \int d\omega d\omega' \, \delta^+(M^2_R-Ql_+ )\,
\nonumber\\[0ex] 
&&\hspace*{-3cm} 
\times \,\frac{1}{\omega\omega'}\,\delta^+(M^2_L-Ql_-) 
\nonumber \,\frac{\alpha_s}{4\pi} \,\frac{\theta(\omega)
\delta(\omega-\omega') }{e^{-\epsilon\gamma_E}\Gamma(1-\epsilon)}
\,\bigg[ \delta(l_- ) \,\theta(\omega-l_+ ) \theta(l_+ )\,
\omega \left(\frac{l_+ \omega}{\mu^2}\right)^{\!-\epsilon}
\\[0ex] 
&& \hspace{-2cm}- \frac{ \omega^{2}}{1-\epsilon} 
\, \delta(l_+ )\,\delta(l_- -\omega)
\left(\frac{ \omega^2}{\mu^2}\right)^{\!-\epsilon}\bigg]
\nonumber\\
&&\hspace{-4cm} =
\frac{\alpha_s}{4\pi}\,
\frac{2C_F}{Q^2} \,\frac{f(\epsilon)}{e^{-\epsilon\gamma_E}\Gamma(1-\epsilon)}
\bigg[ \frac{1}{\epsilon}\delta^+(M^2_L)
\left(\frac{({M^2_R})^2}{{Q}^2\mu^2}\right)^{\!-\epsilon}
- \,\delta^+(M^2_R )\,\frac{ 1}{1-\epsilon}
\left(\frac{ (M^2_L)^2}{Q^2\mu^2}\right)^{\!-\epsilon}\bigg],\quad\;
\end{eqnarray}
where the pole in $\epsilon$ arises from the logarithmic divergence
of the integral $\int_{l_+}^\infty d\omega/\omega$ 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
\begin{eqnarray}
\label{eq:BtypeTree}
\frac{1}{\sigma_0}\frac{d^2\sigma }{dM^2_L\,dM^2_R}
|_{\scriptsize \begin{array}{l}
$\rm B--type$,\\[-0.1cm] \rm tree \end{array}}
&=& \frac{2C_F}{Q^2 } \,f(\epsilon) \int dl_+\,dl_- \int_0^1\, 
dr dr'\,\delta^+(M^2_L-Ql_-) \delta(l_+)\delta(l_-)
\nonumber\\
&&\hspace{-3.5cm} \times\,\bigg\{
\left(\frac{1}{r}\frac{1}{r'}+\frac{1}{\bar{r}}\frac{1}{\bar{r}'} 
\right) 
+
\left(\frac{1}{r}\frac{1}{\bar{r}'}+\frac{1}{\bar{r}}\frac{1}{r'} 
\right)\frac{\epsilon}{1-\epsilon}
\bigg\}
\,\delta(r -r' ) 
\,\frac{\alpha_s}{4\pi}
\frac{ r \bar{r} }{\Gamma(1-\epsilon)}
\bigg(\frac{(M^2_R-Ql_+)\, {r} \bar{r} }{\mu^2e^{\gamma_E}} 
\bigg)^{\!-\epsilon}
\nonumber\\[0.2cm]
&&\hspace{-3.5cm} = \frac{\alpha_s}{4\pi}\,
\frac{4C_F}{Q^2 }\,f(\epsilon)\,\bigg\{
-\frac{1}{\epsilon} + 
\frac{\epsilon}{(1-\epsilon)^2 }\,\bigg\} \, \delta^+(M^2_L)
 \bigg(\frac{M^2_R}{\mu^2e^{ \gamma_E}} \bigg)^{\!-\epsilon}
\frac{\Gamma (2-\epsilon)}{\Gamma (2-2 \epsilon)}\,.
\end{eqnarray}
In this case, the singularity arises from a logarithmic divergence
in the integral over the momentum fraction $r$ as $r\to 0$ and 
$r\to 1$. Adding both terms and taking the limit $\epsilon \to 0$, 
we find
\begin{eqnarray}
\label{eq:treeresultA+B}
\frac{1}{\sigma_0}\frac{d^2\sigma }{dM^2_L\,dM^2_R} &=&
\frac{\alpha_s C_F}{\pi}\,\frac{1}{Q^2 }
\,\bigg\{\delta^+(M^2_L) \left[\ln\frac{Q^2}{M_R^2} -1\right] 
- \delta^+(M^2_R)\bigg\}\,.
\end{eqnarray}
The poles cancel in the sum of both the A and B-type contributions
in equations \eqref{eq:AtypesqTree} and \eqref{eq:BtypeTree}
(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 two-hemisphere invariant mass distribution to the
thrust distribution according to \eqref{eq:taudist}, we reproduce
the coefficient $-\alpha_s C_F/\pi$ of the 
$\ln \tau$ term in \cite{Moult:2016fqy}.

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
$\epsilon \to 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 $\omega, \omega'$
in \eqref{eq:AtypesqTree} in the limit of 
large $\omega, \omega'$ (where the divergence arises) is given by
\begin{eqnarray} 
\label{eq:AtypesqTreeAsy}
&&\frac{\alpha_s}{4\pi}\,\frac{2 C_F}{Q}\, f(\epsilon)\int dl_+dl_-\,
\delta^+(M^2_R-Ql_+ )\delta^+(M^2_L-Ql_-)
\delta(l_- ) \,\theta(l_+ )\, \,
\nonumber\\ 
&& \hspace*{-0cm} \times\frac{\delta(\omega-\omega')}{\omega\omega'}\,
\nonumber  \,\frac{1}{e^{-\epsilon\gamma_E}\Gamma(1-\epsilon)}
\,
\omega \left(\frac{l_+ \omega}{\mu^2}\right)^{\!-\epsilon}
\nonumber\\ 
&&= \,\frac{\alpha_s}{4\pi}\,\frac{2 C_F}{Q^2}\,\delta^+(M^2_L)\, 
\frac{\delta(\omega-\omega')}{\omega\omega'}\,
\frac{f(\epsilon)}{e^{-\epsilon\gamma_E}\Gamma(1-\epsilon)}
\,
\omega \left(\frac{M_R^2\omega}{Q\mu^2}\right)^{\!-\epsilon}.
\end{eqnarray}
In the B-type term \eqref{eq:BtypeTree}, the divergence of the 
convolution arises from small momentum fraction. The 
integrand of the $r, r'$ integrals in the limit of small $r$ 
reads 
\begin{eqnarray}\label{eq:BtypeTreeAsy}
&&\frac{\alpha_s}{4\pi}
\,\frac{2 C_F}{Q^2 } f(\epsilon)\int dl_+dl_-\,
\,\delta^+(M^2_L-Ql_-) \delta(l_+)\delta(l_-)
\nonumber\\
&&\times\,\frac{\delta(r -r' )}{r r'} 
\frac{1}{\Gamma(1-\epsilon)}\,r\,
\bigg(\frac{(M^2_R-Ql_+)\, {r} }{\mu^2e^{\gamma_E}} 
\bigg)^{\!-\epsilon}
\nonumber\\ 
&&= \frac{\alpha_s}{4\pi}
\,\frac{2 C_F}{Q^2 }\,\delta^+(M^2_L)\,
\frac{\delta(r -r' )}{r r'} 
\frac{ f(\epsilon)}{\Gamma(1-\epsilon)}\,r\,
\bigg(\frac{M^2_Rr }{\mu^2e^{\gamma_E}} 
\bigg)^{\!-\epsilon}\,,
\end{eqnarray}
which is identical to the previous expression provided we 
identify $r=\omega/Q$, $r'=\omega'/Q$. A similar agreement 
holds for the identical soft anti-quark contribution to the 
A-type term and the $r\to 1$ limit of the B-type 
term \eqref{eq:BtypeTree}. From the heuristic discussion of 
Sec.~\ref{sec:heuristic} it is evident that this agreement is 
not a coincidence, but holds to all orders in the expansion 
in $\alpha_s$.


\section{Endpoint factorization}
\label{sec:refactorization}

Eqs.~\eqref{eq:Atypefactformula}, \eqref{eq:Btypefactformula} 
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 $\alpha_s$ 
by evaluating the convolutions of the $d$-dimensional 
hard, (anti-)collinear and soft function. The limit 
$\epsilon\to 0$ must be taken at the end, because the 
convolutions in the A-type term are divergent when 
$\omega$ and $\omega^\prime$ both tend to infinity. In the B-type 
term this happens when $r$ and $r^\prime$ both approach 0 or 1.

This form is not quite sufficient to sum the logarithms 
to all orders in $\alpha_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 $\epsilon\to 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 \cite{Beneke:2008pi}, which involves SCET and 
non-relativistic QCD, and Higgs decay to two photons 
\cite{Liu:2020wbn}, which is a SCET$_{\rm II}$ process. 
These are amplitude-factorization problems. The present 
case of gluon thrust is somewhat different as it is a 
SCET$_{\rm I}$ problem, and involves inclusive collinear and 
soft functions defined at cross-section level, similar 
to \cite{Beneke:2020ibj}. 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\to\gamma\gamma$. 

\subsection{Soft-collinear limit of the B1 matching 
coefficients}

An important ingredient in this discussion 
is the factorization property of the coefficient function of 
$q\bar q g$ 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_i^{\rm B1}$ matching 
coefficients in this limit, which in turn arises because the 
intermediate quark or anti-quark goes on-shell, see 
Fig.~\ref{fig:QCD}. 

By its original definition, the matching coefficient is a 
single-scale, hard function. For $r\to 0$ (soft quark) 
and $r\to 1$ (soft anti-quark), it becomes a two-scale object, 
which can itself be factorized according to 
\cite{Beneke:2020ibj}
\begin{align}
& \hskip0.5cm\displaystyle
C_1^{\rm B1}(Q^2,r) = C^{\rm A0}(Q^2)\times 
\frac{D^{\rm B1}(r Q^2)}{r} 
+\mathcal{O}(r^0)\,,
\label{eq:B1fact1}
\end{align}
where the two scales $Q^2$ and $r Q^2$ are now separated 
into the hard matching coefficient of the leading-power 
A0 SCET current, and a new coefficient $D^{\rm B1}(r Q^2)$, 
which depends only on the endpoint-scale $\sqrt{r} Q$. 
We recall that $C_1^{\rm B1}(Q^2,r)$ is regular in the 
soft anti-quark limit  $r\to 1$. An analogous factorization 
holds for the second B1 operator in the limit $r\to 1$:
\begin{align}
& \hskip0.5cm\displaystyle
C_2^{\rm B1}(Q^2,r) = -C^{\rm A0}(Q^2)\times 
\frac{D^{\rm B1}(\bar{r} Q^2)}{\bar r} +\mathcal{O}(\bar{r}^0)\,.
\label{eq:B1fact2}
\end{align}
Since the two cases are related by $r\leftrightarrow \bar{r}$ 
(plus a global minus sign), the following is phrased for the 
soft-quark limit only. A graphical representation of 
\eqref{eq:B1fact1} is shown in Fig.~\ref{fig:Cb1refact}. 

\begin{figure}
\begin{center}
\includegraphics[width=0.7\textwidth]{figs/CB1toDB1fact}
\caption{
Graphical representation of the factorization \eqref{eq:B1fact1} 
of the B1 matching coefficient in the soft-quark limit.
\label{fig:Cb1refact}}
\end{center}
\end{figure}

For $r\to 0$ the outgoing quark represented by 
$\bar\chi_{\bar c}$ field in the second line of 
\eqref{eq:hardmatching} becomes soft-anticollinear. The 
soft-collinear coefficient $D^{\rm B1}(p^2)$, which appears in 
\eqref{eq:B1fact1}, \eqref{eq:B1fact2}, can be defined by matching 
the time-ordered product \cite{Beneke:2020ibj} 
\begin{equation}
\langle g^a_{c}(p_{c}) 
q_{\overline{sc}}(p_{\overline{sc}})| 
\int d^4x \,T\big\{\bar{\chi}_{c}(0),\mathcal{L}_{\xi q}(x) 
\big\} \,|0\rangle=g_s
 \bar{u}(p_{\overline{sc}}) t^a\slashed{\varepsilon}_{c\perp}(p_c) 
\frac{i n_+ p_c}{p^2} \frac{\slashed n_-}{2} \,D^{\rm B1}(p^2)\,,
\label{eq:DB1def}
\end{equation}
with $p^2=(p_{c}+p_{\overline{sc}})^2=r Q^2$ and
\begin{align}
\mathcal{L}_{\xi q}(x)=\bar{q}_{\overline{sc}}(x_{-})
\slashed{\mathcal{A}}_{c\perp}(x)\chi_{c}(x)+\mbox{h.c.}
\end{align}
the leading soft-quark Lagrangian \cite{Beneke:2002ph,Beneke:2002ni} 
with $\bar{q}_s$ replaced by $\bar{q}_{\overline{sc}}$. 
This relation becomes evident by comparing the second diagram in 
Fig.~\ref{fig:SCET} to the right-hand side in 
Fig.~\ref{fig:Cb1refact}. From its definition in terms of a 
time-ordered product it is clear that the matching coefficient 
$D^{\rm B1}(p^2)$ is a universal function that will appear 
in different processes involving soft quark emission.

The leading double-logarithmic resummation of 
$D^{\rm B1}(p^2)$ to all orders in $\alpha_s$ has been derived  
in \cite{Beneke:2020ibj}. It exhibits the all-order 
$(C_A-C_F)^n$ colour coefficient, which seems to be characteristic 
for soft quark 
emission \cite{Vogt:2010cv,Moult:2019uhz,Beneke:2020ibj}. 
Interestingly, the same coefficient 
enters the $ggH$ amplitude with a bottom-quark loop and was 
recently computed at the two-loop level \cite{Liu:2021mac}, 
together with its one-loop evolution kernel, which was obtained 
from renormalization-group invariance of the full $ggH$ amplitude. 

The $D^{\rm 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\to 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 \cite{Beneke:2021pkl}, see 
also \cite{Liu:2020wbn}. We refer to App.~\ref{sec:DB1} for 
this derivation, which results in 
\begin{equation}
D^{\rm B1}(p^2) = 1 + \frac{\alpha_s}{4\pi}\,
\big(C_F- C_A \big)
\bigg(\frac{2}{\epsilon^2} -1 - \frac{\pi^2}{6} \bigg)
\left(\frac{\mu^2}{-p^2-i\varepsilon}\right)^{\!\epsilon}+\mathcal{O}(\alpha_s^2).
\label{eq:DB11loop}
\end{equation}
\begin{equation}
\frac{d}{d\ln\mu}D^{\rm B1}(p^2) = \int_0^\infty d\hat p^2
\,\gamma_D(\hat p^2,p^2)D^{\rm B1}(\hat p^2)\,,
\label{eq:DBQrge}
\end{equation}
with 
\begin{eqnarray}
\gamma_D(\hat p^2,p^2) &=& 
\frac{\alpha_s(C_F-C_A)}{\pi}\,\delta(\hat p^2-p^2)
\ln\left(\frac{\mu^2}{-p^2-i\varepsilon}\right)\nonumber\\
  && + \,\frac{\alpha_s}{\pi} \left(\frac{C_A}{2}-C_F\right)
p^2 \left[\frac{\theta(\hat p^2 -p^2)}
{\hat p^2 (\hat p^2 -p^2)}
+ \frac{\theta(p^2-\hat p^2 )}{p^2(p^2-\hat p^2 )}\right]_+  
\,,\quad
\label{eq:DB1andim}
\end{eqnarray}
in agreement with \cite{Liu:2021mac}.

When the gluon field is replaced by a photon, 
the corresponding abelian version of $D^{\rm B1}(p^2)$ 
is identical to the jet function, which appears at 
LP in the factorization of the 
radiative semi-leptonic $B$ decay $B\to\gamma\ell\nu$ 
\cite{Bosch:2003fc} and in the NLP endpoint factorization of the 
photonic Higgs decay $H\to\gamma\gamma$ amplitude through 
a bottom-quark loop \cite{Liu:2020wbn}. In the abelian 
case, the two-loop evolution has been inferred from 
renormalization-group consistency of the $B\to\gamma\ell\nu$ 
observable \cite{Liu:2020ydl}. The direct computation of 
the one-loop evolution kernel of this jet function can be 
found in \cite{Bodwin:2021epw}.

\subsection{Endpoint factorization consistency conditions}

As discussed in Sec.~\ref{sec:heuristic} and checked explicitly above 
at $\mathcal{O}(\alpha_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_-\cdot p_{\bar c} = r Q, r'Q \to 0$ in the B-type term, 
matches the limit $n_-\cdot k = \omega,\omega' 
\to \infty$ of the corresponding soft momentum component in 
the A-type term.\footnote{There is a similar matching 
of the limit $(1-r)Q,(1-r')Q\to 0$ with the 
soft {\em anti}-quark contribution to the A-type term. 
This contribution can be treated as below, after 
change variables from $r^{(\prime)}\to 
1-r^{(\prime)}$.}

We start from the expressions \eqref{eq:AtypeLaplace}, 
\eqref{eq:BtypeLaplace} in Laplace space, in which the 
integrand 
of the $\omega,\omega'$ 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:
\begin{eqnarray}
&&\frac{2C_F}{Q^2} \,f(\epsilon) \,
C^{\rm B1*}_{1}(Q^2,r')C^{\rm B1}_{1}(Q^2,r) \,
\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)\,
\widetilde{\mathcal{J}}_{c}^{(g)}(s_L)\,\widetilde{S}^{(g)}(s_R,s_L)\Big|_{
r^{(\prime)}=\frac{\omega^{(\prime)}}{Q}\to \,0}
\nonumber\\
&& \hspace*{0.5cm} 
= \frac{2C_F}{Q} \,f(\epsilon)\,|C^{\rm A0}(Q^2)|^2 \,
\widetilde{\mathcal{J}}_{\bar c}^{(\bar q)}(s_R)
\,\times\,
\bigg\{\,\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
\,\widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')
\nonumber\\[0.1cm]
&& \hspace*{1cm} +\;\,
\widetilde{\!\!\widehat{\mathcal{J}}}_{\!c}(s_{L},
\omega,\omega')
\,\,\widetilde{\!\widehat{S}}_{\rm NLP}(s_R,s_L,\omega,\omega')\,
\bigg\}\Big|_{\omega,\omega^{\prime} \to\infty}\,.
\label{eq:reff1}
\end{eqnarray} 
The left-hand side of this equation 
represents the leading asymptotics of the 
B-type term \eqref{eq:BtypeLaplace} 
for $r, r' \to 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\to 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' \to \bar r,\bar{r}^\prime$ 
and $C^{\rm B1}_{1} \to C^{\rm B1}_{2}$ applies to the 
limit $r,r'\to 1$. The right-hand side of \eqref{eq:reff1} 
contains two terms in the curly brackets, which 
correspond to different colour structure in the jet 
and soft function, \eqref{eq:defjetNLP} and 
\eqref{eq:NLPsoftfndef}, respectively. We shall now show that 
endpoint factorization consistency implies that the 
hatted term  in \eqref{eq:reff1} is subleading as 
$\omega,\omega^\prime\to \infty$, and can be dropped.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.8\textwidth]{figs/jetfn_1PR}
\caption{Factorization of the radiative collinear function 
$\mathcal{J}_{c}(p^2,\omega,\omega')$, which appears in the 
A-type term, in the limit $\omega, \omega^\prime\to\infty$.
\label{fig:jet1PR}}
\end{center}
\end{figure}


We recall from \eqref{eq:defjetNLP} that 
$\mathcal{J}_{c}(p^2,\omega,\omega')$ is an inclusive 
radiative jet function of the non-local, time-ordered product 
operator defined in \eqref{eq:nonlocaljetfield}. 
Graphically, we may represent it by the diagram to the 
left of the arrow in Fig.~\ref{fig:jet1PR}, 
where the dashed external lines to the left and right 
of the cut denotes the {\em soft} momentum of the 
composite field $\not\hspace*{-1.5mm}\mathcal{A}_{\perp c}\chi_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/(\omega\omega^\prime)$ for 
large $\omega,\omega^\prime$ requires that the diagram 
is ``one-collinear-particle reducible'', as displayed to 
the right of the arrow in Fig.~\ref{fig:jet1PR}. Then 
colour conservation implies that the collinear line must 
be a colour-octet, that is a gluon, and hence only the 
unhatted colour term \eqref{eq:defjetNLP}, which defines 
$\mathcal{J}_{c}(p^2,\omega,\omega')$ 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.~\ref{fig:jet1PI}.
Endpoint factorization consistency requires that the 
large $\omega,\omega^\prime$ 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 \eqref{eq:hardmatching}. 
Consistency therefore requires that the asymptotic behaviour 
takes the diagrammatic form shown in Fig.~\ref{fig:jet1PR}. 
This proves not only that the hatted jet function is 
irrelevant for endpoint factorization, but further 
that $\mathcal{J}_{c}(p^2,\omega,\omega')$ 
factorizes into the product\footnote{The collinear 
scale on the left-hand side is $n_+ p \,\omega^{(\prime)} 
=Q  \omega^{(\prime)}$, but we suppressed the factor of 
$Q$ in the respective arguments of 
$\mathcal{J}_{c}\left(p^2,\omega,\omega'\right)$.}
\begin{eqnarray}
\mbox{(I)} &\hskip0.5cm 
\mathcal{J}_{c}\left(p^2,\omega,\omega'\right) = 
\mathcal{J}_{c}^{(g)}(p^2) \,
\displaystyle \frac{D^{\rm B1}(\omega Q)}{\omega}
  \frac{{D^{\rm B1}}^*(\omega^\prime Q)}{\omega^\prime} 
+\mathcal{O}\!\left(\frac{1}{\omega^{(\prime)}}\right)\,,
\label{eq:EPfactI}
\end{eqnarray}
where the function $D^{\rm B1}(p^2)$ is the same as the one  
that appears in the factorization of the hard B1 operator 
coefficient \eqref{eq:B1fact1}. 

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.3\textwidth]{figs/jetfn_1PI}
\caption{Example of a one-collinear-particle reducible diagram, 
which does not contribute to the large-$\omega^{(\prime)}$ 
behaviour of $\mathcal{J}_{c}(p^2,\omega,\omega')$.
\label{fig:jet1PI}}
\end{center}
\end{figure}

Eq.~\eqref{eq:EPfactI} is the first endpoint factorization 
consistency condition. Inserting the Laplace-transformed version 
of this relation together with 
\eqref{eq:B1fact1} for the B1 hard function 
into \eqref{eq:reff1} immediately implies 
a similar consistency relation for the NLP soft quark 
function, 
\begin{eqnarray}
\mbox{(II)} & \hskip-1.5cm 
Q\,\widetilde{\mathcal{J}}^{(\bar{q})}_{\bar{c}}(s_R)\,
\widetilde{S}_{\rm NLP} \left(s_R,s_L,\omega,\omega'\right)\Big|_{\omega^{(\prime)} \to\infty}
\nonumber\\
& \hskip2.5cm \displaystyle
= 
\widetilde{\mathcal{J}}^{q\bar q(8)}_{\bar{c}}\left(s_R,r,r'\right)
\widetilde{S}^{(g)}(s_R,s_L)\Big|_{r^{(\prime)}=\omega^{(\prime)}/Q\to 0}\,.
\end{eqnarray}
The same identity holds 
with $r,r^\prime\to \bar{r},\bar{r}^\prime$. Relations (I) and (II) formalize the heuristic picture developed 
in Sec.~\ref{sec:heuristic}. For large $\omega$, the soft 
quark field in $S_{\rm NLP}$ becomes anti-collinear, 
$\bar{q}_s Y_{n_-} \to \bar{\chi}_{\bar c}$, hence it moves 
from $S_{\rm NLP}$ to the anti-collinear function 
$\mathcal{J}^{q\bar q(8)}_{\bar{c}}$. Removing  $\bar{q}_s Y_{n_-}$ 
from  $S_{\rm NLP}$ leaves the leading-power soft function 
$S^{(g)}$, while adding it as $\bar{\chi}_{\bar c}$ to 
$\mathcal{J}^{(\bar{q})}_{\bar{c}}$ turns the anti-quark 
jet function into $\mathcal{J}^{q\bar q(8)}_{\bar{c}}$. Overall, 
this results in $S_{\rm NLP} \mathcal{J}^{(\bar{q})}_{\bar{c}} \to 
S^{(g)}\mathcal{J}^{q\bar q(8)}_{\bar{c}}$, which is relation 
(II). At the same time, the quark fields in the A-type 
collinear function \eqref{eq:nonlocaljetfield} become highly 
off-shell, which removes them from $\mathcal{J}_{c}(p^2,\omega,\omega^\prime)$, leaving 
only the collinear gluon, and consequently $C^{\rm A0}$ turns 
into $C^{\rm B1}$. Thus, $|C^{\rm A0}|^2 \mathcal{J}_{c} \to 
|C^{\rm B1}|^2 \mathcal{J}_{c}^{(g)}$, which is relation (I).

\subsection{Endpoint factorization formula}
\label{sec:endpointfact}

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 \cite{Liu:2019oav} 
to denote the asymptotic behaviours of the various functions. 
The precise definitions are as follows: In functions of 
$\omega, \omega^\prime$, rescale $\omega \to \kappa 
\omega$, $\omega^\prime \to \kappa 
\omega^\prime$ and take $\kappa\to \infty$. Then 
\begin{eqnarray}
&& \llbracket 
S_{\rm NLP} \left(l_+,l_-,\omega,\omega'\right)
\rrbracket \equiv S_{\rm NLP} \left(l_+,l_-,\omega,\omega'\right)|_{
\mathcal{O}(\kappa^0)}, \label{eq:Sasym}\\
&&\llbracket \mathcal{J}_{c}(p^2,\omega,\omega^\prime)
\rrbracket \equiv \mathcal{J}_{c}(p^2,\omega,\omega^\prime)|_{
\mathcal{O}(\kappa^{-2})}\,,
\end{eqnarray}
where $\delta(\omega-\omega')$ counts as $\kappa^{-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^\prime$, rescale $r \to r \kappa$, $r^\prime \to 
r^\prime \kappa$ and take $\kappa\to 0$, {\em or} the 
corresponding rescaling is applied to $\bar r$, $\bar{r}^\prime$. 
Which of the two is meant, will be indicated by the subscript 
0 or 1 on the double bracket. Then 
\begin{eqnarray}
&& \llbracket C^{\rm B1}_1(Q^2,r) \rrbracket_0 \equiv 
C^{\rm B1}_1(Q^2,r)|_{\mathcal{O}(\kappa^{-1})}\,,
\\[0.1cm]
&& \llbracket C^{\rm B1}_2(Q^2,r) \rrbracket_1 \equiv 
C^{\rm B1}_2(Q^2,r)|_{\mathcal{O}(\kappa^{-1})}\,,
\end{eqnarray}
while $\llbracket C^{\rm B1}_1(Q^2,r) \rrbracket_1 = 
\llbracket C^{\rm B1}_2(Q^2,r) \rrbracket_0 = 0$. With this 
definition the right-hand sides of these equations are 
given by \eqref{eq:B1fact1} and \eqref{eq:B1fact2}, respectively. 
Finally, we have 
\begin{equation}
\llbracket \mathcal{J}^{q\bar q(8)}_{\bar{c}}(p^2,r,r^\prime)
\rrbracket_0 = 
\llbracket \mathcal{J}^{q\bar q(8)}_{\bar{c}}(p^2,r,r^\prime)\rrbracket_1 
\equiv \mathcal{J}^{q\bar q(8)}_{\bar{c}}(p^2,r,r^\prime)|_{
\mathcal{O}(\kappa^0)}
\end{equation}
due to the symmetry $r^{(\prime)}\leftrightarrow \bar{r}^{(\prime)}$.
In the $d$-dimensional expressions there are 
$\epsilon$-dependent powers of $\kappa$, which turn into 
logarithms in the renormalized functions. Thus 
$\epsilon$ counts as an 
infinitesimal variable for the purpose of endpoint 
$\kappa$ power-counting.

To implement the rearrangement of endpoint-singular terms 
we start with the scaleless integral
\begin{eqnarray}
&&\frac{2 C_F}{Q} \,f(\epsilon)\,|C^{\rm A0}(Q^2)|^2 \widetilde{\mathcal{J}}^{(\bar q)}_{\bar{c}}(s_R) 
\widetilde{\mathcal{J}}_{c}^{(g)}(s_L) \nonumber\\
&&\times \int_0^\infty d\omega d\omega'\,
\frac{D^{\rm B1}(\omega Q)}{\omega}
\frac{D^{\rm B1*}(\omega'Q)}{\omega'} 
\left \llbracket \widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')\right 
\rrbracket,\qquad
\label{eq:scaleless_integral}
\end{eqnarray}
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.~\eqref{fig:overlap}: 
(1) $\omega$ {\em and} $\omega'$ smaller than 
an endpoint factorization parameter $\Lambda$ 
(integral $I_1$), and the complement region $I_2$. 
(2)  $\omega$ {\em or} $\omega'$ smaller $\Lambda$ 
and the complement region. Since the mixed regions 
where one variable is larger and the other smaller than 
$\Lambda$ 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 $\omega,\omega^\prime > \Lambda$ and the 
double-bracket asymptotic behaviour can be used for  
functions of $\omega,\omega'$ in the A-type term.\footnote{
While in version (1), it could be used for functions 
of $r, r'$ in the B-type term with identification 
$r^{(\prime)}=\omega^{(\prime)}/Q$. The resulting 
expressions for this version 
are provided in App.~\ref{version1endpintfact}.} 
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 $\Lambda$. However, as long as no approximations 
are made, the $\Lambda$ dependence cancels exactly between 
the two terms.

\begin{figure}
\begin{center}
\includegraphics[width=0.85\textwidth]{figs/subtractionversions}
\caption{\label{fig:overlap}
Version (1) [left] and (2) [right] of the split of 
\eqref{eq:scaleless_integral} into $I_1+I_2$ according to 
the correspondingly indicated regions in the 
$\omega-\omega'$ plane as described in the text 
below \eqref{eq:scaleless_integral}. In the overlap 
region in green the asymptotic behaviour of the A- 
and B-type term must agree.}
\end{center}
\end{figure}

We begin with the A-type contribution 
\eqref{eq:AtypeLaplace}, which accounts for soft-quark 
emission, to the factorization theorem and subtract from it the complement region $I_2$, $\omega,\omega^\prime > \Lambda$,  of the integral~(\ref{eq:scaleless_integral}), resulting in
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{\rm A-type} &=& 
\frac{2C_F}{Q} \,f(\epsilon)\,|C^{\rm A0}(Q^2)|^2 \,
\widetilde{\mathcal{J}}_{\bar c}^{(\bar q)}(s_R)
\,\int_0^\infty d\omega d\omega'\,
\nonumber\\
&&\hspace*{-3cm}\times\,
\bigg\{\,\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
\,\widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')
\nonumber\\[-0.15cm]
&&\hspace*{-2.2cm}
-\,\theta(\omega-\Lambda)\theta(\omega'-\Lambda)
\widetilde{\mathcal{J}}_{c}^{(g)}(s_L) 
\frac{D^{\rm B1}(\omega Q)}{\omega}
\frac{D^{\rm B1*}(\omega'Q)}{\omega'} 
\left \llbracket \widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')\right 
\rrbracket
\nonumber\\
&&\hspace*{-2.2cm}
+\;\,\widetilde{\!\!\widehat{\mathcal{J}}}_{\!c}(s_{L},
\omega,\omega') 
\,\,\widetilde{\!\widehat{S}}_{\rm NLP}(s_R,s_L,\omega,\omega')\,\bigg\}\,. \quad
\label{eq:Atype_subtracted_general}
\end{eqnarray}
For $\omega,\omega'\to\infty$, we use the endpoint 
factorization condition (I) from \eqref{eq:EPfactI} 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 (\ref{eq:Sasym}). Thus the 
first and second term in the integrand cancel. The third 
term is subleading in this limit, and  the convolution 
integrals over $\omega$ and $\omega'$ are now endpoint-finite, as the logarithmically divergent part of the integrand in the large $\omega$ and $\omega'$ region has been removed. 
To make this explicit, we use  \eqref{eq:EPfactI} in 
\eqref{eq:Atype_subtracted_general}, and obtain 
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{\rm A-type} &=& 
\frac{2C_F}{Q} \,f(\epsilon)\,|C^{\rm A0}(Q^2)|^2 \,
\widetilde{\mathcal{J}}_{\bar c}^{(\bar q)}(s_R)
\,\int_0^\infty d\omega d\omega'\,
\nonumber\\
&&\hspace*{-3cm}\times\,
\bigg\{\,\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
\,\widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')
\nonumber\\[0cm]
&&\hspace*{-2.2cm}
-\,\theta(\omega-\Lambda)\theta(\omega'-\Lambda)\,
\llbracket\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
 \rrbracket
\llbracket \widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')
\rrbracket
\nonumber\\
&&\hspace*{-2.2cm}
+\;\,\widetilde{\!\!\widehat{\mathcal{J}}}_{\!c}(s_{L},
\omega,\omega') 
\,\,\widetilde{\!\widehat{S}}_{\rm NLP}(s_R,s_L,\omega,\omega')\,\bigg\}\,. \quad
\label{eq:Atype_subtracted_general2}
\end{eqnarray}

If we further choose $\Lambda\gg 1/s_{R},1/s_L$, we may simplify 
the previous equation to 
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{\rm A-type} &=& 
\frac{2C_F}{Q} \,f(\epsilon)\,|C^{\rm A0}(Q^2)|^2 \,
\widetilde{\mathcal{J}}_{\bar c}^{(\bar q)}(s_R)
\,\int_0^\infty d\omega d\omega'\,
\nonumber\\
&&\hspace*{-2.5cm}\times\,
\bigg\{\,\Big[1-\theta(\omega-\Lambda)\theta(\omega'-\Lambda)\Big]\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
\,\widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')
\nonumber\\
&&\hspace*{-1.7cm}
+\;\,\widetilde{\!\!\widehat{\mathcal{J}}}_{\!c}(s_{L},
\omega,\omega') 
\,\,\widetilde{\!\widehat{S}}_{\rm NLP}(s_R,s_L,\omega,\omega')\bigg\}\nonumber\\
&&\hspace*{-3cm}=
\,\frac{2C_F}{Q} \,f(\epsilon)\,|C^{\rm A0}(Q^2)|^2 \,
\widetilde{\mathcal{J}}_{\bar c}^{(\bar q)}(s_R)
\,\int d\omega d\omega'\,\Big[1-\theta(\omega-\Lambda)\theta(\omega'-\Lambda)\Big]
\nonumber\\[0.0cm]
&&\hspace*{-2.5cm}\times\,
\bigg\{\,\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
\,\widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')+\,\widetilde{\!\!\widehat{\mathcal{J}}}_{\!c}(s_{L},
\omega,\omega') 
\,\,\widetilde{\!\widehat{S}}_{\rm NLP}(s_R,s_L,\omega,\omega')\bigg\}\,,\qquad
\label{eq:Atype_subtracted}
\end{eqnarray}
where the equality signs hold up to corrections 
of $\mathcal{O}(1/(s_L\Lambda),1/(s_R\Lambda))$. It is then 
understood that in evaluating \eqref{eq:Atype_subtracted}, 
terms suppressed by powers of $\Lambda$ are dropped.
This is simply \eqref{eq:AtypeLaplace} with integration 
region $\omega,\omega^\prime > \Lambda$ removed, see
Fig.~\ref{fig:overlap}. 

The remaining part $I_1$ of the integral (\ref{eq:scaleless_integral}) can now be combined with  $i=i'=1$ part of the 
B-type term (\ref{eq:BtypeLaplace}). 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 \leftrightarrow \bar{r}$. The term $i\neq i'$ is endpoint-finite so it remains unaffected. 

Beginning with the $i=i'=1$ term, we have
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{
\scriptsize \begin{array}{l}
$\rm B--type$\\[-0.1cm] $i=i'=1$\end{array}} 
&=& 
\frac{2 C_F}{Q^2} \, f(\epsilon)\,\bigg[\, 
\widetilde{\mathcal{J}}_{c}^{(g)}(s_L)\,\widetilde{S}^{(g)}(s_R,s_L)
\nonumber\\
&&\hspace*{-1.5cm}\times\,
\int_0^1 dr dr'\,
C^{\rm B1*}_{1}(Q^2,r')C^{\rm B1}_{1}(Q^2,r) 
\,\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)
\,\quad 
\nonumber \\ &&\hspace*{-2.5cm}-\,Q\,|C^{A0}(Q^2)|^2 
\widetilde{\mathcal{J}}^{(\bar q)}_{\bar{c}}(s_R) 
\widetilde{\mathcal{J}}_{c}^{(g)}(s_L) 
\int_0^\infty d\omega d\omega'  \,\big[1-\theta(\omega-\Lambda)\theta(\omega'-\Lambda)\big]
\nonumber\\ &&
\hspace*{-1.5cm}\times\,\frac{D^{\rm B1}(\omega Q)}{\omega}
\frac{D^{\rm B1*}(\omega'Q)}{\omega'} 
\left \llbracket \widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')\right 
\rrbracket\bigg]\,.\qquad 
\label{eq:Btype_subtracted}
\end{eqnarray}
We substitute $\omega = r Q$ and $\omega'=r'Q$ in the subtraction term, and employ (\ref{eq:B1fact1}) to express the $D^{\rm 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 $\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)$ in the asymptotic limit to obtain 
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{
\scriptsize \begin{array}{l}
$\rm B--type$\\[-0.1cm] $i=i'=1$\end{array}} 
&=& 
\frac{2C_F}{Q^2}\,f(\epsilon) \,\widetilde{\mathcal{J}}_{c}^{(g)}(s_L)\,
\widetilde{S}^{(g)}(s_R,s_L)\, \int_0^\infty dr dr' \nonumber\\ 
&&\hspace*{-2.5cm}\times\, \bigg[ \,
\theta(1-r)\theta(1-r')\,
C^{\rm B1*}_{1}(Q^2,r')C^{\rm B1}_{1}(Q^2,r) 
\,\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)
\quad 
\nonumber \\ &&\hspace*{-2cm}-
\big[1-\theta(r-\Lambda/Q)\theta(r'-\Lambda/Q)\big]\,
\nonumber \\ &&\hspace*{-2cm}\times\,
\llbracket C^{\rm B1*}_{1} (Q^2,r')\rrbracket_{0}\, 
\llbracket C^{\rm B1}_{1}(Q^2,r) \rrbracket_{0} \, \llbracket \widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)\rrbracket_0
\bigg]\,,
\label{eq:Btype_subtracted2}
\end{eqnarray}
where now the integrals over $r,r'$ are evaluated from $0$ to $\infty$. The convolution integrals are now convergent, but they must be evaluated after the integrands are combined.
Provided $\Lambda\ll Q$, this expression can be simplified to
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{
\scriptsize \begin{array}{l}
$\rm B--type$\\[-0.1cm] $i=i'=1$\end{array}} 
&=& 
\frac{2C_F}{Q^2}\,f(\epsilon) \,\widetilde{\mathcal{J}}_{c}^{(g)}(s_L)\,
\widetilde{S}^{(g)}(s_R,s_L)
\nonumber\\ 
&&\hspace*{-3cm}\times\,\bigg\{\int_0^1 dr dr' \,\big[1-\theta(\Lambda/Q-r)\theta(\Lambda/Q-r')\big]\,
C^{\rm B1*}_{1}(Q^2,r')C^{\rm B1}_{1}(Q^2,r) 
\,\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)
\quad 
\nonumber \\ &&\hspace*{-2.5cm}-
\int_0^\infty dr dr' \,
\big[\theta(r-\Lambda/Q)\theta(\Lambda/Q-r')+
\theta(\Lambda/Q-r)\theta(r'-\Lambda/Q)\big]\,
\nonumber \\ &&\hspace*{-2cm}\times\,
\llbracket C^{\rm B1*}_{1} (Q^2,r')\rrbracket_{0}\, 
\llbracket C^{\rm B1}_{1}(Q^2,r) \rrbracket_{0} \, \llbracket \widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)\rrbracket_0\,
\bigg\}\,,
\label{eq:Btype_subtracted3}
\end{eqnarray}
up to corrections of $\mathcal{O}(\Lambda/Q)$. 
The first integral in the  curly brackets is 
simply \eqref{eq:BtypeLaplace} with integration 
region $r,r^\prime < \Lambda/Q$ removed, see
Fig.~\ref{fig:overlap}, 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 (\ref{eq:BtypeLaplace})  we change integration variables to $\bar{r}, \bar{r}'$ to map the singular point $r=r'=1$ to $\bar{r}=\bar{r}'=0$. The subtraction term is then obtained after mapping $\omega =\bar{ r} Q$ and $\omega'=\bar{r}'Q$ and following the same steps as before we find
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{
\scriptsize \begin{array}{l}
$\rm B--type$\\[-0.1cm] $i=i'=2$\end{array}} 
&=& 
\frac{2C_F}{Q^2}\, f(\epsilon)\,\widetilde{\mathcal{J}}_{c}^{(g)}(s_L)\,\widetilde{S}^{(g)}(s_R,s_L)\, \int_0^\infty d\bar{r} d\bar{r}' \nonumber\\ 
&&\hspace*{-2.5cm}\times\, \bigg[ \,
\theta(1-\bar r)\theta(1-\bar{r}')\,
C^{\rm B1*}_{2}(Q^2,r')C^{\rm B1}_{2}(Q^2,r) 
\,\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)
\quad 
\nonumber \\ &&\hspace*{-2cm}-
\big[1-\theta(\bar r-\Lambda/Q)\theta(\bar{r}'-\Lambda/Q)\big]\,
\nonumber \\ &&\hspace*{-2cm}\times\,
\llbracket C^{\rm B1*}_{2} (Q^2,r')\rrbracket_{1}\, 
\llbracket C^{\rm B1}_{2}(Q^2,r) \rrbracket_{1} \, \llbracket \widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)\rrbracket_1
\bigg]\,,
\label{eq:Btype_subtractedii2}
\end{eqnarray}
which coincides with \eqref{eq:Btype_subtracted2} up to 
the expected substitutions. Similar simplifications 
hold when $\Lambda\ll Q$ is assumed.
The total endpoint-subtracted B-type is the sum of \eqref{eq:Btype_subtracted2}, \eqref{eq:Btype_subtractedii2} and  
the mixed term $i\not=i'$ in (\ref{eq:B1fact1}). 

By construction, the dependence on $\Lambda$ cancels 
in the sum of the A- and B-type contributions. In 
App.~\ref{version1endpintfact} we provide the 
corresponding version of the 
endpoint-rearranged factorization formulas, when 
\eqref{eq:scaleless_integral} is split according to the 
first version mentioned below this equation, that is, 
for $\omega$ {\em and} $\omega'$ smaller than 
$\Lambda$. 


\section{Renormalization-group equations}
\label{sec:consistency}

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 {\em integrands} of the 
A- and B-type term must be separately independent of the 
scale $\mu$. The endpoint-factorization conditions then 
provide further relations between the anomalous dimensions 
in the respective large-$\omega$ and small-$r$ limits. We then 
obtain the solutions to leading-logarithmic (LL) accuracy, 
which includes running-coupling effects.

\subsection{RGEs for the A-type functions}
\label{sec:AtypeRGE}

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 
\eqref{eq:EPfactI} to derive the evolution of the collinear 
function $\mathcal{J}_c(p^2,\omega,\omega')$ for large 
$\omega,\omega'$ from the known RGEs for the LP gluon jet function 
and the $D^{\rm 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. 

\paragraph{Hard function}
The evolution of the LP hard matching coefficient $C^{\rm A0}$ is well known \cite{Becher:2006mr}. It obeys the RGE
\begin{align}
\frac{ d}{ d\ln\mu}C^{\rm A0}\klammer{Q^2,\mu^2}=\left[C_F\gamma_{\text{cusp}}\klammer{\alpha_s}\ln\frac{-Q^2}{\mu^2}+\gamma_{\rm A0}\klammer{\alpha_s}\right]\,C^{\rm A0}\klammer{Q^2,\mu^2},
\label{eq:CA0RGE}
\end{align}
where $-Q^2$ in the logarithm should be read as $-Q^2-i\varepsilon$. The cusp anomalous dimension $\gamma_{\text{cusp}}$ and the hard non-cusp anomalous dimension $\gamma_{A0}$ are given by
\begin{align}
&\gamma_{\text{cusp}}\klammer{\alpha_s}=4\,\frac{\alpha_s}{4\pi}+\mathcal{O}\klammer{\alpha_s^2},\\
\label{eq:A0 non-cusp}&\gamma_{\rm A0}\klammer{\alpha_s}=-6C_F\frac{\alpha_s}{4\pi}+\mathcal{O}\klammer{\alpha_s^2}. 
\end{align}
The one-loop cusp anomalous dimension sums the leading 
logarithms. Non-cusp terms become relevant only at NLL 
accuracy.\footnote{We make an exception for the non-cusp terms 
proportional to $\beta_0$, for example in \eqref{eq:gluonADM} 
below, which determine the scale of the coupling that appears 
in the tree-level $\mathcal{O}(\alpha_s)$ term. See 
Sec.~\ref{sec:resummation} for further discussion.} 
Dropping them, we obtain the LL solution
\begin{align}
\label{eq: hard function general result only cusp}
C_{\text{LL}}^{\rm A0}\klammer{Q^2,\mu^2}=&\exp\left[2C_FS\klammer{\mu_h,\mu}\right]\klammer{\frac{-Q^2}{\mu_h^2}}^{-C_FA_{\gamma_{\text{cusp}}}\klammer{\mu_h,\mu}}C^{\rm A0}\klammer{Q^2,\mu_h^2},
\end{align}
where $\mu_h\sim Q$ denotes the hard initial scale, such that 
$C^{\rm A0}\klammer{Q^2,\mu_h^2}$ free from large logarithms. The functions $S\klammer{\nu,\mu}$ and $A_{\gamma_i}\klammer{\nu,\mu}$ are defined as \cite{Neubert:2004dd}
\begin{align}
\label{eq:Resummation Definition of S}&S\klammer{\nu,\mu}=-\int_{\alpha_s\klammer{\nu}}^{\alpha_s\klammer{\mu}} d\alpha\,\frac{\gamma_{\text{cusp}}\klammer{\alpha}}{\beta\klammer{\alpha}}\int_{\alpha_s\klammer{\nu}}^{\alpha} d\alpha'\frac{1}{\beta\klammer{\alpha'}}\,,\\
\label{eq:Resummation Definition of A}&A_{\gamma_i}\klammer{\nu,\mu}=-\int_{\alpha_s\klammer{\nu}}^{\alpha_s\klammer{\mu}} d\alpha\,\frac{\gamma_i\klammer{\alpha}}{\beta\klammer{\alpha}}\,.
\end{align}
The QCD beta-function is
\begin{equation}
\beta(\alpha_s) = \frac{ d\alpha_s}{ d\ln\mu}  
= -2\beta_0\,\frac{\alpha_s^2}{4\pi}
+\mathcal{O}\klammer{\alpha_s^3}, \qquad 
\beta_0 = 11-\frac{2n_l}{3}\,.
\end{equation}

\paragraph{Anti-collinear function}

The anti-collinear function appears already at LP. In Laplace space, it obeys the local RGE \cite{Becher:2006mr}
\begin{align}
\frac{d}{d\ln\mu}\widetilde{\mathcal{J}}_{\bar{c}}^{\klammer{\bar{q}}}\klammer{s,Q^2,\mu^2}=-\left[2C_F\gamma_{\text{cusp}}\klammer{\alpha_s}\ln\frac{Q}{se^{\gamma_E}\mu^2}+\gamma_{\mathcal{J}^{\klammer{\bar{q}}}}\klammer{\alpha_s}\right]\widetilde{\mathcal{J}}_{\bar{c}}^{\klammer{\bar{q}}}\klammer{s,Q^2,\mu^2},
\end{align}
where the non-cusp term is given by
\begin{align}
\gamma_{\mathcal{J}^{\klammer{\bar{q}}}}\klammer{\alpha_s}=-6C_F\frac{\alpha_s}{4\pi}+\mathcal{O}\klammer{\alpha_s^2}.
\end{align}
At LL accuracy, the RGE is solved by
\begin{align}
\widetilde{\mathcal{J}}_{\bar{c}, \text{LL}}^{\klammer{\bar{q}}}\klammer{s,Q^2,\mu^2}=\exp\left[-4C_FS\klammer{\mu_{\bar{c}},\mu}\right]\klammer{\frac{Q}{se^{\gamma_E}\mu_{\bar{c}}^2}}^{2C_FA_{\gamma_{\text{cusp}}}\klammer{\mu_{\bar{c}},\mu}}\widetilde{\mathcal{J}}_{\bar{c}}^{\klammer{\bar{q}}}\klammer{s,Q^2,\mu_{\bar{c}}^2},
\end{align}
where $\mu_{\bar{c}}$ denotes the anti-collinear initial scale.

\paragraph{Collinear function}

The collinear function  $\mathcal{J}_c(p^2,\omega,\omega')$ is a new NLP object and its anomalous dimension is currently unknown. In the asymptotic region of large $\omega, \omega'$, it factorizes into the LP gluon jet function and two $D^{B1}$ matching coefficients as shown in \eqref{eq:EPfactI}. We therefore derive its asymptotic evolution from the LP gluon jet function,  which renormalizes in Laplace space as \cite{Becher:2009th}
\begin{align}
\frac{d}{d\ln\mu}\widetilde{\mathcal{J}}_c^{(g)}(s,Q^2,\mu^2)=
-\left[2C_A\gamma_{\text{cusp}}\klammer{\alpha_s}\ln\frac{Q}{se^{\gamma_E}\mu^2}+\gamma_{\mathcal{J}^{\klammer{g}}}\klammer{\alpha_s}\right]\,\widetilde{\mathcal{J}}_{c}^{(g)}(s,Q^2,\mu^2),
\label{eq:gluonjetRGE}
\end{align}
with the  gluon non-cusp term given by
\begin{align}
\gamma_{\mathcal{J}^{\klammer{g}}}\klammer{\alpha_s}=-2\beta_0\frac{\alpha_s}{4\pi}+\mathcal{O}\klammer{\alpha_s^2}, 
\label{eq:gluonADM}
\end{align}
and the evolution of the $D^{\rm B1}$ coefficient. Making the cusp 
term explicit, \eqref{eq:DBQrge}, \eqref{eq:DB1andim} or 
\eqref{eq:DB1RGEinomega}, \eqref{eq:DB1RGEinomegagamma} 
imply the non-local RGE 
\begin{eqnarray}
\frac{ d}{ d\ln\mu}D^{\rm B1}(\omega,\mu^2)
&=& -\klammer{C_F-C_A}\gamma_{\text{cusp}}\klammer{\alpha_s}\ln\frac{-\omega Q-i\varepsilon}{\mu^2} \,D^{\rm B1}(\omega,\mu^2)
\nonumber\\
&& +\,\int_0^\infty d\hat{\omega}\,\hat{\gamma}_D\klammer{\hat{\omega},\omega} D^{\rm B1}(\hat{\omega},\mu^2),
\label{eq:DB1evolution}
\end{eqnarray}
where the non-cusp anomalous dimension $\hat{\gamma}_D\klammer{\hat{\omega},\omega}$ is given by the second line of \eqref{eq:DB1RGEinomegagamma}. 
The first endpoint factorization consistency condition \eqref{eq:EPfactI} now implies the following asymptotic RGE for the collinear function $\mathcal{J}_c$,
\begin{eqnarray}
\frac{d}{d\ln\mu}\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s,\omega,\omega',Q^2,\mu^2}\right\rrbracket 
&=&-\bigg[2C_A\gamma_{\text{cusp}}\klammer{\alpha_s}\ln\frac{Q}{se^{\gamma_E}\mu^2}
\nonumber\\
&& \hspace*{-3cm} +\,\klammer{C_F-C_A}\gamma_{\text{cusp}}\klammer{\alpha_s}\klammer{\ln\frac{\omega Q}{\mu^2}+\ln\frac{\omega' Q}{\mu^2}}\bigg]\,
\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s,\omega,\omega',Q^2,\mu^2}\right\rrbracket
\nonumber\\
&& \hspace*{-4cm}-\gamma_{\mathcal{J}^{\klammer{g}}}\klammer{\alpha_s}\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s,\omega,\omega',Q^2,\mu^2}\right\rrbracket
+\int_0^\infty d\hat{\omega}\,\frac{\hat{\omega}}{\omega}\,\hat{\gamma}_D\klammer{\hat{\omega},\omega}\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s,\hat{\omega},\omega',Q^2,\mu^2}\right\rrbracket
\nonumber\\
&& \hspace*{-4cm}+\,\int_0^\infty d\hat{\omega}'\,\frac{\hat{\omega}'}{\omega'}\,\hat{\gamma}^*_D\klammer{\hat{\omega}',\omega'}\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s,\omega,\hat{\omega}',Q^2,\mu^2}\right\rrbracket.
\label{eq:Atype_collinearRGE}
\end{eqnarray}
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 $\beta_0$) and obtain the solution
\begin{align}
\left\llbracket\widetilde{\mathcal{J}}_{c,\text{LL}}\klammer{s,\omega,\omega',Q^2,\mu^2}\right\rrbracket=
\nonumber&\frac{\alpha_s\klammer{\mu_c}}{\alpha_s\klammer{\mu}}\exp\left[-4C_AS\klammer{\mu_c,\mu}-4\klammer{C_F-C_A}S\klammer{\mu_{c\Lambda},\mu}\right]\\
\nonumber&\hspace*{-3cm}\times\klammer{\frac{Q}{se^{\gamma_E}\mu_c^2}}^{2C_A A_{\gamma_{\text{cusp}}}\klammer{\mu_c,\mu}}\klammer{\frac{\omega Q}{\mu_{c\Lambda}^2}}^{\klammer{C_F-C_A} A_{\gamma_{\text{cusp}}}\klammer{\mu_{c\Lambda},\mu}}\\
&\hspace*{-3cm}\times\klammer{\frac{\omega' Q}{\mu_{c\Lambda}^2}}^{\klammer{C_F-C_A} A_{\gamma_{\text{cusp}}}\klammer{\mu_{c\Lambda},\mu}}\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s,\omega,\omega',Q^2,\mu_c^2,\mu_{c\Lambda}^2}\right\rrbracket\,.
\end{align}
We introduce two initial scales, $\mu_c$ and $\mu_{c\Lambda}$, 
connected to the two different cusp logarithms in 
\eqref{eq:Atype_collinearRGE}, allowing for the possibility to choose 
different initial scales for the gluon jet function and 
$D^{\rm B1}$ coefficient, which will be used in Sec.~\ref{sec:logsandscales}.


\paragraph{Soft function}

The soft quark function $S_{\rm NLP}(l_+,l_-,\omega,\omega')$ 
appears first at NLP and its complete RGE is currently unknown 
even at lowest order. We can obtain its asymptotic evolution 
for large $\omega, \omega'$ from RGE consistency. For 
fixed values of $\omega, \omega'$ 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,
\begin{eqnarray}
\frac{ d}{ d\ln\mu}\left\llbracket\widetilde{S}_{\text{NLP}}(s_R,s_L,\omega,\omega',\mu^2)\right\rrbracket
&=&\bigg[2C_A\gamma_{\text{cusp}}(\alpha_s)\ln\frac{1}{s_Le^{\gamma_E}s_Re^{\gamma_E}\mu^2}
\nonumber\\
&& \hspace*{-5cm} 
+\,\klammer{C_F-C_A}\gamma_{\text{cusp}}(\alpha_s)\klammer{\ln\frac{\omega }{s_Re^{\gamma_E}\mu^2}+\ln\frac{\omega'}{s_Re^{\gamma_E}\mu^2}}\bigg] 
\left\llbracket\widetilde{S}_{\text{NLP}}(s_R,s_L,\omega,\omega',\mu^2)\right\rrbracket
\nonumber\\[0.1cm]
&& \hspace*{-5cm}+\,\Big[\gamma_{\mathcal{J}^{\klammer{\bar{q}}}}(\alpha_s)+\gamma_{\mathcal{J}^{(g)}}(\alpha_s)-2\gamma_{A0}\klammer{\alpha_s}\Big]\left\llbracket\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega',\mu^2}\right\rrbracket
\nonumber\\
&& \hspace*{-5cm}-\,\int_0^\infty d\hat{\omega}\,\frac{\omega}{\hat{\omega}}\,\hat{\gamma}_D(\omega,\hat{\omega})\left\llbracket\widetilde{S}_{\text{NLP}}(s_R,s_L,\hat{\omega},\omega',\mu^2)\right\rrbracket
\nonumber\\
&& \hspace*{-5cm}-\,\int_0^\infty d\hat{\omega}'\,\frac{\omega'}{\hat{\omega}'}\,\hat{\gamma}^*_D(\omega',\hat{\omega}')\left\llbracket\widetilde{S}_{\text{NLP}}(s_R,s_L,\omega,\hat{\omega}',\mu^2)\right\rrbracket .
\end{eqnarray}
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
\begin{eqnarray}
\left\llbracket\widetilde{S}_{\text{NLP,LL}}(s_R,s_L,\omega,\omega',\mu^2)\right\rrbracket
&=&\frac{\alpha_{s}(\mu)}{\alpha_s(\mu_s)}\exp\left[4C_AS\klammer{\mu_s,\mu}+4\klammer{C_F-C_A}S\klammer{\mu_{s\Lambda},\mu}\right]
\nonumber\\
&&\hspace*{-4cm}\times\,
\klammer{\frac{1}{s_Le^{\gamma_E}s_Re^{\gamma_E}\mu_s^2}}^{\!-2C_AA_{\gamma_{\text{cusp}}}\klammer{\mu_s,\mu}}
\,\klammer{\frac{\omega }{s_Re^{\gamma_E}\mu_{s\Lambda}^2}}^{\!-\klammer{C_F-C_A}A_{\gamma_{\text{cusp}}}\klammer{\mu_{s\Lambda},\mu}}
\nonumber\\
&&\hspace*{-4cm}\times\,\klammer{\frac{\omega' }{s_Re^{\gamma_E}\mu_{s\Lambda}^2}}^{\!-\klammer{C_F-C_A}A_{\gamma_{\text{cusp}}}\klammer{\mu_{s\Lambda},\mu}}\,
\left\llbracket\widetilde{S}_{\text{NLP}}(s_R,s_L,\omega,\omega',\mu_s^2,\mu_{s\Lambda}^2)\right\rrbracket,
\end{eqnarray}
where again we allow for two separate initial soft scales, 
$\mu_s$ and $\mu_{s\Lambda}$. The lowest-order initial condition is proportional to  $\alpha_s$, which we choose to evaluate at $\mu_s$. Together with the prefactor $\alpha_s(\mu)/\alpha_s(\mu_s)$, which arises from the $\beta_0$-term in $\gamma_{\mathcal{J}^{(g)}}$, this results in an overall factor of $\alpha_s(\mu)$. 
The point of keeping the non-cusp $\beta_0$-term in $\gamma_{\mathcal{J}^{(g)}}$ is that it renders the result independent of the 
initial scale -- choosing to evaluate $\alpha_s$ at $\mu_{s\Lambda}$ instead, we would obtain the same resummed soft function.   

\subsection{RGEs for the B-type functions}
\label{sec:BtypeRGE}

In the B-type contribution, the collinear function and the soft function are known LP objects while the hard functions $C_i^{\rm B1}$ and the anti-collinear function $\mathcal{J}_{\bar c}^{q\bar q(8)}(p^2,r,r')$ are NLP objects. The evolution equation of the hard matching coefficients is of the form discussed in \cite{Beneke:2017ztn,Beneke:2018rbh} 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'\to 0, 1$, 
and therefore we will only need the equations for the 
asymptotic hard and anti-collinear functions. 

\paragraph{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 
\eqref{eq:gluonjetRGE} and \eqref{eq:gluonADM}, respectively.
The LL solution to the RGE reads
\begin{align}
\widetilde{\mathcal{J}}_{c,\text{LL}}^{\klammer{g}}\klammer{s,\mu^2}=\frac{\alpha_s\klammer{\mu_c}}{\alpha_s\klammer{\mu}}\exp\left[-4C_AS\klammer{\mu_c,\mu}\right]\klammer{\frac{Q}{se^{\gamma_E}\mu_c^2}}^{2C_AA_{\gamma_{\text{cusp}}}\klammer{\mu_c,\mu}}\widetilde{\mathcal{J}}_c^{\klammer{g}}\klammer{s,\mu_c^2},
\end{align}
where $\mu_c$ denotes the collinear initial scale. 

\paragraph{Soft function}

The LP two-hemisphere soft function for gluons, which is defined 
with adjoint Wilson lines, obeys the evolution equation
\begin{eqnarray}
\label{eq:LPGluonSoftFunctionADM}
\frac{d}{d\ln\mu}\widetilde{S}^{(g)}(s_R,s_L,\mu^2)=\left[2C_A\gamma_{\text{cusp}}(\alpha_s)\ln\frac{1}{s_L s_Re^{2\gamma_E}\mu^2}+\gamma_{S^{(g)}}(\alpha_s)\right]\widetilde{S}^{\klammer{g}}(s_R,s_L,\mu^2)\;\;
\end{eqnarray}
in Laplace space. 
The non-cusp part of the anomalous dimension is given by \cite{Berger:2010xi}
\begin{align}
\gamma_{S^{(g)}}(\alpha_s)=0+\mathcal{O}(\alpha_s^2)\,.
\end{align}
At LL accuracy, the solution to the RGE reads
\begin{align}
\widetilde{S}_{\text{LL}}^{\klammer{g}}(s_R,s_L,\mu^2)=&\,
\exp\left[4C_AS\klammer{\mu_s,\mu}\right]\klammer{\frac{1}{s_L s_Re^{2\gamma_E}\mu_s^2}}^{\!-2C_AA_{\gamma_{\text{cusp}}}\klammer{\mu_s,\mu}}\\
\nonumber&\times\widetilde{S}^{\klammer{g}}(s_R,s_L,\mu_s^2),
\end{align}
where $\mu_s$ denotes the soft initial scale.

\paragraph{Hard function} 
Following the calculations of \cite{Beneke:2017ztn,Beneke:2018rbh}, the specific anomalous dimension matrix for 
the two coefficients $C^{\rm B1}_i$, $i=1,2$
needed here can be inferred from App.~\ref{app:asympDB1RGE}. For 
the flavour non-singlet case, which is assumed here, $\Delta_F$ 
can be set to zero there, and the two coefficients evolve 
independently,  
\begin{eqnarray}
\frac{ d}{ d\ln\mu}C^{\rm B1}_i(Q^2,r) &=& \int_0^1  d\hat r\, \gamma_{\rm B1}(\hat r,r)C^{\rm B1}_i(Q^2,\hat r)\,,
\end{eqnarray}
with the same anomalous dimension (see~\eqref{eq:gamma}) 
\begin{eqnarray}
\gamma_{\rm B1}(r,s) &=& -\frac{\alpha_s}{\pi}\,\delta(r-s)\left[C_A\ln\left(\frac{\mu^2}{-Q^2}\right)-\frac{C_A}{2}\ln(r\bar r)+\frac{3C_F}{2}\right]+2\gamma(r,s)\,,\qquad
\end{eqnarray}
where
\begin{eqnarray}
\gamma(r,s) &=& \frac{\alpha_s}{2\pi}\left(\frac{C_A}{2}-C_F\right)\Bigg[ \frac{\bar r}{\bar s}\left(\left(\frac{\theta(r-s)}{r-s}\right)_+ + \theta(r-s)\right)\nonumber\\
&& {} + \frac{r}{ s}\left(\left(\frac{\theta(s-r)}{s-r}\right)_+ + \theta(s-r)\right) \Bigg]\,.
\end{eqnarray} 

The non-cusp part of this RGE contains a hidden large 
logarithm of $r$ ($\bar r$), which must be extracted and 
combined with the cusp part when one wants to sum the hard 
function for small $r$ or $\bar{r}$, see 
App.~\ref{app:asympDB1RGE}. This 
results in the evolution equation \eqref{eq:DBQrge} for 
the $D^{\rm B1}$ coefficient. The evolution equation 
for the asymptotic $C^{\rm B1}_i$ coefficients then 
follows from \eqref{eq:B1fact1}, \eqref{eq:B1fact2}. 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^{\rm B1}$ coefficient. Combining 
\eqref{eq:CA0RGE} and \eqref{eq:DB1evolution}, we find 
\begin{eqnarray}
\frac{ d}{ d\ln\mu}\left\llbracket C^{\rm B1}_1\klammer{Q^2,r,\mu^2}\right\rrbracket_0
&=&\bigg[\,C_F\gamma_{\text{cusp}}\klammer{\alpha_s}\ln\frac{-Q^2}{\mu^2}
\nonumber\\
&&\hspace*{-3cm}-\,\klammer{C_F-C_A}\gamma_{\text{cusp}}\klammer{\alpha_s}\ln\frac{-rQ^2}{\mu^2}+\gamma_{\rm A0}\klammer{\alpha_s}\bigg]\left\llbracket C^{\rm B1}_1\klammer{Q^2,r,\mu^2}\right\rrbracket_0
\nonumber\\
&&\hspace*{-3cm}+\,Q\int_{0}^\infty d\hat{r}\,\frac{\hat{r}}{r}\,
\hat{\gamma}_D\klammer{\hat{r}Q,rQ}\left\llbracket C^{\rm B1}_1\klammer{Q^2,\hat{r},\mu^2}\right\rrbracket_0\,.
\end{eqnarray}
The LL solution reads 
\begin{eqnarray}
\left\llbracket C_1^{\rm B1}\klammer{Q^2,r,\mu^2}\right\rrbracket_0
&=&\exp\left[2C_FS\klammer{\mu_h,\mu}-2\klammer{C_F-C_A}S\klammer{\mu_{h\Lambda},\mu}\right]
\nonumber\\
&&\hspace*{-2cm}\times
\left(\frac{-Q^2}{\mu_h^2}\right)^{\!-C_FA_{\gamma_{\text{cusp}}}\klammer{\mu_h,\mu}}\klammer{\frac{-rQ^2}{\mu_{h\Lambda}^2}}^{\klammer{C_F-C_A}A_{\gamma_{\text{cusp}}}\klammer{\mu_{h\Lambda},\mu}}
\nonumber\\
&& \hspace*{-2cm}\times\left\llbracket C_1^{\rm B1}\klammer{Q^2,r,\mu_h^2,\mu_{h\Lambda}^2}\right\rrbracket_0,\quad
\label{eq:CB1LLsol}
\end{eqnarray}
where as in the case of similar A-type RGE solutions, we allowed for two different initial scales, $\mu_h$ and $\mu_{h\Lambda}$, 
connected with the different cusp logarithms. 

For each of the two endpoint limits $r\rightarrow0,1$, we only 
need one of the two B1 coefficients since the other one is regular. 
The RGE equation and solution for 
$C_2^{\rm B1}\klammer{Q^2,r,\mu^2}$ is obtained from the above 
by replacing $r\to \bar{r}$, and similarly for $\hat{r}$ including 
the integration measure $ d\hat{r}$.

\paragraph{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 RGE-invariant on its own, which allows us to derive the RGE and anomalous dimension of the anti-collinear 
function in the form
\begin{eqnarray}
\frac{ d}{ d\ln\mu}\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s_R,r,r'}&=&\bigg[-2C_A\gamma_{\text{cusp}}\klammer{\alpha_s}\ln\frac{Q}{s_Re^{\gamma_E}\mu^2}
\nonumber\\
&& \hspace*{-3cm}
+\,\gamma_{\mathcal{J}^{\klammer{g}}}(\alpha_s)-\gamma_{S^{\klammer{g}}}\klammer{\alpha_s}\bigg]\,\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s_R,r,r'}\nonumber\\
&&\hspace*{-3cm}
 -\,\int_0^1 d\hat{r}\,\hat{\gamma}_{\rm B1}\klammer{r,\hat{r}}\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s_R,\hat{r},r'}-\int_0^1 d\hat{r}'\,\hat{\gamma}^*_{\rm B1}\klammer{r',\hat{r}'    }\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s_R,r,\hat{r}'}\,,\qquad
\end{eqnarray}
where the non-cusp part of the B1 anomalous dimension $\hat{\gamma}_{\rm B1}$ is given by
\begin{align}
\hat{\gamma}_{\rm B1}\klammer{r,s}=\frac{\alpha_s}{\pi}\delta\klammer{r-s}\left[\frac{C_A}{2}\ln\klammer{r\bar{r}}-\frac{3C_F}{2}\right]+2\gamma\klammer{r,s}\,.
\end{align}
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
\begin{eqnarray}
\frac{ d}{ d\ln\mu}\left\llbracket\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s,r,r',\mu^2}\right\rrbracket_0&=&-\bigg[2C_F\gamma_{\text{cusp}}\klammer{\alpha_s}\ln\frac{Q}{se^{\gamma_E}\mu^2}\nonumber\\
&&\hspace*{-3cm}-\klammer{C_F-C_A}\gamma_{\text{cusp}}\klammer{\alpha_s}\left(\ln\frac{rQ}{se^{\gamma_E}\mu^2}+\ln\frac{r'Q}{se^{\gamma_E}\mu^2}\right)\nonumber\\
&&\hspace*{-3cm}+\left[2\gamma_{A0}(\alpha_s)+\gamma_{{S}^{\klammer{g}}}(\alpha_s)-\gamma_{\mathcal{J}^{\klammer{g}}}(\alpha_s)\right]\bigg]\left\llbracket\widetilde{\mathcal{J}}^{q\bar{q}\klammer{8}}_{\bar{c}}\klammer{s,r,r',\mu^2}\right\rrbracket_0\nonumber\\
&&\hspace*{-3cm}-Q\int_0^\infty d\hat{r}\,\frac{r}{\hat{r}}\,\hat{\gamma}_D\klammer{rQ,\hat{r}Q}\left\llbracket\widetilde{\mathcal{J}}^{q\bar{q}\klammer{8}}_{\bar{c}}\klammer{s,\hat{r},r',\mu^2}\right\rrbracket_0\nonumber\\
&&\hspace*{-3cm}-Q\int_0^\infty d\hat{r}'\,\frac{r'}{\hat{r}'}\,\hat{\gamma}^*_D\klammer{r'Q,\hat{r}'Q}\left\llbracket\widetilde{\mathcal{J}}^{q\bar{q}\klammer{8}}_{\bar{c}}\klammer{s,r,\hat{r}',\mu^2}\right\rrbracket_0\,.
\end{eqnarray}
At LL accuracy, the asymptotic RGE is solved by
\begin{eqnarray}
\left\llbracket\widetilde{\mathcal{J}}_{\bar{c},\text{LL}}^{q\bar{q}\klammer{8}}\klammer{s,r,r',\mu^2}\right\rrbracket_0&=&\frac{\alpha_s(\mu)}{\alpha_s(\mu_{\bar{c}})}
\exp\left[-4C_FS\klammer{\mu_{\bar{c}},\mu}+4\klammer{C_F-C_A}S\klammer{\mu_{\bar{c}\Lambda},\mu}\right]\nonumber\\
&&\hspace*{-3cm}\times\,\klammer{\frac{Q}{se^{\gamma_E}\mu_{\bar{c}}^2}}^{\!2C_FA_{\gamma_{\text{cusp}}}\klammer{\mu_{\bar{c}},\mu}}\klammer{\frac{rQ}{se^{\gamma_E}\mu_{\bar{c}\Lambda}^2}}^{\!-\klammer{C_F-C_A}A_{\gamma_{\text{cusp}}}\klammer{\mu_{\bar{c}\Lambda},\mu}}\nonumber\\
&&\hspace*{-3cm}
\times\,\klammer{\frac{r'Q}{se^{\gamma_E}\mu_{\bar{c}\Lambda}^2}}^{\!-\klammer{C_F-C_A}A_{\gamma_{\text{cusp}}}\klammer{\mu_{\bar{c}\Lambda},\mu}}\left\llbracket\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s,r,r',\mu_{\bar{c}}^2,\mu_{\bar{c}\Lambda}^2}\right\rrbracket_0,
\end{eqnarray}
where $\mu_{\bar{c}}$ and $\mu_{\bar{c}\Lambda}$ denote the 
anti-collinear initial scales. As for the A-type soft function, we may evaluate the factor of $\alpha_s$ in the fixed-order initial condition either at $\mu_{\bar{c}}$ or $\mu_{\bar{c}\Lambda}$, 
and $\mu_{\bar{c}}$ has been chosen above. This choice does not affect the overall resummed anti-collinear function since the ratio of couplings in front always guarantees a global factor of $\alpha_s\klammer{\mu}$. Once again the RGE and solution 
for the anti-collinear function for $r\to 1$, $\llbracket\widetilde{\mathcal{J}}_{\bar{c},\text{LL}}^{q\bar{q}\klammer{8}}\klammer{s,r,r',\mu^2}\rrbracket_1$, are obtained from the above by replacing 
 $r\to \bar{r}$, and similarly for $\hat{r}$ including 
the integration measure $ d\hat{r}$ as well as for the 
corresponding primed quantities.


\section{Resummation} 
\label{sec: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
\eqref{eq:Atype_subtracted_general2} and 
\eqref{eq:Btype_subtracted2} for the A- and B-type term, 
respectively. In \eqref{eq:Atype_subtracted_general2}, 
we rewrite
\begin{eqnarray}
&&\widetilde{\mathcal{J}}_c\klammer{s_L,\omega,\omega'}\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega'}
\nonumber\\
&&\hspace*{0.5cm}=\,
\widetilde{\mathcal{J}}_c\klammer{s_L,\omega,\omega'}\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega'}
-\sigma(\omega,\omega')\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s_L,\omega,\omega'}\right\rrbracket\left\llbracket\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega'}\right\rrbracket
\nonumber\\
&&\hspace*{1cm}+\,
\sigma(\omega,\omega')\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s_L,\omega,\omega'}\right\rrbracket\left\llbracket\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega'}\right\rrbracket,
\end{eqnarray}
where $\sigma(\omega,\omega')$ is an auxiliary function, 
which equals unity whenever $\omega^{(\prime)}s_R, 
\omega^{(\prime)}s_L \gg 1$ and which goes to zero 
for  $\omega^{(\prime)}\to 0$. We find 
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds_R ds_L}|_{\rm A-type}&=&
\frac{2C_F}{Q}\,\left|C^{\rm A0}\klammer{Q^2}\right|^2\widetilde{\mathcal{J}}^{\klammer{\bar{q}}}_{\bar{c}}\klammer{s_R}
\int_0^\infty d\omega d\omega'
\times\Bigg\{
\nonumber\\
&&\hspace*{-3cm}
\left[\sigma(\omega,\omega')-\theta\klammer{\omega-\Lambda}\theta\klammer{\omega'-\Lambda}\right]\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s_L,\omega,\omega'}\right\rrbracket\left\llbracket\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega'}\right\rrbracket
\nonumber\\
&&\hspace*{-3cm}+\,\bigg[\widetilde{\mathcal{J}}_c\klammer{s_L,\omega,\omega'}\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega'}
-\sigma(\omega,\omega')\left\llbracket\widetilde{\mathcal{J}}_c\klammer{s_L,\omega,\omega'}\right\rrbracket\!\left\llbracket\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega'}\right\rrbracket\bigg]
\nonumber\\
&&\hspace*{-3cm}+\,\widetilde{\widehat{\mathcal{J}}}_c\klammer{s_L,\omega,\omega'}\widetilde{\widehat{S}}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega'}\Bigg\}\,.
\label{eq:Atyperearranged}
\end{eqnarray}
In this expression, the first line in curly brackets is the 
subtraction term, which contains the endpoint divergence, 
if $\Lambda\to\infty$, and therefore results in a 
logarithmically enhanced contribution from the 
$\omega,\omega'$ convolution. The next term in square 
brackets resembles a ``plus-distribution'' type term -- 
this combination has no support for large $\omega,\omega'$ 
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 $\sigma$ 
is necessary to eliminate a spurious singularity 
for small $\omega,\omega'$, which arises if one wants 
to integrate the subtraction and plus-distribution-like 
term separately, and which 
cancels between the two.

Similarly, we write \eqref{eq:Btype_subtracted2} in the form
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds_R ds_L}|_{
\scriptsize \begin{array}{l}
$\rm B--type$\\[-0.1cm] $i=i'=1$\end{array}} &=&
\frac{2C_F}{Q^2}\,\widetilde{J}_c^{\klammer{g}}\klammer{s_L}\widetilde{S}^{\klammer{g}}\klammer{s_R,s_L}\times\Bigg\{
\nonumber\\
&&\hspace*{-2cm}\int_0^\infty dr dr'\,\left[\theta\klammer{1-r}\theta\klammer{1-r'}-1+\theta\klammer{r-\frac{\Lambda}{Q}}\theta\klammer{r'-\frac{\Lambda}{Q}}\right]
\nonumber\\
&&\hspace*{-1.5cm}\times\left\llbracket C_1^{\rm B1*}\klammer{Q^2,r'}\right\rrbracket_0 \left\llbracket C_1^{\rm B1}\klammer{Q^2,r}\right\rrbracket_0 \left\llbracket \widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s_R,r,r'}\right\rrbracket_0
\nonumber\\
&&\hspace*{-2cm}+
\int_0^1 dr dr'\,\bigg[C_1^{\rm B1*}\klammer{Q^2,r'}C_1^{\rm B1}\klammer{Q^2,r}\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s_R,r,r'}
\nonumber\\[-0.15cm]
&&\hspace*{-1.5cm}-\left\llbracket C_1^{\rm B1*}\klammer{Q^2,r'}\right\rrbracket_0 \left\llbracket C_1^{\rm B1}\klammer{Q^2,r}\right\rrbracket_0 \left\llbracket \widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s_R,r,r'}\right\rrbracket_0\bigg]\,\Bigg\}\,.
\label{eq:Btyperearranged}
\end{eqnarray}
The first integral in curly brackets represents the 
subtraction term, which contributes to the leading 
logarithms and cancels the $\Lambda$-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^{(\prime)}
\leftrightarrow \bar{r}^{(\prime)}$. The $i\not= 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 non-asymptotic 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-to-leading logarithms 
associated with the running coupling in the leading 
terms.

\subsection{Counting of logarithms and proper 
choice of scales}
\label{sec:logsandscales}

Standard Sudakov exponentiation puts an observable into 
the form 
\begin{equation}
I \sim F(\alpha_s) \exp\left(\ln\tau \,g_0(\lambda) + 
g_1(\lambda) +\alpha_s g_2(\lambda)+\ldots\right)\,,
\end{equation}
where $\lambda = \alpha_s\ln\tau$ counts as 
$\mathcal{O}(1)$, while $\tau\ll 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 
$\mathcal{O}(\alpha_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 {\em exponent}.

The {\em integrands} of \eqref{eq:Atyperearranged}, 
\eqref{eq:Btyperearranged} are of this form when the 
solutions to the RGEs from Sec.~\ref{sec:consistency} 
are used. However, one now has to take integrals of 
the form
\begin{equation}
\int_{\tau Q}^\Lambda \frac{d\omega}{\omega}\,I(\omega),
\quad 
\int^1_{\Lambda/Q} \frac{dr}{r}\,I(r),
\label{eq:logints}
\end{equation}
and the assumption underlying standard Sudakov resummation 
that $\omega =\mathcal{O}(\tau)$ and $r=\mathcal{O}(1)$ 
is not justified, since either $\Lambda\gg\tau Q$ 
or $\Lambda\ll Q$, or both. The fact that the integrals 
are logarithmic for large $\omega$ and small $r$ implies 
two modifications of the standard resummation scheme.

First, the logarithmic integrals combine to an extra 
logarithm promoting their log-counting by one order. 
Thus, inserting the LL-resummed RGE functions in the 
subtraction terms in \eqref{eq:Atyperearranged}, 
\eqref{eq:Btyperearranged} 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 
$\omega$ 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  \eqref{eq:Atyperearranged}, 
\eqref{eq:Btyperearranged}.

\begin{figure}
\begin{center}
\includegraphics[width=0.65\textwidth]{figs/initialconditionresummation}
\caption{
Graphical illustration of the resummation of the initial condition 
of the B1 operator for $r\ll 1$ (left) and the LL-accurate 
approximate procedure employed here.
\label{fig:ICres}}
\end{center}
\end{figure}

Second, the fact that $\omega =\mathcal{O}(\tau)$ 
and $r=\mathcal{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^{\rm B1}_1(Q^2,r,\mu^2)$. Standard RGE summation 
implies that in the initial condition, one chooses 
the scale $\mu_h\sim Q$, at which there are no large 
logarithms and evolves to $\mu\ll Q$. However, the 
initial condition for $\mu_h\sim Q$ may still contain 
an infinite series of $\ln r$ terms. While these 
normally count as $\mathcal{O}(1)$, they lead to a 
breakdown of the perturbative expansion in integrals 
such as \eqref{eq:logints}, which receive contributions 
from $r\ll 1$. These logarithms are summed by evolving 
with asymptotic evolution kernel for the $D^{\rm B1}$ 
coefficient from $\mu^2\sim r Q^2$ to $\mu_h^2$. At LL 
accuracy, 
combining this evolution with the standard evolution from 
$\mu_h$ to $\mu$, amounts to the expression 
\eqref{eq:CB1LLsol} with split initial scales -- 
$\mu_h\sim Q$ in $C^{\rm A0}$, and the running scale 
$\mu_{h\Lambda}^2\sim r Q^2$ for the part that arises 
from the $D^{\rm B1}$ coefficient, whose natural 
scale is $r Q^2$. The procedure is illustrated in 
Fig.~\ref{fig:ICres}. A similar reasoning applies 
to the choice of initial conditions for the anti-collinear 
function $\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s_R,r,r'}$ and the functions 
$\widetilde{\mathcal{J}}_c\klammer{s_L,\omega,\omega',Q^2,\mu^2}$, 
$\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega',\mu^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 
\begin{align}
&\mu_h^2\sim Q^2\,,&&\mu_c^2\sim \frac{Q}{s_L}\,,&&\mu_{\bar{c}}^2\sim \frac{Q}{s_R}\,,&&\mu_s^2\sim \frac{1}{s_Ls_R}\,,
\nonumber\\
&\mu_{h\Lambda}^2\sim rQ^2\,,&&\mu_{c\Lambda}^2\sim \omega Q\,,&&\mu_{\bar{c}\Lambda}^2\sim \frac{rQ}{s_R}\,,&&\mu_{s\Lambda}^2\sim \frac{\omega}{s_R}\,.
\label{eq:initialscales}
\end{align}
To ensure the cancellation of the rearrangement scale 
$\Lambda$ between the A- and B-type term, the integrands 
of the subtractions terms must match, which requires
\begin{equation}
\label{eq:intermediatescaleconstraints}
\mu_{s\Lambda}\sim\mu_{\bar{c}\Lambda}\quad\text{and}
\quad\mu_{h\Lambda}\sim\mu_{c\Lambda}\,,
\end{equation}
which is satisfied by \eqref{eq:initialscales} with 
the identification $rQ\leftrightarrow\omega$. 

We now insert the resummed functions from Sec.~\ref{sec:consistency}, 
evolved from their respective initial scales to a common scale $\mu$,  
into the {\em integrands} of \eqref{eq:Atyperearranged},  
\eqref{eq:Btyperearranged}, and simplify the expressions employing 
the identities  
\begin{eqnarray}
\label{eq:ASidentities}
&& A_{\gamma_i}\klammer{\mu_1,\mu_2}+A_{\gamma_i}\klammer{\mu_2,\mu_3}=A_{\gamma_i}\klammer{\mu_1,\mu_3}\,,\\
&&S\klammer{\mu_1,\mu_2}+S\klammer{\mu_2,\mu_3}=S\klammer{\mu_1,\mu_3}+\ln\frac{\mu_1}{\mu_2}\,A_{\gamma_{\text{cusp}}}\klammer{\mu_2,\mu_3}\,,
\end{eqnarray} 
and the special case $\mu_3=\mu_1$. 
We also abbreviate $A\equiv A_{\gamma_{\text{cusp}}}$. 
We find 
\begin{eqnarray}\label{eq:LLintegrandA}
&& \left|C^{\rm A0}(Q^2)\right|^2\widetilde{\mathcal{J}}_{\bar{c}}^{\klammer{\bar{q}}}\klammer{s_R,Q^2}
\left\llbracket\widetilde{\mathcal{J}}_{c}\klammer{s_L,\omega,\omega'}\right\rrbracket\left\llbracket\widetilde{S}_{\text{NLP}}\klammer{s_R,s_L,\omega,\omega'}\right\rrbracket
\nonumber\\[0.1cm]
&&\hspace*{0.5cm} 
=\,\exp\left[4C_FS\klammer{\mu_h,\mu_{\bar{c}}}+4C_AS\klammer{\mu_s,\mu_c}+4\klammer{C_F-C_A}S\klammer{\mu_{s\Lambda},\mu_{c\Lambda}}\right]
\nonumber\\[0.1cm]
&&\hspace*{1cm} 
\times\klammer{\frac{Q^2}{\mu_h^2}}^{\!-2C_FA\klammer{\mu_h,\mu_{\bar{c}}}}\klammer{\frac{1}{s_Ls_Re^{2\gamma_E} \mu_s^2}}^{\!-2C_A A\klammer{\mu_s,\mu_c}}\klammer{\frac{\omega }{s_R e^{\gamma_E}\mu_{s\Lambda}^2}}^{\!-2\klammer{C_F-C_A}A\klammer{\mu_{s\Lambda},\mu_{c\Lambda}}}
\nonumber\\
&&\hspace*{1cm} \times\klammer{s_R e^{\gamma_E} Q}^{2C_FA\klammer{\mu_{c\Lambda},\mu_{\bar{c}}}+2C_A A\klammer{\mu_c,\mu_{c\Lambda}}}\left[\frac{\alpha_s(\mu_{c})}{4\pi}\,\delta\klammer{\omega-\omega'}\frac{1}{\omega s_R}\right]
\end{eqnarray}
for the A-type integrand, and for the B-type, we find
\begin{eqnarray}\label{eq:LLintegrandB}
	&&\left\llbracket C_{1}^{\rm B1}\klammer{Q^2,r}\right\rrbracket_0\left\llbracket C_{1}^{\rm B1*}\klammer{Q^2,r'}\right\rrbracket_0\widetilde{\mathcal{J}}_{c}^{\klammer{g}}\klammer{s_L}\left\llbracket\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\klammer{s_R,r,r'}\right\rrbracket_0\widetilde{S}^{\klammer{g}}\klammer{s_R,s_L}\nonumber\\[0.1cm]
	&&\hspace*{0.5cm} 
	=\,\exp\left[4C_FS\klammer{\mu_h,\mu_{\bar{c}}}
	+4C_AS\klammer{\mu_s,\mu_c}
	+4\klammer{C_F-C_A}S\klammer{\mu_{\bar{c}\Lambda},\mu_{h\Lambda}}\right]
	\nonumber\\[0.1cm]
	&&\hspace*{1cm} \times\klammer{\frac{Q^2}{\mu_h^2}}^{\!-2C_F A\klammer{\mu_h,\mu_{\bar{c}}}}
	\klammer{\frac{1}{s_Ls_R e^{2\gamma_E}\mu_s^2}}^{\!-2C_A A\klammer{\mu_s,\mu_c}}\klammer{\frac{rQ}{s_R e^{\gamma_E}\mu_{\bar c\Lambda}^2}}^{\!-2\klammer{C_F-C_A}A\klammer{\mu_{\bar{c}\Lambda},\mu_{h\Lambda}}}
	\nonumber\\
	&&\hspace*{1cm}\times\klammer{s_R e^{\gamma_E} Q}^{2C_F A\klammer{\mu_{h\Lambda},\mu_{\bar{c}}}+2C_A A\klammer{\mu_c,\mu_{h\Lambda}}}\left[\frac{\alpha_s(\mu_c)}{4\pi}\delta\klammer{r-r'}\frac{Q}{rs_R}\right].
\end{eqnarray}
As it should be with all functions evaluated consistently to LL accuracy the above expressions are manifestly independent of the arbitrary renormalization scale $\mu$. 
Note that if we would not have included the terms proportional to $\beta_0$ in the non-cusp parts of the anomalous dimensions, the only change would be that the overall prefactor of $\alpha_s$ in the square bracket in~\eqref{eq:LLintegrandA}
would be evaluated at scale $\mu_s$ (or $\mu_{s\Lambda}$) instead of $\mu_c$. Similarly, in~\eqref{eq:LLintegrandB} the prefactor would be evaluated at  $\mu_{\bar{c}}$ (or $\mu_{\bar{c}\Lambda}$) instead of $\mu_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 $\beta_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 $\widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}\to \widetilde{\mathcal{J}}_{\bar{c}}^{q\bar{q}\klammer{8}}/\alpha_s$ and $C_{1}^{\rm B1}\to C_{1}^{\rm B1}\times\sqrt{\alpha_s}$, that corresponds to a reshuffling of the overall factor of $\alpha_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 \eqref{eq:intermediatescaleconstraints}, both expressions 
are exactly equal up to an overall factor of $Q$ with the 
identification $\omega^{(\prime)}=r^{(\prime)} Q$.
We thus sum the subtraction terms from
\eqref{eq:Atyperearranged},  \eqref{eq:Btyperearranged}, 
and change variables from $r$ to $\omega/Q$ in the 
B-type contribution \eqref{eq:Btyperearranged}. We also 
choose the auxiliary function to be of the form 
$\sigma(\omega,\omega) = \theta(\omega-\sigma)$, where 
$\sigma$ is of order $1/s_R$, or $1/s_L$. The two integrals 
to be added then combine as 
\begin{equation}
\int_0^\infty\frac{d\omega}{\omega}
\Big[\,\theta(\omega-\sigma)-\theta(\omega-\Lambda)\,\Big] 
+\int_{\Lambda/Q}^1 \frac{dr}{r} 
=\int_\sigma^Q\frac{d\omega}{\omega}\,,
\end{equation}
and we obtain 
\begin{eqnarray}
	\frac{1}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds_R ds_L}|_{\rm LL}&=&
	2\cdot \frac{2 C_F}{Q s_R}\frac{\alpha_s(\mu_c)}{4\pi}\,
	\,\exp\left[4C_FS\klammer{\mu_h,\mu_{\bar{c}}}
	+4C_AS\klammer{\mu_s,\mu_c}\right]
	\nonumber\\
	&&\hspace*{-2cm}
	\times\klammer{\frac{Q^2}{\mu_h^2}}^{\!-2C_F A\klammer{\mu_h,\mu_{\bar{c}}}}
	\klammer{\frac{1}{s_Ls_R e^{2\gamma_E}\mu_s^2}}^{\!-2C_A A\klammer{\mu_s,\mu_c}}
	\nonumber\\
	&&\hspace*{-2cm}
	\times\int_\sigma^Q\frac{d\omega}{\omega}\,
	\exp\left[4\klammer{C_F-C_A}S\klammer{\mu_{s\Lambda},\mu_{h\Lambda}}\right]
	\klammer{\frac{\omega }{s_R e^{\gamma_E} \mu_{s\Lambda}^2}}^{\!-2\klammer{C_F-C_A}A\klammer{\mu_{s\Lambda},\mu_{h\Lambda}}}
	\nonumber\\
	&&\hspace*{-0.8cm}\times\,\klammer{s_R e^{\gamma_E} Q}^{2C_FA\klammer{\mu_{h\Lambda},\mu_{\bar{c}}}+2C_A A\klammer{\mu_c,\mu_{h\Lambda}}}
	\,,
	\label{eq:mainLLresult}
\end{eqnarray}
where $\mu_{h\Lambda}^2\sim \omega Q$ and 
$\mu_{s\Lambda}^2\sim \omega/s_R$ according to
\eqref{eq:initialscales}. 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
\begin{equation}
\frac{1}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds} = 
\frac{1}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds_R ds_L}|_{s_R=s_L=s},
\end{equation}
which follows from \eqref{eq:taudist} and the definition of 
the Laplace transformation.

\subsection{Double-logarithmic limit}

The double-logarithmic limit corresponds to evaluating 
the renormalization group functions at leading order 
in an expansion in $\alpha_s$:
\begin{eqnarray}
A(\nu,\mu)&=&-\frac{\alpha_s}{2\pi}\ln\frac{\mu^2}{
\nu^2}\,,
\\
S(\nu,\mu)&=&-\frac{\alpha_s}{8 \pi}\ln^2\frac{\mu^2}{
\nu^2}\,.
\end{eqnarray}
Setting scales as prescribed by \eqref{eq:initialscales}, 
the LL result \eqref{eq:mainLLresult} simplifies to 
\begin{eqnarray}
\label{eq:dbllog}
\frac{1}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds_R ds_L}&=&
\frac{\alpha_s C_F}{\pi}\frac{1}{Q s_R}
\,e^{-\frac{\alpha_s C_A}{\pi}\ln\klammer{s_L e^{\gamma_E} Q}\ln\klammer{s_R e^{\gamma_E} Q}}\int_{\sigma}^Q\frac{ d\omega}{\omega}\,
\klammer{\frac{\omega}{Q}}^{\!\frac{\alpha_s}{\pi}\klammer{C_F-C_A}\ln\klammer{s_R e^{\gamma_E} Q}
}
\nonumber\\
&=&
\frac{C_F}{C_F-C_A}\frac{1}{Q s_R  \ln(s_R e^{\gamma_E} Q)}
\exp\,\bigg[-\frac{\alpha_s C_A}{\pi}\ln\klammer{s_L e^{\gamma_E} Q}\ln\klammer{s_R e^{\gamma_E} Q}\bigg]
\nonumber\\
&&\times\bigg\{
1-\klammer{\frac{\sigma}{Q}}^{\!\frac{\alpha_s}{\pi}\klammer{C_F-C_A}\ln\klammer{s_R e^{\gamma_E} Q}} \bigg\}\,.
\end{eqnarray}
The scale of $\alpha_s$ remains undetermined in the 
double-logarithmic approximation.

To obtain the thrust distribution we set $s_R=s_L=s$ and 
choose $\sigma=1/(s e^{\gamma_E})$. 
For the inverse Laplace-transformation from $s$ to $\tau=1-T$ we use that the Laplace transform of $\tau^a$ is $\Gamma(1+a)/(sQ)^{1+a}$. Taylor-expanding this relation in $a$ gives 
\begin{equation}
  \frac{1}{sQ}\ln^n\left(\frac{1}{s e^{\gamma_E} Q}\right) \rightarrow \ln^n\tau-\frac{\pi^2}{12}n(n-1)\ln^{n-2}\tau+{\cal O}(\ln^{n-3}\tau)\,,
\end{equation}
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 $s e^{\gamma_E}\to 1/(Q\tau)$ resulting in 
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{ d\sigma}{ d\tau}|_{\rm DL} &=&
\frac{C_F}{C_F-C_A}\frac{1}{\ln(1/\tau)}
\,e^{-\frac{\alpha_s C_A}{\pi}\ln^2\tau}\,
\left\{
1-e^{-\frac{\alpha_s}{\pi}\klammer{C_F-C_A}\ln^2\tau} 
\right\}\,.
\label{eq:doublelogresult}
\end{eqnarray}
This agrees with the previous 
result \cite{Moult:2019uhz,Beneke:2020ibj}, which, 
however, was obtained through $d$-dimensional factorization 
and consistency relations.\footnote{The overall sign of (E.13) 
in \cite{Beneke:2020ibj} is given incorrectly. The signs of all 
previous equations leading to (E.13) are, however, correct.}

\subsection{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
\begin{align}
&A\klammer{\nu,\mu}=\frac{2}{\beta_0}\ln\frac{\alpha_s\klammer{\mu}}{\alpha_s\klammer{\nu}},\\
\nonumber&S\klammer{\nu,\mu}=\frac{4\pi}{\beta_0^2\alpha_s\klammer{\nu}}\left[1-\frac{\alpha_s\klammer{\nu}}{\alpha_s\klammer{\mu}}-\ln\frac{\alpha_s\klammer{\mu}}{\alpha_s\klammer{\nu}}\right],
\end{align}
and the strong coupling $\alpha_s$ at a scale $\nu$ is 
given by
\begin{align}
\alpha_s\klammer{\nu}=\frac{\alpha_s(\mu)}{1+\frac{\beta_0}{2\pi}\alpha_s\klammer{\mu}\ln\frac{\nu}{\mu}}
\label{eq:onelooprunningalphas}
\end{align}
in terms of the coupling at a reference scale $\mu$. 
In the following $\alpha_s$ without argument 
refers to $\alpha_s(\mu)$.

These expressions are now to be used in 
\eqref{eq:mainLLresult}, which upon 
implementing \eqref{eq:initialscales} (with ``$\sim$'' replaced 
by ``='' and $s_i$ by $s_i e^{\gamma_E}$)
reads\footnote{In an abuse of notation, we here 
give the scale arguments of $\alpha_s$, $A$ and $S$ as 
$\mu^2$ rather than $\mu$.}
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds_R ds_L}|_{\rm LL}&=&
\frac{\alpha_s(Q/(s_L e^{\gamma_E})) C_F}{\pi}\,\frac{1}{Q s_R}
\nonumber\\
&&\hspace*{-2cm}
\times\,
\exp\left[4C_F S(Q^2,\frac{Q}{s_R e^{\gamma_E}})
+4C_A S(\frac{1}{s_L s_R e^{2\gamma_E}},\frac{Q}{s_L e^{\gamma_E}})\right]
\nonumber\\
&&\hspace*{-2cm}
\times\int_\sigma^Q\frac{d\omega}{\omega}\,
\exp\left[-4\klammer{C_F-C_A}S(\omega Q,\frac{\omega}{s_R e^{\gamma_E}})\right]
\nonumber\\
&&\hspace*{-1.2cm}\times\,
\klammer{s_R e^{\gamma_E} Q}^{
2C_F A\klammer{\omega/s_R  e^{\gamma_E},Q/s_R  e^{\gamma_E}}
+2C_A A\klammer{Q/s_L  e^{\gamma_E},\omega/s_R  e^{\gamma_E}}}
\,.
\label{eq:LLresult}
\end{eqnarray}
This is our main result for the resummed two-hemisphere 
invariant mass distribution in the gluon-thrust region. 
The integral over $\omega$ can be performed only numerically. 

An approximate analytic expression can be obtained by expanding 
the one-loop running coupling to linear order in $\beta_0$ in 
the exponent and prefactor, resulting in 
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds_R ds_L}|_{\rm LL}
&=& 
\frac{\alpha_s C_F}{\pi}\,\frac{1}{Q s_R}\,
\,\bigg(1-\frac{\beta_0\alpha_s}{4\pi}
\ln\frac{Q}{s_L e^{\gamma_E} \mu^2}\bigg)
\nonumber\\
&&\hspace*{-1.5cm}\times\,
\exp\Bigg[-\frac{\alpha_s C_A}{\pi}\ln\klammer{s_L e^{\gamma_E} Q}\ln\klammer{s_R e^{\gamma_E} Q}
\bigg(1-\frac{\beta_0\alpha_s}{4\pi}
\ln\frac{Q}{\sqrt{s_R s_L} e^{\gamma_E} \mu^2}\bigg)\,\Bigg]
\nonumber\\
&&\hspace*{-1.5cm}\times\,\int_{\sigma}^Q\frac{ d\omega}{\omega}
\klammer{\frac{\omega}{Q}}^{\!\frac{\alpha_s}{\pi}\klammer{C_F-C_A}\ln\klammer{s_R e^{\gamma_E} Q}
\,\big(1-\frac{\beta_0\alpha_s}{4\pi}
\ln\frac{Q\sqrt{\omega}}{\sqrt{s_R e^{\gamma_E}}\mu^2}\big)}
\nonumber\\[0.1cm]
&&\hspace*{-2cm}=\,\frac{C_F}{Q s_R}\exp\Bigg[-\frac{\alpha_s}{\pi}C_A\ln\klammer{s_L e^{\gamma_E} Q}\ln\klammer{s_R e^{\gamma_E} Q}\klammer{1-\frac{\beta_0\alpha_s}{4\pi}\ln\frac{Q}{\sqrt{s_L s_R} e^{\gamma_E}\mu^2}}\Bigg]
\nonumber\\
&&\hspace*{-1.5cm}\times\,
\frac{\sqrt{2\pi}\left(1-\frac{\beta_0\alpha_s}{4\pi}\ln\frac{Q}{s_L e^{\gamma_E} \mu^2}\right)}
{\sqrt{\beta_0}\sqrt{C_A-C_F}\sqrt{\ln(s_R e^{\gamma_E} Q})}
\klammer{s_R e^{\gamma_E}  Q}^{\frac{2\klammer{C_F-C_A}}{\beta_0}\klammer{1-\frac{\beta_0\alpha_s}{4\pi}\ln\frac{Q^{\frac{3}{2}}}{\sqrt{s_R e^{\gamma_E}}\mu^2}}^2}\nonumber\\[0.2cm]
&&\hspace*{-1.5cm}\times\,\Bigg[\,
\mathrm{erfi}\klammer{
\sqrt{2 \,\frac{C_A-C_F}{\beta_0} \ln(s_R e^{\gamma_E} Q)}
\,\bigg(1-\frac{\beta_0\alpha_s}{4\pi}\ln\frac{\sigma\sqrt{Q}}{\sqrt{s_R e^{\gamma_E}}\mu^2}\bigg)}\nonumber\\
&&\hspace*{-0.8cm}-\,\mathrm{erfi}\klammer{
\sqrt{2 \,\frac{C_A-C_F}{\beta_0} \ln(s_R e^{\gamma_E} Q)}
\,\bigg(1-\frac{\beta_0\alpha_s}{4\pi}\ln\frac{Q^{\frac{3}{2}}}{\sqrt{s_R e^{\gamma_E}}\mu^2}\bigg)}\Bigg]\,,
\label{eq:linearbeta0LLresults}
\end{eqnarray}
where 
\begin{equation}
\mathrm{erfi}(x) = -\frac{2 i}{\sqrt{\pi}}\int_0^{i x} dt\,e^{-t^2}\,.
\end{equation}
Once again, the thrust distribution in Laplace-space is obtained 
by identifying $s_R=s_L\equiv s$, in which case we also 
choose $\sigma=1/(s e^{\gamma_E}) $. 

\begin{figure}[p]
\begin{center}
\includegraphics[width=0.7\textwidth]{figs/taudistribution}\\[0.7cm]
\includegraphics[width=0.7\textwidth]{figs/laplacesdistribution}\\[0.2cm]
\caption{\label{fig:numerics}
Upper panel: Gluon thrust distribution 
$\frac{1}{\sigma_0}\frac{ d\sigma}{ d\tau}$ in the double-logarithmic 
approximation including scale variation (light-blue band from 
\eqref{eq:doublelogresult} and light-green band from the numerical 
inverse Laplace transform of \eqref{eq:dbllog}), and the 
LL  approximation (red line). 
The tree-level result is shown in dotted-blue.
Lower panel: Laplace-space gluon thrust distribution  (times $Q s$),
$\frac{Q s}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds}$, 
in the LL approximation (red line), allowing for a variation of the 
initial scales as described in the text (light-red band). Also shown is the 
truncation of the LL result to linear order in $\beta_0$ from 
\eqref{eq:linearbeta0LLresults} (blue), but the difference is hardly 
visible.}
\end{center}
\end{figure}

To display the effect of resummation and of running coupling 
effects, we show in the upper panel of 
Fig.~\ref{fig:numerics} the gluon-thrust 
distribution $\frac{1}{\sigma_0}\frac{ d\sigma}{ d\tau}$ in 
the leading $\mathcal{O}(\alpha_s)$ approximation 
(blue-dotted line) and in the resummed double-logarithmic 
approximation \eqref{eq:doublelogresult} (light-blue band). 
We adopt $Q=m_Z=91.1876\,$GeV corresponding to final states from 
hadronic $Z$-boson decay, and $\alpha_s(m_Z)=0.1179$. The coupling 
is evolved with one-loop accuracy and $n_l=5$ according to 
\eqref{eq:onelooprunningalphas}. In the leading-order approximation 
we choose $\alpha_s$ at the collinear scale $\sqrt{\tau} Q$. 
Since the scale of $\alpha_s$ is not fixed in the double-logarithmic 
approximation, we vary it from the soft scale $\tau 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 \eqref{eq:dbllog} in Laplace space. 
To this end, we compute 
\begin{equation}
\frac{1}{\sigma_0}\frac{ d\sigma}{ d\tau} = \frac{Q}{2\pi i}
\int_{c-i\infty}^{c+i\infty} \!\! ds\, e^{s\tau Q}\,
\frac{1}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds_R ds_L}|_{s_R=s_L=s},
\label{eq:inverseLaplacedist}
\end{equation}
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 \cite{Catani:1996yz}.
We also set 
$\sigma = 1/(s e^{\gamma_E})$ in \eqref{eq:dbllog}. We observe 
a significant difference between the two bands, which must 
be attributed to formally next-to-leading logarithmic effects 
in the relation between Laplace and momentum space. 

Finally, the 
upper plot shows (red line) the 
full leading-logarithmic result \eqref{eq:LLresult} 
(with the $\omega$ integral 
evaluated numerically for $\sigma = 1/(s e^{\gamma_E})$) after 
taking the inverse Laplace transformation 
\eqref{eq:inverseLaplacedist} 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 $\tau\approx 0.03$, 
since the DL result at the soft scale 
becomes sensitive to the strong coupling regime 
at smaller values $\tau$. 
We also recall that perturbative 
resummation requires 
$Q\tau\gg \Lambda_{\rm QCD}\approx 0.5\,$GeV and $\tau\ll 1$.

We can gain some insight on the importance of intrinsic 
{\em next}-to-leading logarithms by varying in the LL result the 
various matching scales around the 
values adopted in \eqref{eq:LLresult}. For this purpose, we 
vary the three pairs of scales $(\mu_h, \mu_{h\Lambda})$, 
$(\mu_c,\mu_{\bar c})$, $(\mu_s, \mu_{s\Lambda})$ 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,\frac{1}{2},\frac{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 
$\frac{Q s}{\sigma_0}\frac{\widetilde{ d\sigma}}{ ds}$ in the 
lower panel of Fig.~\ref{fig:numerics} as the light-red 
band around the red curve (LL) that represents 
\eqref{eq:LLresult}. For comparison the tree-level (LO) and 
linear-$\beta_0$ truncation \eqref{eq:linearbeta0LLresults} 
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 $\mathcal{O}(\alpha_s)$ anomalous 
dimension of 
$\mathcal{J}_c(p^2,\omega,\omega')$ and $S_{\rm NLP}(l_+, l_-,
\omega,\omega')$  
for general $\omega^{(\prime)}$, which require 
two-loop calculations beyond the scope of this work, as 
well as the  $\mathcal{O}(\alpha_s^2)$ cusp parts of the anomalous 
dimensions of the NLP objects.


\section{Conclusion}
\label{sec:conclusion}

The lack of convergence of the convolution integrals appearing in subleading-power factorization theorems is currently the biggest obstacle to resummation of the power-suppressed large logarithmic corrections in collider physics. This applies in particular to the classic 
$1\to 2$ and $2\to 1$ processes, such as event shapes in 
$e^+ e^-$ annihilation,  deep-inelastic scattering for $x\to 1$, 
and the Drell-Yan threshold. 

In earlier work \cite{Beneke:2020ibj} on off-diagonal deep-inelastic scattering for $x\to 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
\cite{Moult:2019uhz,Beneke:2020ibj}, 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 sub\-leading-power soft and jet 
functions defined at the cross-section level. The above-mentioned 
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$_{\rm I}$ factorization theorems and hence opens the path 
to systematic next-to-leading-power 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  \cite{Liu:2019oav,Liu:2020wbn}, 
which however is a SCET$_{\rm II}$ process, where factorization 
applies to the amplitude. For the $1\to 2$ and $2\to 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 \cite{Beneke:2018gvs,Beneke:2019oqx}, yet their 
divergence as $\epsilon\to 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$_{\rm 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 \cite{Moult:2019uhz,Beneke:2020ibj} and evaluated the numerical impact of the new 
LL corrections. Our main technical result for gluon thrust is 
\eqref{eq:LLresult}, 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$_{\rm 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 $\mathcal{O}(\alpha_s)$, i.e., at the two-loop level 
(counting phase-space integrals as loops). 
For the {\em diagonal} processes, the leading NLP logarithms 
have been considered for the thrust distribution and the 
DY threshold in \cite{Moult:2018jjd,Beneke:2018gvs,Bahjat-Abbas:2019fqa,Beneke:2019mua}, 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 \cite{Beneke:2019oqx}. 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. 

\subsubsection*{Acknowledgements} 

This work has been supported in part by the Excellence Cluster 
ORIGINS funded by the Deutsche Forschungsgemeinschaft 
(DFG, German Research Foundation) under Ger\-many's Excellence Strategy --EXC-2094 --390783311. S.J. is supported by the UK Science and Technology Facilities Council (STFC) under grant ST/T001011/1. R.S. is supported by the United States Department of Energy under Grant Contract DE-SC0012704. L.V. is supported by Fellini Fellowship for Innovation at INFN, funded by the European Union's Horizon 2020 research programme under the Marie Sk\l{}odowska-Curie Cofund Action, grant agreement no. 754496. J.W. was supported in part by the National Natural Science Foundation of China (No.~12005117) 
and the Taishan Scholar Foundation of Shandong province (tsqn201909011).
Figures were drawn with \texttt{Jaxodraw}~\cite{Binosi:2008ig}. 


\begin{appendix}

\section{\boldmath $C^{\rm B1}$ and $D^{\rm B1}$}
\label{sec:DB1}

In this appendix we provide explicit results for the
hard matching coefficient $C^{\rm B1}_{i}(Q^2,r)$ associated to
the SCET operators defined in the second line of 
\eqref{eq:hardmatching} and \eqref{eq:BtypeDirac}. In 
momentum space, 
\begin{equation}
J^{B1,\mu}_i(r) = \bar\chi_{\bar c}( r)\Gamma_i^{\mu\nu}\mathcal{A}_{\nu\perp c}\chi_{\bar c}(\bar r),\quad i=1,2\,,
\label{eq:B1momspace}
\end{equation}
where $r$ and $ \bar r=1-r$ denote the momentum fractions of the anti-collinear quark and anti-quark, respectively, and
\begin{equation}
  \Gamma_1^{\mu\nu} = \frac{\slashed{n}_-}{2}\gamma_\perp^\nu\gamma_\perp^\mu,\qquad
  \Gamma_2^{\mu\nu} = \frac{\slashed{n}_-}{2}\gamma_\perp^\mu\gamma_\perp^\nu\,.  
\end{equation}
We then obtain the one-loop correction and one-loop evolution equation 
of the refactorization coefficient $D^{\rm B1}(p^2)$ by taking 
the $r\to 0$ limit of the corresponding results for  
$C^{\rm B1}_{1}(Q^2,r)$ and \eqref{eq:B1fact1},
\begin{eqnarray}
\label{eq:DB1def2}
C^{\rm B1}_1(Q^2,r) &=& C^{\rm A0}(Q^2)\,\frac{D^{\rm B1}(rQ^2)}{r}
+\text{endpoint-regular}\,.
\end{eqnarray}
Here endpoint-regular contributions diverge at most logarithmically 
in the limit $r\to 0$.


\subsection{\boldmath One-loop results}
\label{sec:CB1}

Allowing for the one-loop correction to \eqref{eq:CB11tree}, 
we write 
\begin{equation}\label{CB1loopgeneralform}
C_1^{\rm B1}(Q^2,r) =
\frac{1}{r} \left(1+\frac{\alpha_s}{4\pi} \,\Delta^{\rm B1}(Q^2,r)
+\mathcal{O}(\alpha_s^2)\right).
\end{equation}
We recall that charge conjugation implies 
$C_2^{\rm B1}(Q^2,r) = -C_1^{\rm B1}(Q^2,\bar{r})$, 
hence it suffices to 
discuss the case $i=1$. A standard matching calculation of 
the $\gamma^*\to q\bar{q}g$ amplitude results in 
\begin{eqnarray}\label{deltaB1full} \nonumber
\Delta^{\rm B1}(Q^2,r) &=&
\bigg[\bigg(C_F-\frac{C_A}{2}\bigg)
\bigg(- \frac{2}{\epsilon^2} + \frac{2}{\bar r \epsilon}
+ \frac{4 r}{\bar r} +\frac{\pi^2}{6} \bigg)
+C_F \frac{2}{\bar r} \bigg]
\bigg(-\frac{\mu^2}{Q^2}\bigg)^{\epsilon} \\ \nonumber
&&\hspace{-2.0cm}+\,\bigg[C_F \bigg(\frac{2}{\epsilon^2} - \frac{2 r}{\bar r \epsilon}
- \frac{1 + 5 r}{\bar r} -\frac{\pi^2}{6} \bigg)
+C_A \bigg(- \frac{2}{\epsilon^2} + \frac{r}{\bar r \epsilon}
+ \frac{1 + r}{\bar r} + \frac{\pi^2}{6}\bigg)\bigg]
\bigg(-\frac{\mu^2}{r Q^2}\bigg)^{\epsilon} \\
&&\hspace{-2.0cm}+\,\bigg[C_F \bigg(- \frac{5}{\epsilon} -10 \bigg)
+C_A \bigg(- \frac{1}{\epsilon^2} + \frac{1}{\epsilon} + \frac{\pi^2}{12}\bigg)\bigg]
\bigg(-\frac{\mu^2}{\bar r Q^2}\bigg)^{\epsilon}
+\mathcal{O}(\epsilon)\,.
\end{eqnarray}
The limit $r \to 0$ is of particular interest, and results in
\begin{equation}
\Delta^{\rm B1}(Q^2,r) \stackrel{r\to 0}{=} C_F \, \Delta^{A0}
\bigg(-\frac{\mu^2}{Q^2}\bigg)^{\epsilon}
+\big(C_F- C_A \big)
\bigg(\frac{2}{\epsilon^2} -1 - \frac{\pi^2}{6} \bigg)
\bigg(-\frac{\mu^2}{r Q^2}\bigg)^{\!\epsilon}
+\mathcal{O}(\epsilon,r),
\end{equation}
which shows that the B1 coefficient depends on two scales, 
$Q^2$ and $rQ^2$. We can identify 
\begin{equation}
\Delta^{A0} = - \frac{2}{\epsilon^2} - \frac{3}{\epsilon}
-8 +\frac{\pi^2}{6} + \mathcal{O}(\epsilon)
\end{equation}
with the one-loop contribution to the leading-power 
hard matching coefficient coefficient $C^{\rm A0}$:
\begin{equation}
C^{\rm A0}(Q^2) = 1 + \frac{\alpha_s C_F }{4\pi} \Delta^{A0}
\left(-\frac{\mu^2}{Q^2}\right)^{\!\epsilon} +\mathcal{O}(\alpha_s^2).
\end{equation}
Comparing to \eqref{eq:B1fact1}, we find 
\begin{equation}
D^{\rm B1}(p^2) = 1 + \frac{\alpha_s}{4\pi}\,
\big(C_F- C_A \big)
\bigg(\frac{2}{\epsilon^2} -1 - \frac{\pi^2}{6} \bigg)
\left(-\frac{\mu^2}{p^2}\right)^{\!\epsilon} +\mathcal{O}(\alpha_s^2).
\label{eq:DB11loopapp}
\end{equation}
For completeness, we check that the $r\to 1$ limit
of $\Delta^{\rm B1}(Q^2,r)$ does not contain $1/\bar{r}$ terms:
\begin{eqnarray} \nonumber
\Delta^{\rm B1}(Q^2,r) &\stackrel{r\to 1}{=}&
\bigg[C_F \bigg(\frac{2}{\epsilon} -1 \bigg)
+C_A \bigg(-\frac{1}{\epsilon^2} - \frac{1}{\epsilon}
+2 + \frac{\pi^2}{12} \bigg) \bigg]
\bigg(-\frac{\mu^2}{Q^2}\bigg)^{\epsilon} \\
&&\hspace{-2.0cm}+\,\bigg[C_F \bigg(- \frac{5}{\epsilon} -10 \bigg)
+C_A \bigg(- \frac{1}{\epsilon^2} + \frac{1}{\epsilon} + \frac{\pi^2}{12}\bigg)\bigg]
\bigg(-\frac{\mu^2}{\bar r Q^2}\bigg)^{\epsilon}
+\mathcal{O}(\epsilon,\bar r)\,.
\label{eq:DeltaB1xto1limit}
\end{eqnarray}


\subsection{\boldmath Derivation of the evolution equation of the 
asymptotic refactorization coefficient $D^{\rm B1}$}
\label{app:asympDB1RGE}


Our goal is to show that the RGE \eqref {eq:DBQrge} for the 
matching coefficient $D^{\rm B1}(r Q^2)$ can be derived 
from the standard evolution equation of the hard-matching 
coefficient $C_i^{\rm B1}(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~\cite{Beneke:2019kgv}.

The renormalization of SCET operators of the B1 type has 
been studied in detail 
in~\cite{Beneke:2017ztn,Beneke:2018rbh}, except for 
the fermion-number zero case, to which the operators 
$J^{B1,\mu}_i(r)$ in \eqref{eq:B1momspace} belong.
The anomalous dimension matrix is determined by computing the $\overline{\text{MS}}$ renormalization factors 
$Z_{ij}^{\mu\nu}(r,s)=\delta_{ij}g^{\mu\nu}_\perp\delta(r-s)+\delta Z_{ij}^{\mu\nu}(r,s)$ defined via
\begin{equation}
  J^{B1,\mu}_{i,\text{ren}}(r) = \int ds\, Z_{ij}^{\mu\nu}(r,s) J^{B1,\nu}_{j,\text{bare}}(s)\,.
\end{equation}
Following the lines of~\cite{Beneke:2017ztn,Beneke:2018rbh}, we find at one-loop order $\delta Z_{ij}^{\mu\nu}(r,s)=g^{\mu\nu}_\perp \delta Z_{ij}(r,s)$ with 
\begin{eqnarray}
  \delta Z_{11}(r,s) &=& \delta Z_{22}(r,s) = -\frac{\alpha_s}{2\pi}\delta(r-s)\left[\frac{C_A}{\epsilon^2}+\frac{C_A}{\epsilon}\ln\left(\frac{\mu^2}{-Q^2-i\varepsilon}\right)-\frac{C_A}{2\epsilon}\ln(r\bar r)+\frac{3C_F}{2\epsilon}\right]\nonumber\\
  && {} +\frac{\gamma(r,s)}{\epsilon}\,,\\
  \delta Z_{12}(r,s) &=& \delta Z_{21}(r,s) = \frac{\alpha_s}{2\pi}\Delta_F\frac{2T_F}{\epsilon}r\bar r \,,
\end{eqnarray}
where
\begin{eqnarray}\label{eq:gamma}
  \gamma(r,s) &=& \frac{\alpha_s}{2\pi}\,\Bigg\{ \left(\frac{C_A}{2}-C_F\right)\Bigg[ \frac{\bar r}{\bar s}\left(\left(\frac{\theta(r-s)}{r-s}\right)_+ + \theta(r-s)\right)\nonumber\\
  && {} + \frac{r}{ s}\left(\left(\frac{\theta(s-r)}{s-r}\right)_+ + \theta(s-r)\right) \Bigg]
   +2\Delta_F T_F r\bar r\Bigg\} \,.
\end{eqnarray}
The momentum fractions have support in the interval $0\leq r,s\leq 1$.
Furthermore, $T_F=1/2$, and $\Delta_F=1\;(0)$ in the flavour-singlet (non-singlet) projection of the anti-collinear 
$\bar\chi_{\bar c}\chi_{\bar c}$ fields in $J^{B1,\mu}_i$. The expression 
in square brackets defines a symmetric distribution, which applies 
to integration over test functions of $r$ or $s$, 
with definition 
\begin{eqnarray}
\int dr\, f(r) \,[D(r-s)]_+ &=& 
\int dr\, D(r-s) \,(f(r)-f(s))\,,
\\
\int ds\, [D(r-s)]_+ \,g(s) &=& 
\int ds\, D(r-s) \,(g(s)-g(r))\,.
\end{eqnarray} 
The definition applies to any integration range, which is left unspecified here. In~\eqref{eq:gamma} the integrations are limited to
$0\leq r,s\leq 1$, but we will also consider the case $0\leq r,s\leq \infty$ below.

The anomalous dimension can then be obtained as
\begin{equation}
  {\bf \Gamma} = {\bf Z}\frac{d}{d\ln\mu}{\bf Z}^{-1}\,,
\end{equation}
giving the RGEs
\begin{equation}
  \frac{d}{d\ln\mu} C^{\rm B1}_j(Q^2,s) = \int_0^1 dr \, C^{\rm B1}_i(Q^2,r) \Gamma_{ij}(r,s) \,.
\end{equation}
Here we already used that the current renormalization for the specific operators $J^{\rm B1,\mu}_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^\text{B1}_i(Q^2,r)$ with $\delta Z_{ij}(r,s)$ agrees with the
divergent part of the explicit one-loop result for $C^\text{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\to 0$ ($\bar r\to 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^\text{B1}$. 

The all-order identity $C^{\rm B1}_2(Q^2,r) = 
- C^{\rm B1}_1(Q^2,\bar r)$ is consistent with the relations 
$\delta Z_{11}(r,s)=\delta Z_{22}(r,s)=\delta Z_{22}(\bar r,\bar s)$ and $\delta Z_{12}(r,s)=\delta Z_{21}(r,s)=\delta Z_{21}(\bar r,\bar s)$. Furthermore it can be checked that $\delta Z_{12}(r,s)$ does not contribute to the renormalization of $D^{\rm B1}(r Q^2)$, due to the independence on $s$ (see below).
Therefore it is sufficient to consider  $Z(r,s)\equiv \delta(r-s)+\delta Z_{11}(r,s)$ in the following.

The matching coefficient $C^{\rm B1}_1(Q^2,r)$ for the problem at hand features a power-like endpoint singularity $\propto 1/r$ for $r\to 0$, that may be
modulated by logarithmic corrections, and is captured by the $D^{\rm B1}$ contribution in~\eqref{eq:B1fact1}.
In order to extract the renormalization factor for $D^{\rm B1}$, we are interested in the convolution with $\Gamma_{11}(r,s)$ (or equivalently $Z(r,s)$ at order $\alpha_s$).
The result of the convolution can be split into an endpoint divergent piece $\propto 1/s$, and a regular part. The former describes the renormalization of $D^{\rm B1}$, and the latter yields a mixing of $D^{\rm B1}$ into the endpoint-finite part of $C^{\rm 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
\begin{equation}\label{eq:fzF}
  \int dr \,\frac{f(r)}{r} Z(r,s) = \frac{F(s)}{s}+\text{endpoint-regular}\,.
\end{equation}
Here $f(r)$ is an arbitrary test function defined for $r\geq 0$, that depends logarithmically on $r$, and represents the $D^{\rm 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
\begin{equation}
  F(s) = \int dr f(r) \,Z_{A0A0}\times Z_{D^{\rm B1}D^{\rm B1}}(r,s)\,,
\end{equation}
with a renormalization factor $Z_{D^{\rm B1}D^{\rm B1}}(r,s)$, and 
$Z_{A0A0}$ the renormalization constant for the leading-power 
coefficient $C^{\rm A0}(Q^2)$, being given at order $\alpha_s$ by
\begin{equation}
  Z_{A0A0}=1-\frac{\alpha_s}{2\pi}\left[\frac{C_F}{\epsilon^2}+\frac{C_F}{\epsilon}\ln\left(\frac{\mu^2}{-Q^2-i\varepsilon}\right)+\frac{3C_F}{2\epsilon}\right]\,.
\end{equation}
The strategy for obtaining $Z_{D^{\rm B1}D^{\rm B1}}$ is to apply the method of regions \cite{Beneke:1997zp} to~\eqref{eq:fzF}, and expanding the integrand such that a homogeneous scaling in the limit $r,s\to 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 $\gamma(r,s)$,
\begin{eqnarray}\label{eq:fg}
  \int dr \,\frac{f(r)}{r} \gamma(r,s) &=& \frac{\alpha_s}{2\pi}\left(\frac{C_A}{2}-C_F\right)  \Bigg[ \int_s^1 dr  \left(\frac{\frac{\bar r}{\bar s} \frac{f(r)}{r}-\frac{f(s)}{s}}{r-s} + \frac{\bar r}{\bar s} \frac{f(r)}{r}\right)\nonumber\\
  && \hspace*{-2cm} + \int_0^s dr  \left(\frac{\frac{r}{ s}\frac{f(r)}{r}-\frac{f(s)}{s}}{s-r} + \frac{r}{s}\frac{f(r)}{r}\right) \Bigg]
  + \frac{\alpha_s}{\pi}\Delta_F T_F \int_0^1 dr\, r\bar r \,\frac{f(r)}{r}\,.\quad
\end{eqnarray}
The integrals in the second line exist due to the explicit factor $r$ in the second line of~\eqref{eq:gamma}, 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\to 0$.
The first integral in the second line of~\eqref{eq:fg} 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 $\int_0^s dr f(r)/s$. 
In addition, the integral in the ``singlet'' term proportional to $\Delta_F$ is independent of $s$ and can therefore not contribute to $F(s)$ in~\eqref{eq:fzF}.

Let us now turn to the first line of~\eqref{eq:fg}. Expanding in the region of small $r,s$ allows us to replace $\bar r/\bar s\to 1$. Next, we see that the integral $\int_s^1 dr f(r)/r$ contributes only to the endpoint-finite part
on the right-hand side of~\eqref{eq:fzF}, and can therefore be dropped in the present discussion. However, the first term in the round bracket in the first line of~\eqref{eq:fg} yields contributions that
scale as $1/s$. To isolate them, we replace the upper integration boundary $1$ by a Heaviside function $\theta(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 $+\infty$. At this point, we use that the test function $f(r)$ depends logarithmically on $r$. It can therefore be extended over the interval $0\leq r <\infty$, in accordance with the method-of-region expansion.
However, this procedure is too naive, since it implies that the integral diverges for $r\to\infty$. This problem is linked to the observation that the respective integral yields contributions that are logarithmically enhanced as $\ln(s)/s$
for $s\to 0$. They can in turn be isolated by the rewriting
\begin{equation}
  \int_s^1 dr  \,\frac{f(r)/r-f(s)/s}{r-s} = \int_s^1 dr  \, 
\frac{f(r)-f(s)}{r(r-s)} + \frac{f(s)}{s}\ln(s)\,.
\end{equation}
After this rewriting we can extend the integration range of the remaining integral to $s\leq r<\infty$, yielding the desired homogeneous scaling.\footnote{We note the similarity of this 
procedure to the extraction of the ``asymptotic kernel'' for 
the QED-generalized light-cone distribution amplitude of 
a light meson in \cite{Beneke:2021pkl}, which involves 
a colour-singlet but electrically charged operator.}

Note that the relation above could equivalently be written at the level of plus-distributions, 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
\begin{eqnarray}
  \int dr \,\frac{f(r)}{r} \gamma(r,s)\Big|_{\propto \ln()/s} &=& \frac{\alpha_s}{2\pi}\left(\frac{C_A}{2}-C_F\right)  \Bigg[ \int_s^\infty dr  \left(\frac{f(r)-f(s)}{r(r-s)} \right)+ \frac{f(s)}{s}\ln(s)\nonumber\\
  && {} +  \int_0^s dr  \left(\frac{f(r)-f(s)}{s(s-r)} \right) \Bigg]\,.
\end{eqnarray}
By comparing with~\eqref{eq:fzF} we can read off the renormalization coefficient of $D^{\rm B1}$,
\begin{eqnarray}
  Z_{D^{\rm B1}D^{\rm B1}}(r,s) &=& \frac{\alpha_s(C_F-C_A)}{2\pi}\,\delta(r-s)\left[\frac{1}{\epsilon^2}+\frac{1}{\epsilon}\ln\left(\frac{\mu^2}{-rQ^2-i\varepsilon}\right)\right]+\frac{\hat\gamma(r,s)}{\epsilon}\,,\quad
\end{eqnarray}
where
\begin{eqnarray}
  \hat\gamma(r,s) &=& \frac{\alpha_s}{2\pi} \left(\frac{C_A}{2}-C_F\right) 
s \left[\frac{\theta(r-s)}{r(r-s)} 
  + \frac{\theta(s-r)}{s(s-r)}\right]_+\,,
\end{eqnarray}
and the plus distributions are defined with respect to the interval $0\leq r,s<\infty$ here.
The anomalous dimension is, at order $\alpha_s$, given by
\begin{equation}
  \Gamma_{D^{\rm B1}D^{\rm B1}}(r,s) = -\frac{d}{d\ln \mu}Z_{D^{\rm B1}D^{\rm B1}}(r,s) = 2\epsilon Z_{D^{\rm B1}D^{\rm B1}}(r,s)-\frac{\partial Z_{D^{\rm B1}D^{\rm B1}}}{\partial\ln\mu}\,,
\end{equation}
which gives
\begin{eqnarray}
  \Gamma_{D^{\rm B1}D^{\rm B1}}(r,s) &=&\frac{\alpha_s(C_F-C_A)}{\pi}\,\delta(r-s)\ln\left(\frac{\mu^2}{-rQ^2-i\varepsilon}\right)
   + 2\hat\gamma(r,s)\,.
\end{eqnarray}
We can express the result equivalently in terms of $\omega$ and $\hat\omega$ using $s=\omega/Q, r=\hat \omega/Q$.
This yields
\begin{equation}
  \frac{d}{d\ln\mu}D^{\rm B1}(\omega) = \int_0^\infty \!d\hat\omega 
\,\gamma_D(\hat\omega,\omega)D^{\rm B1}(\hat\omega)\,,
\label{eq:DB1RGEinomega}
\end{equation}
with
\begin{eqnarray}
  \gamma_D(\hat\omega,\omega) &=& \frac{\alpha_s(C_F-C_A)}{\pi}\delta(\omega-\hat\omega)\ln\left(\frac{\mu^2}{-\omega Q-i\varepsilon}\right)\nonumber\\
  && +\,\frac{\alpha_s}{\pi} \left(\frac{C_A}{2}-C_F\right)
 \omega \left[\frac{\theta(\hat\omega -\omega)}{\hat\omega (\hat\omega -\omega)}
  + \frac{\theta(\omega-\hat\omega )}{\omega(\omega-\hat\omega )}\right]_+\,.
\label{eq:DB1RGEinomegagamma}
\end{eqnarray}
This result is consistent with Eq.~(3.2) in~\cite{Liu:2021mac} obtained by a different method, when rewriting $\hat\omega\to x \omega$ and $\omega Q\to 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 $\gamma_D$ in \eqref{eq:DB1RGEinomega}.



\section{Alternative version of the endpoint-finite 
factorization formula}
\label{version1endpintfact}

In this appendix consider the alternative split of the integration regions mentioned in Sec.~\ref{sec:endpointfact}, i.e. we present the factorization formula for the version (1), 
corresponding to $\omega$ {\em and} $\omega'$ smaller than 
an endpoint factorization parameter $\Lambda$ 
(integral $I_1$), and the complement region $I_2$, see Fig.~\ref{fig:overlap}.

Subtracting the complement region $I_2$, the A-type 
term takes the from 
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{\rm A-type} &=& 
\frac{2C_F}{Q} \,f(\epsilon)\,|C^{\rm A0}(Q^2)|^2 \,
\widetilde{\mathcal{J}}_{\bar c}^{(\bar q)}(s_R)
\,\int d\omega d\omega'\,
\nonumber\\
&&\hspace*{-3cm}\times\,
\bigg\{\,\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
\,\widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')
\nonumber\\
&&\hspace*{-2cm}- 
\left[1-\theta(\Lambda-\omega)\theta(\Lambda-\omega')\right]
\llbracket \widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega') \rrbracket
\left \llbracket S_{\rm NLP}(s_R,s_L,\omega,\omega')\right 
\rrbracket
\nonumber\\
&&\hspace*{-2cm}+\;\,\widetilde{\!\!\widehat{\mathcal{J}}}_{\!c}(s_{L},
\omega,\omega') 
\,\,\widetilde{\!\widehat{S}}_{\rm NLP}(s_R,s_L,\omega,\omega')\,\bigg\}\,. \quad
\label{eq:Atype_subtracted_not_expanded}
\end{eqnarray}
If we further assume that $\Lambda \gg 1/s_R, 1/s_L$, then we can keep only the leading terms in $1/(s_R\Lambda)$ and $1/(s_L\Lambda)$ and the previous equation simplifies to
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{\rm A-type} &=& 
\frac{2C_F}{Q} \,f(\epsilon)\,|C^{\rm A0}(Q^2)|^2 \,
\widetilde{\mathcal{J}}_{\bar c}^{(\bar q)}(s_R)
\,\int d\omega  d\omega'\,
\nonumber\\
&&\hspace*{-3cm}\times\,\bigg\{\,
\theta(\Lambda-\omega)\theta(\Lambda-\omega') \Big[\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
\,\widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')
\nonumber\\
&&\hspace*{-2.3cm}+\;\,\widetilde{\!\!\widehat{\mathcal{J}}}_{\!c}(s_{L},
\omega,\omega') 
\,\,\widetilde{\!\widehat{S}}_{\rm NLP}(s_R,s_L,\omega,\omega')\,\Big]
\\
&&\hspace*{-2.7cm}+\,
\big[\theta(\omega- \Lambda)\theta(\Lambda-\omega')+\theta(\omega'- \Lambda)\theta(\Lambda-\omega)\big]
\;\Big[\,\widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega')
\,\widetilde{S}_{\rm NLP}(s_R,s_L,\omega,\omega')
\nonumber\\
&&\hspace*{-2.3cm}\,- 
\llbracket \widetilde{\mathcal{J}}_{c}(s_{L},\omega,\omega') \rrbracket
\left \llbracket S_{\rm NLP}(s_R,s_L,\omega,\omega')\right 
\rrbracket +\;\,\widetilde{\!\!\widehat{\mathcal{J}}}_{\!c}(s_{L},
\omega,\omega') 
\,\,\widetilde{\!\widehat{S}}_{\rm NLP}(s_R,s_L,\omega,\omega')\Big]\,\bigg\}\,. \quad
\nonumber
\label{eq:Atype_subtracted_not_expanded2}
\end{eqnarray}

The remaining part $I_1$ of the scaleless integral \eqref{eq:scaleless_integral} must be combined withe the B-type term. To make notation more concise, we introduce $r_\Lambda \equiv \Lambda/Q$. The B-type term for $i=i'=1$ is 
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{ \begin{subarray}{c}\text{B-type}\\
 \text{ i=i'=1}\end{subarray}} 
&=& 
\frac{2C_F}{Q^2}\,f(\epsilon) \,\mathcal{J}_{c}^{(g)}(s_L)\,S^{(g)}(s_R,s_L)\, \int_0^{\infty} dr dr' \nonumber\\ &&\hspace*{-3cm}\times\,
 \bigg\{\,
\theta(1-r)\theta(1-r')C^{\rm B1*}_{1}(Q^2,r')C^{\rm B1}_{1}(Q^2,r) 
\,\mathcal{J}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)
\quad 
\nonumber \\ &&\hspace*{-2.3cm}-
\left[\theta(r_\Lambda-r)\theta(r_\Lambda-r')\right]
\llbracket C^{\rm B1*}_{1} (Q^2,r')\rrbracket_{0} \llbracket C^{\rm B1}_{1}(Q^2,r) \rrbracket_{0} \, \llbracket \mathcal{J}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)\rrbracket
\,\bigg\}\,.\qquad 
\label{eq:Btype_subtracted_11_expanded}
\end{eqnarray}
If we choose $r_\Lambda \ll 1$, then the B-type term simplifies to
\begin{eqnarray}
\frac{1}{\sigma_0}\frac{\widetilde{d\sigma}}{ds_R ds_L}|_{ \begin{subarray}{c}\text{B-type}\\
 \text{ i=i'=1}\end{subarray}} &=& 
\frac{2C_F}{Q^2} \,f(\epsilon)\,\mathcal{J}_{c}^{(g)}(s_L)\,S^{(g)}(s_R,s_L)\, \int_0^1 \!dr \int_0^1\!dr' 
\left[1-\theta(r_\Lambda-r)\theta(r_\Lambda-r')\right]
\nonumber\\ &\times&\,
 \bigg\{\,
C^{\rm B1*}_{1}(Q^2,r')C^{\rm B1}_{1}(Q^2,r) 
\,\mathcal{J}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)
\nonumber\\
&&\hspace*{0.5cm}-\,
\llbracket C^{\rm B1*}_{1} (Q^2,r')\rrbracket_{0} \llbracket C^{\rm B1}_{1}(Q^2,r) \rrbracket_{0} \, \llbracket \mathcal{J}_{\bar{c}}^{q\bar{q}(8)}(s_R,r,r^\prime)\rrbracket\,
\bigg\}
\label{eq:Btype_subtracted_11_expanded2}
\end{eqnarray}
up to corrections of 
$\mathcal{O}(\Lambda/Q)$. Just as for version (2) given in the main text, we can obtain the $i=i'=2$ case by $C_1^{\rm B1(*)}\to C_2^{\rm B1(*)} $, $\llbracket \ldots \rrbracket_0\to \llbracket \ldots \rrbracket_1$ and $r^{(')} \to \bar{r}^{(')}$.

\end{appendix}


