aainstitutetext: Department of Physics, Korea University, Seoul 02841,
Koreabbinstitutetext: Institute of Convergence Fundamental Studies and School
of Liberal Arts, Seoul National University of Science and Technology, Seoul
01811, Korea

#  $\bm{N}$-jettiness for muon jet pairs in electroweak high-energy processes

Junegone Chay<EMAIL_ADDRESS>b Taehyun Kwon<EMAIL_ADDRESS>
###### Abstract

We study the $N$-jettiness in the electroweak high-energy process for the
final muon jet pairs, $e^{-}e^{+}\rightarrow\mu^{+}\
\mathrm{jet}+\mu^{-}\mathrm{jet}$. Compared to QCD, the main difference is
that there exist additional gauge nonsinglet contributions in the weak
interaction, which make the factorization more elaborate. Especially the
nonsinglet contributions arise due to the Block-Nordsieck violation in
electroweak processes, which yields the Sudakov logarithms and the rapidity
divergence. They change the evolution of the factorized parts considerably in
the $N$-jettiness. There are two possible channels, initiated from the gauge
bosons $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, and from the electrons
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, where
$\ell$ denotes the weak doublet. The latter was discussed previously, and we
complete the analysis by studying the first. The factorization for
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ can be proceeded in a similar
way as in the factorization for
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, and
the result exhibits a rich structure. The new ingredients in this study
consist of the $W$ beam functions, and the complex color structure of the soft
functions and the hard functions. The resummation of the large logarithms is
performed by solving the renormalization group equations with respect to the
renormalization scale and the rapidity scale. In the numerical analysis, we
confine to the SU(2) weak gauge interaction, and the numerical results are
presented for both channels at next-to-leading-logarithmic accuracy including
the singlet and the nonsinglet contributions. The nonsinglet contributions
turn out to be appreciable in the 2-jettiness.

###### Keywords:

jettiness, muon jet pairs in electroweak processes, Block-Nordsieck violation,
rapidity divergence, resummation

##  1 Introduction

In high-energy collisions, energetic particles are produced as collimated
beams of particles, in the form of jets. The interplay between the constituent
partons participating in the hard scattering and the resultant hadrons forming
jets after hadronization is intertwined by the perturbative and the
nonperturbative effects of quantum chromodynamics (QCD). Since the strong
interaction is involved at every stage in the scattering process with
disparate energy scales, it is difficult to disentangle the interlocking
aspects of the strong interaction.

The proof of the factorization theorems in which the hard, collinear and soft
parts are separated in high-energy scattering has been one of the biggest
challenges in QCD. In a factorized process, the hard part describes the
contribution from the large energy of order $Q$. The collinear part takes care
of the energetic, collinear particles, which include incoming beams and final-
state jets. The soft part depicts the soft emissions interspersed between
collinear directions.

Soft-collinear effective theory (SCET) Bauer:2000ew ; Bauer:2000yr ;
Bauer:2001ct ; Bauer:2001yt  has set a new stage in proving the factorization
for various jet observables in collider physics. SCET is an effective theory,
in which collinear, soft modes are selected as the relevant degrees of
freedom. The hard degrees of freedom are integrated out to yield the hard
function. The collinear modes in different lightcone directions are decoupled,
and by redefining the collinear fields with the soft Wilson lines, the soft
sector is also decoupled at the Lagrangian level. As a result, the
factorization can be established more transparently in SCET than in full QCD.

In SCET, the phase space is divided into the collinear and the soft regions.
The invariant masses of their modes may be the same ($\mathrm{SCET_{II}}$) or
different ($\mathrm{SCET_{I}}$) depending on the physics we are interested in.
One of the important objectives in partitioning the phase space is to describe
the physics with a single scale in each kinematic region, and perform the
resummation of the large logarithms by solving the renormalization (RG)
equation.

We step up a gear to consider what happens in electroweak processes at
extremely high-energy of order 10 TeV. Still, the effect of the strong
interaction is obviously dominant at LHC. But some electroweak process can be
observed in high-energy electron-positron colliders such as CEPC
dEnterria:2016sca , ILC Djouadi:2007ik , FCC-ee Abada:2019zxq , and CLIC
Charles:2018vfv . As an example, the process
$e^{-}e^{+}\rightarrow\mu^{-}\mu^{+}$ can be described in analogy with the QCD
process $q\overline{q}\rightarrow q^{\prime}\overline{q}^{\prime}$. The
incoming “partons” in an electron and a positron possess certain energy
fractions, and they emit gauge bosons to be far off-shell before they
participate in the hard scattering. This process is described by the weak beam
functions. The energetic partons undergo a hard scattering and the final-state
particles are observed in terms of jets, described by the jet functions. And
the soft function describes the emission of soft particles between the
collinear particles.

We may suspect that it is just another copy of QCD by merely replacing the
gauge group, and nothing in particular can be gained. We claim that it is not
true. The main difference lies in the fact that hadrons appear only as color
singlets in QCD, but weak doublets such as the electrons or the neutrinos can
be observed. Because only the color singlets contribute to the observables in
QCD, the beam functions and the jet functions are obtained by taking the
matrix elements of the relevant operators between the hadronic states (color
singlets) or between the vacuum. However, both singlet and nonsinglet
contributions are involved in all these factorized components in electroweak
processes, and it makes the analysis more intriguing. We can cast this issue
on the imaginary QCD, in which there is no confinement with free quarks and
gluons. Then color nonsinglet contributions affect the jet observables and we
can ask whether the factorization still works in a consistent fashion in the
presence of these extra contributions. Of course, real QCD does not work that
way, but we can ask the same question toward the weak interaction.

There are nonsinglet contributions in the weak interaction to the beam, jet,
and soft functions as well as singlet contributions. More importantly, because
the initial or final states can have fixed charges, Sudakov logarithms can
arise from the non-cancellation of the electroweak logarithms, which is known
as the Block-Nordsieck violation in electroweak processes Ciafaloni:2000df ;
Ciafaloni:2001vt ; Ciafaloni:2006qu ; Manohar:2014vxa . Note that the Sudakov
logarithms from virtual and real contributions in QCD cancel in inclusive
processes. Due to this cancellation, it yields the Dokshitzer-Gribov-Lipatov-
Altarelli-Parisi (DGLAP) evolution of the parton distribution functions (PDF)
in QCD Lipatov:1974qm ; Gribov:1972ri ; Altarelli:1977zs ; Dokshitzer:1977sg .
The non-cancellation of the Sudakov logarithms in weak interaction affects the
renormalization behavior. For example, the nonsinglet PDF satisfies nontrivial
RG equations due to the Sudakov logarithm, while the singlet PDF still
satisfies the DGLAP equation. If we observe more exclusive, or differential
quantities with kinematic constraints on the real emissions, there is only a
partial cancellation between the virtual and real contributions, resulting in
additional Sudakov logarithms. Exclusive jet cross sections, or jet shape
observables such as $N$-jettiness are such examples. Then there appear Sudakov
logarithms even in the singlet case.

Furthermore, there is the issue of rapidity divergence Chiu:2011qc , which
arises in SCET due to the dissection of the phase space. If we add the
contributions of the factorized parts from all the phase spaces, the rapidity
divergence cancels. It is why there is no rapidity divergence in QCD because
the phase space is not divided. But aside from the fact that there is no
rapidity divergence in the singlet contributions, the rapidity divergence in
the electroweak nonsinglet contribution survives. It was first pointed out in
ref. Manohar:2018kfx  that the electroweak nonsinglet PDF has the rapidity
divergence. Because of the rapidity divergence, the double RG evolutions with
respect to the rapidity scale, as well as the conventional renormalization
scale should be solved for both the collinear functions and the soft
functions.

In this paper, we consider the $N$-jettiness in electroweak processes, in
which the Sudakov logarithms appear both in the singlet and the nonsinglet
contributions, while the rapidity logarithm appears in the nonsinglet
contributions. These result in part from the Block-Nordsieck violation, and in
part from the kinematic constraint to extract the $N$-jettiness. The
$N$-jettiness in electroweak processes at extreme high energies is analyzed
for the scattering process in which the final dijets involve a muon and an
anti-muon respectively. In order to address the main differences between QCD
and the weak interaction, we consider the $SU(2)$ weak interaction only. The
extension to the Standard Model is an interesting issue, but we will not
pursue it here. It is enough to consider SU(2) in considering the features
from the nonsinglet contributions. And we select the 2-jettiness because it
involves four lightlike directions, in which the dependence of the various
lightcone directions appears in the hard and soft functions. It offers a
nontrivial check to see the independence of the renormalization and the
rapidity scales in the 2-jettiness.

Before we proceed, we need to point out that the electron here means the
electron with a cloud of gauge bosons and all the possible leptons. It is
analogous to the proton in QCD in the sense that the proton contains gluons
and all the possible quarks (and antiquarks). The confusion may arise because
there is no distinctive terminology to distinguish the electron cloud, and the
electron itself, while a quark describes a parton in the proton in QCD. From
now on, the electron $e$ corresponds to the proton in QCD, while the electron
doublet $\ell_{e}$ describes the partons.

With this terminology in mind, we consider the 2-jettiness in the process
$e^{-}e^{+}\rightarrow\mu^{-}$-jet \+ $\mu^{+}$-jet. The final states consist
of the dijets in which there is at least a muon and antimuon in each jet. In
our previous paper Chay:2021arz , the factorization for the channel
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ has
been analyzed to show the characteristics of the nonsinglet nature in the weak
interaction. In order to analyze the 2-jettiness for the $\mu^{-}\mu^{+}$
jets, the remaining channel $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$
should be included, which is the main issue of the paper. In the factorization
for this channel, we need the additional input of the gauge-boson beam
functions, and the hard and soft functions associated with this channel.
Especially the color structure in the soft and hard functions is more
involved, which can be expressed in terms of the $3\times 3$ matrices
representing the space of the operators in the channel. After performing the
computation at next-to-leading order (NLO), we estimate the contributions from
the singlets and the nonsinglets with the resummed results at next-to-leading
logarithmic accuracy (NLL) for these two channels.

The structure of the paper is as follows: In section 2, we construct the
effective operators relevant to the channel
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ in SCET. We choose a basis of
operators proportional to $\delta^{ab}$, $if^{abc}$ and $d^{abc}$, where
$f^{abc}$ and $d^{abc}$ are the structure constants of the SU($N$) group. This
choice of the basis is significant because the symmetric structure constants
$d^{abc}$ vanish in the SU(2) weak interaction, and it is separated from the
beginning. In section 3, the factorization theorem for
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ is established. By including
the previous result of the factorization for
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, the
whole factorization for the muon-pair dijets is completed. We briefly review
the rapidity divergence in SCET, and discuss how to regulate the rapidity
divergence in section 4. Each factorized part is computed at NLO. In section
5, the beam functions and the PDFs for the gauge bosons with the matching
coefficients are presented. For completeness, the muon semi-inclusive jet
functions are quoted from ref. Chay:2021arz . In section 6, the hard functions
for $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ are presented, based on the
result from QCD Kelley:2010fn , and the soft function for this process is
given in section 7. In section 8, all the anomalous dimensions are collected,
and the RG evolution of the $N$-jettiness is discussed. In section 9, we
present a numerical analysis of the 2-jettiness near threshold for the SU(2)
gauge interaction. In section 10, we give a conclusion and an outlook. In
appendix A, the detailed computation of the beam functions for the gauge
bosons and the PDFs at NLO is explained. In appendix B, the tree-level soft
matrices are listed.

##  2 Effective operators for $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$

The effective operators responsible for
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ are given as

$O_{I}^{\mu\alpha\beta}=\overline{\ell}_{3L}T_{I}^{ab}\gamma^{\mu}\ell_{4L}\mathcal{B}_{1\perp}^{a\alpha}\mathcal{B}_{2\perp}^{b\beta},$
(1)

where the collinear-gauge invariant lepton fields $\ell_{n}$ and the gauge
fields $\mathcal{B}_{\perp}^{\mu}$ are defined as

$\ell_{n}(x)=W_{n}^{\dagger}(x)\xi_{n}(x),\ \
\mathcal{B}_{n\perp}^{\mu}(x)=\frac{1}{g}[W_{n}^{\dagger}(x)iD_{n\perp}^{\mu}W_{n}(x)].$
(2)

Here $iD_{n\perp}^{\mu}=\mathcal{P}_{n\perp}^{\mu}+gA_{n\perp}^{\mu}$ is the
covariant derivative. The collinear Wilson line is given as

$W_{n}(x)=\mathrm{P}\exp\Bigl{(}ig\int_{-\infty}^{0}ds\overline{n}\cdot
A_{n}(x+s\overline{n})\Bigr{)}=\sum_{\mathrm{perm.}}\exp\Bigl{[}-g\frac{\overline{n}\cdot
A_{n}(x)}{\overline{n}\cdot\mathcal{P}}\Bigr{]},$ (3)

where P denotes the path ordering along the integration path.

There are three independent operators, and we choose the basis $T_{I}^{ab}$
$(I=1,2,3)$ as

$T_{1}^{ab}=\delta^{ab},\ \ T_{2}^{ab}=if^{abc}t^{c},\ \
T_{3}^{ab}=d^{abc}t^{c},$ (4)

where $t^{a}$ are the SU($N$) generators. Another equivalent basis for $N\geq
3$ can be selected as $\delta^{ab}$, $t^{a}t^{b}$ and $t^{b}t^{a}$. However,
in SU(2) weak interaction, the bases $\delta^{ab}$, $t^{a}t^{b}$ and
$t^{b}t^{a}$ are no longer independent. In other words, $d^{abc}=0$ in SU(2),
and there are only two independent operators. In this case, we disregard
$T_{3}^{ab}$.

The soft interactions are decoupled from the collinear fields by the field
redefinition Bauer:2001yt

$\ell_{n}^{(0)}(x)=Y_{n}^{\dagger}(x)\ell_{n}(x),\ \
\mathcal{B}_{n\perp}^{\mu(0)}(x)=Y_{n}^{\dagger}(x)\mathcal{B}_{n\perp}^{\mu}(x)Y_{n}(x),$
(5)

where the soft Wilson line $Y_{n}(x)$ in the fundamental representation is
given as

$Y_{n}(x)=\mathrm{P}\exp\Bigl{(}ig\int_{-\infty}^{0}dsn\cdot
A_{us}^{c}t^{c}(x+sn)\Bigr{)}=\sum_{\mathrm{perm.}}\exp\Bigl{[}-g\frac{n\cdot
A_{us}(x)}{n\cdot\mathcal{P}}\Bigr{]}.$ (6)

For $\mathcal{B}_{n\perp}^{\mu a}$, we can employ the adjoint representation
$\mathcal{Y}$ for the soft Wilson line Bauer:2001yt  as

$\mathcal{B}_{n\perp}^{\mu a(0)}=\mathcal{Y}_{n}^{ab}\mathcal{B}_{n\perp}^{\mu
b},$ (7)

where the soft Wilson line $\mathcal{Y}$ is obtained by replacing $t^{c}$ by
the generators $\mathcal{T}^{c}$ in the adjoint representation with
$(\mathcal{T}^{c})^{ab}=-if^{cab}$. And the relation between $\mathcal{Y}$ and
$Y$ is given as
$\mathcal{Y}_{n}^{ab}=\mathrm{tr}[Y_{n}^{\dagger}t^{b}Y_{n}t^{a}]$. From now
on, we use the fields after the decoupling and we drop the superscript (0) for
simplicity.

With the field redefinition, the operators in eq. (1) are written as

$O_{I\mu}^{\alpha\beta}=\overline{\ell}_{3L}Y_{3}^{\dagger}\gamma_{\mu}T_{I}^{ab}Y_{4}\ell_{4L}\mathcal{Y}_{1}^{ac}\mathcal{B}_{1\perp}^{c\alpha}\mathcal{Y}_{2}^{bd}\mathcal{B}_{2\perp}^{d\beta}.$
(8)

In terms of these operators, we can write the effective Lagrangian as

$\mathcal{L}_{\mathrm{eff}}=-i\sum_{I}D^{\mu}_{I\alpha\beta}O^{\alpha\beta}_{I\mu}+\textrm{hermitian\
conjugate},$ (9)

where $D^{\mu}_{I\alpha\beta}$ are the Wilson coefficients, which are obtained
by integrating out the degrees of freedom of order $Q$.

We consider the $N$-jettiness, which is defined as Stewart:2010tn ;
Jouttenus:2011wh

$\mathcal{T}_{N}=\sum_{k}\mathrm{min}\ \Bigl{\\{}\frac{2q_{i}\cdot
p_{k}}{\omega_{i}}\Bigr{\\}},$ (10)

where $i$ runs over 1, 2 for the beams, and $3,\cdots,N+2$ for the final-state
jets. Here $q_{i}$ are the reference momenta of the beams and the jets with
the normalization factors $\omega_{i}=\overline{n}_{i}\cdot q_{i}$, and
$p_{k}$ are the momenta of all the measured particles in the final state.

$q_{1,2}^{\mu}=\frac{1}{2}z_{1,2}E_{\mathrm{cm}}n_{1,2}^{\mu}=\frac{1}{2}\omega_{1,2}n_{1,2}^{\mu},\
q_{i}^{\mu}=\frac{1}{2}\omega_{i}n_{i}^{\mu},(i=3,\cdots,N+2).$ (11)

Here $z_{1,2}$ are the momentum fractions of the beams. The lightcone vectors
$n_{1}$ and $n_{2}$ for the beams are aligned to the $z$ direction,
$n_{1}^{\mu}=(1,0,0,1),\ n_{2}^{\mu}=(1,0,0,-1)$, and $n_{i}$ are the
lightcone vectors specifying the jet directions.

The 2-jettiness from the channel
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ can be written as Bauer:2008jx

$\Bigl{(}\frac{d\sigma}{d\mathcal{T}_{2}}\Big{)}_{W}=\frac{1}{2s}\int
d^{4}x\sum_{X}\langle I|\mathcal{L}_{\mathrm{eff}}(x)|X\rangle\langle
X|\mathcal{L}_{\mathrm{eff}}(0)|I\rangle\delta\Bigl{(}\mathcal{T}_{2}-g(I,X)\Bigr{)},$
(12)

where $|I\rangle$ represents the initial state, $|X\rangle$ denotes the final
state, and the sum over $X$ includes the phase space. The function $g(I,X)$
extracts the jettiness from the states $|I\rangle$ and $|X\rangle$. In SCET,
the final states $|X\rangle$ consist of the $n_{i}$-collinear states
$|X_{i}\rangle$ and the soft states $|X_{s}\rangle$. Since the
$n_{i}$-collinear particles do not interact with each other, and the soft
particles are decoupled from the collinear sectors, the final states
$|X\rangle$ in the Hilbert space consist of the tensor product of the
collinear states $|X_{i}\rangle$ and the soft states $|X_{s}\rangle$ as

$|X\rangle=|X_{1}\rangle\otimes|X_{2}\rangle\otimes|X_{3}\rangle\otimes|X_{4}\rangle\otimes|X_{s}\rangle.$
(13)

In order to obtain the 2-jettiness, we multiply $O_{I\mu}^{\alpha\beta}$ and
the hermitian conjugate
$O_{J\mu^{\prime}}^{\dagger\alpha^{\prime}\beta^{\prime}}$ with their
respective Wilson coefficients and implement the measurement of the
2-jettiness. Since the result looks quite complicated, we first consider the
product of the operators, for simplicity. It is written as

$\displaystyle\sum_{IJ}D_{J\alpha^{\prime}\beta^{\prime}}^{\mu^{\prime}*}\Bigl{(}O_{J\mu^{\prime}}^{\alpha^{\prime}\beta^{\prime}}(x)\Bigr{)}^{\dagger}D_{I\alpha\beta}^{\mu}O_{I\mu}^{\alpha\beta}(0)$
$\displaystyle=\sum_{IJ}D_{J\alpha^{\prime}\beta^{\prime}}^{\mu^{\prime}*}D_{I\alpha\beta}^{\mu}\Bigl{[}(\overline{\ell}_{4})_{j_{4}}^{\mu_{4}}\Bigl{(}Y_{3}^{\dagger}T_{J}^{a^{\prime}b^{\prime}}Y_{3}\Bigr{)}_{j_{4}j_{3}}(\gamma_{\mu^{\prime}})_{\mu_{4}\mu_{3}}(\ell_{3})_{j_{3}}^{\mu_{3}}\mathcal{Y}_{1}^{c^{\prime}a^{\prime}}\mathcal{B}_{1\perp}^{c^{\prime}\alpha^{\prime}}\mathcal{Y}_{2}^{b^{\prime}d^{\prime}}\mathcal{B}_{2\perp}^{d^{\prime}\beta^{\prime}}\Bigr{]}(x)$
$\displaystyle\times\Bigl{[}(\overline{\ell}_{3})_{i_{3}}^{\nu_{3}}\Bigl{(}Y_{3}^{\dagger}T_{I}^{ab}Y_{4}\Bigr{)}_{i_{3}i_{4}}(\gamma_{\mu})_{\nu_{3}\nu_{4}}(\ell_{4})_{i_{4}}^{\nu_{4}}\mathcal{Y}_{1}^{ca}\mathcal{B}_{1\perp}^{c\alpha}\mathcal{Y}_{2}^{bd}\mathcal{B}_{2\perp}^{d\beta}\Bigr{]}(0),$
(14)

where $i_{k}$, $j_{k}$ are the gauge indices, while $\mu_{k}$, $\nu_{k}$ are
the Dirac indices. The measurement of the 2-jettiness, $g(I,X)$, will be
implemented after simplifying this expression.

We rearrange the lepton bilinears in each collinear direction as
Manohar:2018kfx ; Chay:2021arz

$\displaystyle\Bigl{(}\overline{\ell}_{n}(x)\Bigr{)}_{\alpha}^{i}\Bigl{(}\ell_{n}(y)\Bigr{)}_{\beta}^{j}$
$\displaystyle=(\mathcal{P}_{L}{n}\\!\\!\\!/)_{\beta\alpha}\Bigl{[}\frac{1}{2N}\delta^{ij}\overline{\ell}_{n}(x)\frac{{\overline{n}}\\!\\!\\!/}{2}\ell_{n}(y)+(t^{c})^{ij}\overline{\ell}_{n}(x)\frac{{\overline{n}}\\!\\!\\!/}{2}t^{c}\ell_{n}(y)\Bigr{]}$
$\displaystyle=(\mathcal{P}_{L}{n}\\!\\!\\!/)_{\beta\alpha}\sum_{c}k_{\ell}^{c}(T^{c})^{ij}C_{\ell}^{c}(x,y),$
(15)

where $\mathcal{P}_{L}=(1-\gamma_{5})/2$ and we express the generators
$(T^{0})^{ij}=\delta^{ij}$ for the singlet and $(T^{c})^{ij}=(t^{c})^{ij}$ for
the nonsinglets with the corresponding collinear operators

$C_{\ell}^{a}(x,y)=\overline{\ell}_{n}(x)\frac{{\overline{n}}\\!\\!\\!/}{2}T^{a}\ell_{n}(y).$
(16)

The factors $k^{a}_{\ell}$ are defined as $k^{0}_{\ell}=1/(2N)$ for the
singlet and $k^{a}_{\ell}=1$ for the nonsinglets.

The product of the collinear gauge fields can be written as

$\displaystyle\mathcal{B}_{n\perp}^{a\alpha}(x)\mathcal{B}_{n\perp}^{b\beta}(y)$
$\displaystyle=\frac{g_{\perp}^{\alpha\beta}}{D-2}\Bigl{[}\frac{\delta^{ab}}{N^{2}-1}\Bigl{(}\delta^{cd}\mathcal{B}_{n\perp}^{c\mu}(x)\mathcal{B}_{n\perp\mu}^{d}(y)\Bigr{)}+\frac{Nd^{abc}}{N^{2}-4}\Bigl{(}d^{dec}\mathcal{B}_{n\perp}^{d\mu}(x)\mathcal{B}_{n\perp\mu}^{e}(y)\Bigr{)}\Bigr{]}$
$\displaystyle\equiv
g_{\perp}^{\alpha\beta}\sum_{m=0}^{N^{2}-1}k_{W}^{m}(G^{m})^{ab}C_{W}^{m}(x,y),$
(17)

where $D$ is the number of the spacetime dimensions. The factors $k_{W}^{m}$
are given as

$k_{W}^{m}=\left\\{\begin{array}[]{ll}\displaystyle\frac{1}{D-2}\frac{1}{N^{2}-1},&m=0,\
(\textrm{singlet}),\\\ \displaystyle\frac{1}{D-2}\frac{N}{N^{2}-4},&\
m=1,\cdots,N^{2}-1,\ (\textrm{nonsinglets}),\end{array}\right.$ (18)

and $(G^{0})^{ab}=\delta^{ab}$ for the singlet and $(G^{m})^{ab}=d^{abm}$ for
the nonsinglets. Note that $d^{abm}=0$ in SU(2), hence there are no nonsinglet
beam functions for the gauge bosons. However, we keep these contributions for
the general SU($N$) gauge theory, and put $N=2$ later in the numerical
analysis. And $C_{W}^{m}(x,y)$ are the operators, which are defined as

$C_{W}^{m}(x,y)=(G^{m})^{de}\mathcal{B}_{n\perp}^{d\mu}(x)B_{n\perp\mu}^{e}(y).$
(19)

Using these relations, eq. (2) can be written as

$\displaystyle\sum_{IJ}D_{J\alpha^{\prime}\beta^{\prime}}^{\mu^{\prime}*}D_{I\alpha\beta}^{\mu}\mathrm{tr}(\mathcal{P}_{L}{n}\\!\\!\\!/_{4}\gamma_{\mu^{\prime}}{n}\\!\\!\\!/_{3}\gamma_{\mu})g_{\perp}^{\alpha\alpha^{\prime}}g_{\perp}^{\beta\beta^{\prime}}\sum_{efgh}k_{W}^{e}k_{W}^{f}k_{\ell}^{g}k_{\ell}^{h}C_{B_{1}}^{e}(x,0)C_{B_{2}}^{f}(0,x)C_{\ell_{3}}^{g}(0,x)C_{\ell_{4}}^{h}(x,0)$
$\displaystyle\times\Bigl{(}\mathcal{Y}_{1}^{c^{\prime}a^{\prime}}\mathcal{Y}_{2}^{b^{\prime}d^{\prime}}(x)\Bigr{)}(G^{e})^{c^{\prime}c}(G^{f})^{dd^{\prime}}\Bigl{(}\mathcal{Y}_{1}^{ca}\mathcal{Y}_{2}^{bd}(0)\Bigr{)}\mathrm{Tr}\,\Bigl{[}\Bigl{(}T^{l}Y_{4}^{\dagger}T_{J}^{a^{\prime}b^{\prime}}Y_{3}(x)\Bigr{)}\Bigl{(}T^{m}Y_{3}^{\dagger}T_{I}^{ab}Y_{4}(0)\Bigr{)}\Bigr{]}.$
(20)

This is the matrix element squared to be employed in expressing the
$N$-jettiness, which is obtained by implementing the measurement of the
$N$-jettiness $g(I,X)$ in eq. (12). In order to establish the factorization,
we define the beam functions and the jet functions as the matrix elements of
the collinear fields. The soft function consists of the soft Wilson lines, and
the hard function is the combination of the Wilson coefficients of the
operators, contracted with the appropriate Dirac structure.

##  3 Factorization of the muon dijet process

Though eq. (2) looks complicated, it can be rearranged to make the
factorization look manifest. The hard function $H_{IJ}$ is obtained by
combining the matching coefficients $D_{I\alpha\beta}^{\mu}$, contracted with
the Dirac structure, and is written as

$H_{IJ}=D_{I\alpha\beta}^{\mu}D_{J\alpha^{\prime}\beta^{\prime}}^{\mu^{\prime}*}g_{\perp}^{\alpha\alpha^{\prime}}g_{\perp}^{\beta\beta^{\prime}}\mathrm{tr}(\mathcal{P}_{L}{n}\\!\\!\\!/_{4}\gamma_{\mu^{\prime}}{n}\\!\\!\\!/_{3}\gamma_{\mu}).$
(21)

The collinear functions are obtained by taking the matrix elements of the
operators $C_{W}^{a}$ between the gauge bosons to yield the beam functions,
and $C_{\ell}^{a}$ between the vacuum to yield the jet function after
inserting the measurement of the jettiness. The beam functions
$B_{W}^{a}(t,z,\mu,\nu)$ are defined as

$B_{W}^{a}(t,z=\omega/P^{-},\mu,\nu)=-\omega\langle
W|B_{n\perp\mu}^{c}\delta(t+\omega
n\cdot\mathcal{P})[\delta(\omega-\overline{n}\cdot\mathcal{P})B_{n\perp}^{d\mu}](G^{a})^{cd}|W\rangle,$
(22)

where $P^{-}=\overline{n}\cdot P$, and $t=\omega\mathcal{T}$. At tree level,
the beam functions are normalized as
$B_{W}^{a(0)}(t,z,\mu,\nu)=\delta(t)\delta(1-z)(G^{a})$.

As introduced in ref. Chay:2021arz , we define the semi-inclusive jet function
as

$\displaystyle J_{l}^{a}(p^{2},\mu,\nu)$
$\displaystyle=\int\frac{d\overline{n}\cdot p_{l}}{\overline{n}\cdot
p_{l}}\int\frac{d^{2}\mathbf{p}_{l}^{\perp}}{\omega}\sum_{X}\mathrm{Tr}\langle
0|\ell_{L}(0)\delta(n\cdot
p+n\cdot\mathcal{P})\delta^{(2)}(p_{\perp}+\mathcal{P}_{\perp})|lX\rangle$
$\displaystyle\times\langle
lX|[\delta(\omega+\overline{n}\cdot\mathcal{P})\overline{\ell}_{L}(0)]\frac{{\overline{n}}\\!\\!\\!/}{2}T^{a}|0\rangle,$
(23)

where the lepton $l$ (muon or antimuon) is specified in the final state. It is
normalized at tree level as
$J_{l}^{a(0)}(p^{2})=\delta(p^{2})\mathrm{Tr}(T^{a}P_{l})$, where $P_{l}$ is
the projection operator to the given lepton $l$. For example, the projection
operator for the muon is given by $P_{\mu}=(1-\sigma^{3})/2$ in SU(2), and
$P_{\nu_{\mu}}=(1+\sigma^{3})/2$ for the muon neutrino.

The soft functions are given by the vacuum expectation values of the soft part
in eq. (2), and are given by

$\displaystyle S_{JI}^{efgh}(\mathcal{T}_{s},\mu,\nu)$ $\displaystyle=\langle
0|\mathrm{Tr}\,\Bigl{[}\Bigl{(}T^{g}Y_{4}^{\dagger}T_{J}^{a^{\prime}b^{\prime}}Y_{3}\mathcal{Y}_{1}^{c^{\prime}a^{\prime}}\mathcal{Y}_{2}^{b^{\prime}d^{\prime}}(G^{e})^{c^{\prime}c}\Bigr{)}(x)$
$\displaystyle\times\delta\Bigl{(}\mathcal{T}_{s}-\sum_{X_{s}}\mathrm{min}(\\{n_{i}\cdot
p_{X_{s}}\\})\Bigr{)}\Bigl{(}(G^{f})^{dd^{\prime}}\mathcal{Y}_{1}^{ca}\mathcal{Y}_{2}^{bd}T^{h}Y_{3}^{\dagger}T_{I}^{ab}Y_{4}\Bigr{)}(0)\Bigr{]}.$
(24)

The detailed procedure of the factorization is delineated in detail in ref.
Chay:2021arz . Here we present the result following the same procedure in the
scheme of $\mathrm{SCET_{I}}$. The factorization of the 2-jettiness from
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ is written as

$\displaystyle\Bigl{(}\frac{d\sigma}{d\mathcal{T}_{2}}\Bigr{)}_{W}$
$\displaystyle=\frac{1}{2s}\int\frac{dz_{1}}{z_{1}}\int\frac{dz_{2}}{z_{2}}\int
dt_{1}dt_{2}\int\frac{d^{4}p_{3}}{(2\pi)^{4}}\frac{d^{4}p_{4}}{(2\pi)^{4}}(2\pi)^{4}\delta^{(4)}\Bigl{(}\frac{\omega_{1}n_{1}}{2}+\frac{\omega_{2}n_{2}}{2}-\frac{\omega_{3}n_{3}}{2}-\frac{\omega_{4}n_{4}}{2}\Bigr{)}$
$\displaystyle\times\int
d\mathcal{T}_{s}\delta\Bigl{(}\mathcal{T}_{2}-\frac{t_{1}}{\omega_{1}}-\frac{t_{2}}{\omega_{2}}-\frac{p_{3}^{2}}{\omega_{3}}-\frac{p_{4}^{2}}{\omega_{4}}-\mathcal{T}_{s}\Bigr{)}\sum_{IJ}H_{IJ}S_{JI}^{efgh}(\mathcal{T}_{s},\mu,\nu)\sum_{efgh}k_{W}^{e}k_{W}^{f}k_{\ell}^{g}k_{\ell}^{h}$
$\displaystyle\times
B_{1}^{e}(t_{1},z_{1},\mu,\nu)B_{2}^{f}(t_{2},z_{2},\mu,\nu)J_{3}^{g}(p_{3}^{2},\mu,\nu)J_{4}^{h}(p_{4}^{2},\mu,\nu).$
(25)

Combining all these ingredients, the 2-jettiness for the final
$\mu^{-}\mu^{+}$ is obtained by adding the contribution from the initial
particles $WW$ and $\ell_{e}\overline{\ell}_{e}$, and is factorized as

$\displaystyle\Bigl{(}\frac{d\sigma}{d\mathcal{T}_{2}}\Bigr{)}$
$\displaystyle=\frac{1}{2s}\int\frac{dz_{1}}{z_{1}}\int\frac{dz_{2}}{z_{2}}\int
dt_{1}dt_{2}\int\frac{d^{4}p_{3}}{(2\pi)^{3}}\frac{d^{4}p_{4}}{(2\pi)^{3}}(2\pi)^{4}\delta^{(4)}\Bigl{(}\frac{\omega_{1}n_{1}}{2}+\frac{\omega_{2}n_{2}}{2}-\frac{\omega_{3}n_{3}}{2}-\frac{\omega_{4}n_{4}}{2}\Bigr{)}$
$\displaystyle\times\int
d\mathcal{T}_{s}\delta\Bigl{(}\mathcal{T}_{2}-\frac{t_{1}}{\omega_{1}}-\frac{t_{2}}{\omega_{2}}-\frac{p_{3}^{2}}{\omega_{3}}-\frac{p_{4}^{2}}{\omega_{4}}-\mathcal{T}_{s}\Bigr{)}\sum_{i=W,e}\sum_{IJ}H_{IJ}^{i}S_{JI}^{i,efgh}(\mathcal{T}_{s},\mu,\nu)$
$\displaystyle\times\sum_{efgh}k_{i}^{e}k_{i}^{f}k_{\mu}^{g}k_{\mu}^{h}B_{i}^{e}(t_{1},z_{1},\mu,\nu)B_{\bar{i}}^{e}(t_{2},z_{2},\mu,\nu)J_{\mu}^{g}(p_{3}^{2},\mu,\nu)J_{\bar{\mu}}^{h}(p_{4}^{2},\mu,\nu),$
(26)

and the corresponding hard and soft functions in each process should be
employed.111Because the contraction of the Dirac structure is different in two
channels, we use $H^{e}=16H$, where $H$ is the hard function for
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ in
ref. Chay:2021arz .

##  4 Rapidity divergence

We briefly explain the rapidity divergence and the method of extracting it by
prescribing the rapidity regulators. The rapidity divergence arises because
the collinear and the soft modes with the same invariant mass reside in
disparate phase spaces. These modes lie on the same hyperbola in the
$p^{+}$-$p^{-}$ plane, and they are distinguished only by their rapidities.
When the soft particles reach the collinear region, the momentum approaches
the region $p^{+}p^{-}\rightarrow\infty$, with $p_{\perp}^{2}$ fixed. In this
limit, an additional divergence is induced, which has different source from
the UV and the IR divergences, hence the rapidity divergence cannot be handled
with the dimensional regularization. The converse is also true. When the
collinear particles reach the soft region, the zero-bin subtraction should be
performed to avoid double counting. And the zero-bin contribution induces the
rapidity divergence. However, no rapidity divergence shows up in the full
theory because there is no separation of the phase space. It means that there
may be rapidity divergence in the collinear and the soft sectors in SCET, but
it is cancelled in the total contribution.

In spite of the cancellation, we need to regulate the rapidity divergence in
each sector. There have been many methods to regulate the rapidity divergence
collins_2011 ; Idilbi:2007ff ; Idilbi:2007yi ; Becher:2011dz ; Chiu:2011qc ;
Chiu:2012ir ; Li:2016axz ; Ebert:2018gsn . Recently, one of the authors has
constructed the soft and collinear rapidity regulators in ref. Chay:2020jzn ,
which correctly yields the angular dependence when different lightcone
directions are involved. We will adopt that prescription to compute the
rapidity divergence here. Let us describe how to prescribe the rapidity
regulators in the collinear and the soft sectors.

In the $n$-collinear sector, the basic idea is to attach a rapidity regulator
of the form $(\nu/\overline{n}\cdot k)^{\eta}$ in the $n$-collinear Wilson
line, where the rapidity divergence arises when $\overline{n}\cdot
k\rightarrow\infty$ with $\mathbf{k}_{\perp}^{2}$ fixed. This prescription is
employed in ref. Chiu:2012ir , in which the collinear rapidity divergence
appears as poles of $\eta$, and $\nu$ is the rapidity scale. The rapidity
divergence is manipulated by modifying the region $\overline{n}\cdot
k\rightarrow\infty$. However, it also modifies the region $\overline{n}\cdot
k\rightarrow 0$, which is unwanted in the rapidity evolution. The undesirable
divergence in this region is cancelled when the zero-bin subtraction is
performed.

The construction of the soft rapidity regulator is different from that in ref.
Chiu:2012ir , and our choice yields the correct angular dependence of the
lightlike directions. In terms of Feynman diagrams, the $n$-collinear Wilson
line is obtained by summing over all the emissions of $n$-collinear gluons
from the sources other than the $n$-collinear particle. For example, if we
consider the back-to-back collinear current
$\overline{\xi}_{\bar{n}}W_{\bar{n}}Y^{\dagger}_{\bar{n}}\gamma^{\mu}Y_{n}W^{\dagger}_{n}\xi_{n}$,
the $n$-collinear Wilson line is obtained by considering the emission of the
$n$-collinear gluons from the $\overline{n}$-collinear source, and it is
exponentiated to produce the $n$-collinear Wilson line.

On the contrary, the soft Wilson line $Y_{n}$ is obtained by considering the
emission of soft gluons from the $n$-collinear source. Note that the sources
of the emitted gluons to construct the Wilson lines are different in the
collinear case ($W_{n}$) and in the soft case ($Y_{n}$). If we require the
consistency to treat the rapidity divergence in the collinear and the soft
sectors, we should track the same source for the rapidity divergence.
Therefore, the soft rapidity regulator in the soft Wilson line, corresponding
to the $n$-collinear gluons in $W_{n}$, should be attached to the
$\overline{n}$-collinear source, that is, $Y_{\bar{n}}$. By generalizing the
procedure, let us consider the collinear current
$\overline{\xi}_{n_{1}}W_{n_{1}}Y_{n_{1}}^{\dagger}\Gamma
Y_{n_{2}}W_{n_{2}}^{\dagger}\xi_{n_{2}}$, which is invariant under collinear
and soft gauge transformations respectively. The modified collinear Wilson
line $W_{n_{1}}$ and the soft Wilson line $Y_{n_{2}}$, with the rapidity
regulators are given as

$\displaystyle W_{n_{1}}$
$\displaystyle=\sum_{\mathrm{perm.}}\exp\Bigl{[}-\frac{g}{\overline{n}_{1}\cdot\mathcal{P}}\Bigl{(}\frac{\nu}{|\overline{n}_{1}\cdot\mathcal{P}|}\Bigr{)}^{\eta}\overline{n}_{1}\cdot
A_{n_{1}}\Bigr{]},$ $\displaystyle Y_{n_{2}}$
$\displaystyle=\sum_{\mathrm{perm.}}\exp\Bigl{[}-\frac{g}{n_{2}\cdot\mathcal{P}}\Bigl{(}\frac{\nu}{|n_{2}\cdot\mathcal{P}|}\frac{n_{1}\cdot
n_{2}}{2}\Bigr{)}^{\eta}n_{2}\cdot A_{us}\Bigr{]},$ (27)

where $\mathcal{P}$ is the operator extracting the momentum. The remaining
Wilson lines $W_{n_{2}}$ and $Y_{n_{1}}$ can be obtained by switching $n_{1}$
and $n_{2}$. The point in selecting the rapidity regulator is to trace the
same emitted gauge bosons both in the collinear and the soft sectors, which
are eikonalized to produce the Wilson lines. Note that the rapidity
divergences from $W_{n_{1}}$ and $Y_{n_{2}}$ have the same origin because the
collinear and soft gauge bosons are emitted from the $n_{2}$-collinear quark
for both of the Wilson lines.

When the rapidity divergence arises in the soft sector, the soft momentum $k$,
in the limit $\overline{n}_{1}\cdot k\rightarrow\infty$ approaches
$k^{\mu}\approx(\overline{n}_{1}\cdot k)n_{1}^{\mu}/2$ and the soft rapidity
regulator becomes

$\Bigl{(}\frac{\nu}{n_{2}\cdot k}\frac{n_{1}\cdot
n_{2}}{2}\Bigr{)}^{\eta}\xrightarrow[\bar{n}_{1}\cdot
k\to\infty]{}\Bigl{(}\frac{\nu}{\overline{n}_{1}\cdot k}\Bigr{)}^{\eta},$ (28)

which has the same form as the collinear rapidity regulator for $W_{n_{1}}$.
Another pair possessing the same source of rapidity divergence is
$W_{n_{2}}^{\dagger}$ and $Y_{n_{1}}^{\dagger}$.

##  5 Beam functions, PDF and semi-inclusive jet functions

The bare operator for the gauge-boson beam function is given as

$\mathcal{O}_{W}^{a,\mathrm{bare}}(t,\omega)=-\omega\theta(\omega)B_{n\perp\mu}^{c}(0)\delta(t+\omega
n\cdot\mathcal{P})[\delta(\omega-\overline{n}\cdot\mathcal{P})B_{n\perp}^{d\mu}(0)](G^{a})^{cd}.$
(29)

This operator is extended from QCD Berger:2010xi  to include the nonsinglets.
The beam function for the gauge boson is given by the matrix element of the
renormalized operator $\mathcal{O}_{W}^{a}(t,\omega,\mu,\nu)$,

$\displaystyle B_{W}^{a}(t,x=\omega/P^{-},\mu,\nu)$ $\displaystyle=\langle
W(P^{-})|\mathcal{O}_{W}^{a}(t,\omega,\mu,\nu)|W(P^{-})\rangle$
$\displaystyle=-\omega\langle W|B_{n\perp\mu}^{c}\delta(t+\omega
n\cdot\mathcal{P})[\delta(\omega-\overline{n}\cdot\mathcal{P})B_{n\perp}^{d\mu}](G^{a})^{cd}|W\rangle.$
(30)

Note that there is no nonsinglet beam function for the SU(2) gauge group in
the weak interaction because $d^{abc}=0$. However, we will proceed for the
general SU($N$) here. The definition of the singlet beam function for the
gauge boson is the same as the gluon beam function in QCD Berger:2010xi .
However, the authors in ref. Berger:2010xi  have computed the gluon beam
function with massless gluons, with no rapidity regulator. We compute the
radiative corrections of the beam function at NLO with the nonzero gauge boson
mass $M$, and the rapidity regulators.

The detailed computation of the beam functions at NLO is relegated to appendix
A, and we present only the results here. The singlet and the nonsinglet
functions in eq. (5) are denoted as $B_{Ws}(G^{0})^{ab}$ and
$B_{Wn}(G^{a})^{bc}$ respectively by extracting the group theory factors. The
singlet beam function for the gauge boson at NLO is given as

$\displaystyle
B_{Ws}^{(1)}(t,x,\mu)=\frac{\alpha}{2\pi}\Bigl{\\{}\delta(1-x)\Bigl{[}C_{A}\Bigl{(}\frac{2}{\epsilon^{2}}\delta(t)-\frac{2}{\epsilon}\frac{1}{\mu^{2}}\mathcal{L}_{0}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}\Bigr{)}+\frac{\beta_{0}}{2}\frac{1}{\epsilon}\Bigr{]}$
$\displaystyle+\delta(t)\Bigl{(}\frac{\beta_{0}}{2}\delta(1-x)+C_{A}P_{WW}(x)\Bigr{)}\ln\frac{\mu^{2}}{M^{2}}+C_{A}\delta(t)\Bigl{(}P_{WW}(x)\ln\frac{1-x}{x}-\frac{\pi^{2}}{6}\delta(1-x)\Bigr{)}$
$\displaystyle+C_{A}\Bigl{[}\delta(1-x)\frac{2}{\mu^{2}}\mathcal{L}_{1}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}+P_{WW}(x)\frac{1}{\mu^{2}}\mathcal{L}_{0}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}\Bigr{]}$
$\displaystyle+C_{A}\delta(t)\Bigl{[}\delta(1-x)\Bigl{(}\frac{31}{18}-\frac{\pi^{2}}{9}-\frac{\pi}{2\sqrt{3}}\Bigr{)}-\Bigl{(}\frac{2(1-x)}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}\Bigr{)}$
$\displaystyle-
P_{WW}(x)\ln(1-x+x^{2})\Bigr{]}+\delta(1-x)\delta(t)\Bigl{[}-\frac{2}{9}n_{f}T_{F}+\frac{1}{2}n_{s}T_{F}\Bigl{(}-\frac{17}{9}+\frac{\pi}{\sqrt{3}}\Bigr{)}\Bigr{]}\Bigr{\\}},$
(31)

where $M$ is the mass of the gauge boson. Compared to QCD, the first three
lines in eq. (5) are the same as those from QCD with the appropriate color
factors after replacing $\ln\mu^{2}/M^{2}$ by the IR pole
$1/\epsilon_{\mathrm{IR}}$. The last two lines are the additional
contributions due to the nonzero $M$. Here $P_{WW}(x)$ is the splitting
function for $W\rightarrow WW$, which is the same as the splitting function
$P_{gg}(x)$ for $g\rightarrow gg$, and is given by

$\displaystyle P_{WW}(x)$
$\displaystyle=2\mathcal{L}_{0}(1-x)x+2\theta(1-x)\Bigl{[}\frac{1-x}{x}+x(1-x)\Bigr{]}$
$\displaystyle=2\theta(1-x)\Bigl{[}\frac{x}{(1-x)_{+}}+\frac{1-x}{x}+x(1-x)\Bigr{]},$
(32)

and the relation $2(1-x+x^{2})^{2}/x=(1-x)P_{WW}(x)$ is used. The first term
$\beta_{0}$ in the beta function is given by

$\beta_{0}=\frac{11}{3}C_{A}-\frac{4}{3}\Bigl{(}\frac{1}{2}\Bigr{)}n_{f}T_{f}-\frac{1}{3}n_{s}T_{f},$
(33)

where $n_{s}$ is the number of complex scalar multiplets in the theory222In
the standard model with the SU(2) weak interaction, there is one scalar
multiplet, but in SU($N$) gauge theory we need more multiplets to attain equal
masses $M$ of the gauge bosons. We will not dwell on the detailed model
building to achieve this here., and the factor 1/2 in front of $n_{f}$ evokes
the fact that only the left-handed fields contribute here.

The nonsinglet beam function for the gauge boson at NLO is given as

$\displaystyle B_{Wn}^{(1)}(t,x,\mu,\nu)=B_{Ws}^{(1)}(t,x,\mu)$
$\displaystyle+\frac{\alpha
C_{A}}{4\pi}\Bigl{\\{}2\delta(1-x)\Bigl{[}\delta(t)\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}\Bigl{(}\frac{1}{\eta}+\ln\frac{\nu}{p^{-}}\Bigr{)}-\frac{1}{\epsilon^{2}}\delta(t)+\frac{1}{\epsilon}\frac{1}{\mu^{2}}\mathcal{L}_{0}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}-\frac{1}{\mu^{2}}\mathcal{L}_{1}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}\Bigr{]}$
$\displaystyle-
P_{WW}(x)\Bigl{(}\delta(t)\ln\frac{\mu^{2}}{M^{2}}+\frac{1}{\mu^{2}}\mathcal{L}_{0}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}\Bigr{)}+\delta(t)\Bigl{(}P_{WW}(x)\ln\frac{x(1-x+x^{2})}{1-x}+\frac{\pi^{2}}{6}\delta(1-x)$
$\displaystyle+2\frac{1-x}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}\Bigr{)}\Bigr{\\}}.$
(34)

The result for the nonsinglet beam function in eq. (5) is new. Note that the
nonsinglet beam function contains the rapidity divergence.

The bare gauge-boson PDF operator is defined as

$\mathcal{Q}_{W}^{a,\mathrm{bare}}(\omega)=-\omega\theta(\omega)B_{n\perp\mu}^{c}(0)[\delta(\omega-\overline{n}\cdot\mathcal{P})B_{n\perp}^{d\mu}(0)](G^{a})^{cd},$
(35)

and the PDF for the gauge boson is given by the matrix element as

$f_{W}^{a}(x=\omega/P^{-},\mu,\nu)=\langle
W(P^{-})|\mathcal{Q}_{W}^{a}(\omega,\mu,\nu)|W(P^{-})\rangle.$ (36)

By performing the operator-product expansion of
$\mathcal{O}_{W}^{a}(t,\omega,\mu,\nu)$, it can be expressed in terms of the
operators $\mathcal{Q}_{W}^{b}(\omega)$ with the matching coefficients
$\mathcal{I}_{ij}^{ab}$ as

$\mathcal{O}_{W}^{a}(t,\omega,\mu,\nu)=\sum_{j,b}\int\frac{d\omega^{\prime}}{\omega^{\prime}}\mathcal{I}_{Wj}^{ab}\Bigl{(}t,\frac{\omega}{\omega^{\prime}},\mu\Bigr{)}\mathcal{Q}_{j}^{b}(\omega^{\prime},\mu,\nu),$
(37)

to leading order in SCET. By taking the matrix elements, we obtain the
relation between the beam function and the PDF as

$B_{W}^{a}(t,x,\mu,\nu)=\sum_{j,b}\int_{x}^{1}\frac{dx^{\prime}}{x^{\prime}}\mathcal{I}_{Wj}^{ab}\Bigl{(}t,\frac{x}{x^{\prime}},\mu\Bigr{)}f_{j}^{b}(x^{\prime},\mu,\nu).$
(38)

The matching coefficients $\mathcal{I}_{ij}^{ab}$ describe the collinear
initial-state radiation and can be computed perturbatively. Note that they are
independent of the rapidity scale $\nu$. Here $i$, $j$ are the indices for
particle species, and $a$, $b$ are the weak indices.

The singlet and the nonsinglet PDFs for the gauge bosons are denoted as
$f_{Ws}(G^{0})^{ab}$ and $f_{Wn}(G^{a})^{bc}$ respectively. At NLO, the
singlet PDF is given as

$\displaystyle f_{Ws}^{(1)}$
$\displaystyle=\frac{\alpha}{2\pi}\Bigl{\\{}\Bigl{(}C_{A}P_{WW}(x)+\frac{\beta_{0}}{2}\delta(1-x)\Bigr{)}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}$
$\displaystyle-
C_{A}\Bigl{(}P_{WW}(x)\ln(1-x+x^{2})+2\frac{1-x}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}\Bigr{)}$
$\displaystyle+\delta(1-x)\Bigl{[}C_{A}\Bigl{(}\frac{31}{18}-\frac{\pi}{2\sqrt{3}}-\frac{\pi^{2}}{9}\Bigr{)}-\frac{2}{9}n_{f}T_{F}+\frac{1}{2}n_{s}T_{F}\Bigl{(}-\frac{17}{9}+\frac{\pi}{\sqrt{3}}\Bigr{)}\Bigr{\\}}.$
(39)

The matching coefficient $\mathcal{I}_{WW}^{s(1)}$ is obtained by comparing
eqs. (5) and (5) as

$\displaystyle\mathcal{I}_{WW}^{s(1)}$ $\displaystyle=\frac{\alpha
C_{A}}{2\pi}\Bigl{\\{}\frac{2}{\mu^{2}}\mathcal{L}_{1}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\delta(1-x)+\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}P_{WW}(x)$
$\displaystyle+\delta(t)\Bigl{[}P_{WW}(x)\ln\frac{1-x}{x}-\frac{\pi^{2}}{6}\delta(1-x)\Bigr{]}\Bigr{\\}}.$
(40)

It is the same as the matching coefficient for the gluon case in QCD except
the group theory factor. (See eq. (2.14) in ref. Berger:2010xi .)

The nonsinglet PDF for the gauge boson at NLO is given as

$\displaystyle f_{W}^{n(1)}$
$\displaystyle=f_{W}^{s(1)}+\frac{\alpha}{2\pi}\frac{C_{A}}{2}\Bigl{[}\delta(1-x)\Bigl{(}\frac{2}{\eta}+2\ln\frac{\nu}{p^{-}}\Bigr{)}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}$
$\displaystyle-\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}P_{WW}(x)+2\frac{1-x}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}$
$\displaystyle+P_{WW}(x)\ln(1-x+x^{2})\Bigr{]}.$ (41)

The matching coefficient for the nonsinglet $\mathcal{I}_{WW}^{n(1)}$ is
proportional to the singlet matching coefficient $\mathcal{I}_{WW}^{s(1)}$,
and its relation is given by

$\mathcal{I}_{WW}^{n(1)}=\frac{1}{2}\mathcal{I}_{WW}^{s(1)}.$ (42)

It is interesting to note that the matching coefficients for the nonsingets
are proportional to those for the singlets both for the gauge bosons and for
the leptons though the proportionality constant is different.
$I_{\ell\ell}^{n(1)}=-I_{\ell\ell}^{s(1)}/(N^{2}-1)$. (See ref. Chay:2021arz
.)

For completeness, we present the lepton beam functions from ref. Chay:2021arz
. We express the singlet and nonsinglet beam functions $B_{\ell}^{0}$ and
$B_{\ell}^{a}$ in terms of $B_{s}$ and $B_{n}$ by extracting and separating
the group theory factors as

$B^{0}_{\ell}(t,x,M,\mu)=B_{s}(t,x,M,\mu)\mathrm{Tr}(T^{0}P_{\ell}),\
B^{a}_{\ell}(t,x,M,\mu)=B_{n}(t,x,M,\mu)\mathrm{Tr}(T^{a}P_{\ell}).$ (43)

The bare beam functions $B_{s}$ and $B_{n}$ are given at NLO as

$\displaystyle B_{s}^{(1)}(t,x,M,\mu)$ $\displaystyle=\frac{\alpha
C_{F}}{2\pi}\Bigl{\\{}\delta(t)\delta(1-x)\Bigl{(}\frac{2}{\epsilon^{2}}+\frac{3}{2\epsilon}+\frac{9}{4}-\frac{\pi^{2}}{2}\Bigr{)}$
(44)
$\displaystyle+\delta(t)\Bigl{[}P_{\ell\ell}(x)\ln\frac{\mu^{2}}{x^{2}M^{2}}+(1+x^{2})\mathcal{L}_{1}(1-x)-(1-x)\theta(x)\theta(1-x)\Bigr{]}$
$\displaystyle+\delta(1-x)\Bigl{[}-\frac{2}{\epsilon}\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}+\frac{2}{\mu^{2}}\mathcal{L}_{1}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\Bigr{]}+(1+x^{2})\mathcal{L}_{0}(1-x)\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\Bigr{\\}},$
$\displaystyle B_{n}^{(1)}(t,x,M,\mu)$ $\displaystyle=B_{s}(t,x,M,\mu)$
$\displaystyle-\frac{\alpha
C_{A}}{4\pi}\Bigl{\\{}-2\delta(t)\delta(1-x)\Bigl{[}\Bigl{(}\frac{1}{\eta}+\ln\frac{\nu}{\omega}\Bigr{)}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}-\frac{1}{\epsilon^{2}}+\frac{\pi^{2}}{12}\Bigr{]}$
$\displaystyle+\delta(t)\Bigl{[}(1+x^{2})\mathcal{L}_{0}(1-x)\ln\frac{\mu^{2}}{x^{2}M^{2}}+(1+x^{2})\mathcal{L}_{1}(1-x)-(1-x)\theta(x)\theta(1-x)\Bigr{]}$
$\displaystyle+\delta(1-x)\Bigl{[}-\frac{2}{\epsilon}\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}+\frac{2}{\mu^{2}}\mathcal{L}_{1}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\Bigr{]}+(1+x^{2})\mathcal{L}_{0}(1-x)\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\Bigr{\\}}.$

Here the splitting function $P_{\ell\ell}(x)$ for $\ell\rightarrow\ell W$ is
the same as the quark splitting function $P_{qq}(x)$, and is given by

$P_{\ell\ell}(x)=P_{qq}(x)=\mathcal{L}_{0}(1-x)(1+x^{2})+\frac{3}{2}\delta(1-x)=\Bigl{[}\theta(1-x)\frac{1+x^{2}}{1-x}\Bigr{]}_{+}.$
(45)

The bare singlet and nonsinglet semi-inclusive jet functions $J_{l}^{0}$ and
$J_{l}^{a}$, with the lepton $l$ in the final state, are given as Chay:2021arz

$J_{l}^{0}(p^{2},M,\mu)=J_{s}(p^{2},M,\mu)\mathrm{Tr}(P_{l}T^{0}),\
J_{l}^{a}(p^{2},M,\mu)=J_{n}(p^{2},M,\mu)\mathrm{Tr}(P_{l}T^{a}),$ (46)

where $P_{l}$ is the projection operator to the lepton $l$. The bare singlet
and the nonsinglet jet functions at NLO are given as

$\displaystyle J_{s}^{(1)}(p^{2},M,\mu)$ $\displaystyle=\frac{\alpha
C_{F}}{2\pi}\Bigl{[}\delta(p^{2})\Bigl{(}\frac{2}{\epsilon^{2}}+\frac{3}{2\epsilon}+\frac{7}{2}-\frac{\pi^{2}}{2}\Bigr{)}-\Bigl{(}\frac{2}{\epsilon}+\frac{3}{2}\Bigr{)}\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{p^{2}}{\mu^{2}}\Bigr{)}+\frac{2}{\mu^{2}}\mathcal{L}_{1}\Bigl{(}\frac{p^{2}}{\mu^{2}}\Bigr{)}\Bigr{]},$
$\displaystyle J_{n}^{(1)}(p^{2},M,\mu)$
$\displaystyle=J_{s}(p^{2},M,\mu)+\frac{\alpha
C_{A}}{2\pi}\Bigl{\\{}\delta(p^{2})\Bigl{[}\frac{1}{\eta}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}-\frac{1}{\epsilon^{2}}+\frac{1}{\epsilon}\ln\frac{\nu}{\omega}+\ln\frac{\mu^{2}}{M^{2}}\ln\frac{\nu}{\omega}$
$\displaystyle+\frac{3}{4}\ln\frac{\mu^{2}}{M^{2}}-\frac{5}{8}+\frac{\pi^{2}}{12}\Bigr{]}+\Bigl{(}\frac{1}{\epsilon}+\frac{3}{4}\Bigr{)}\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{p^{2}}{\mu^{2}}\Bigr{)}-\frac{1}{\mu^{2}}\mathcal{L}_{1}\Bigl{(}\frac{p^{2}}{\mu^{2}}\Bigr{)}\Bigr{\\}}.$
(47)

Like the beam functions, the nonsinglet jet function develops the rapidity
divergence.

##  6 Hard function

The hard functions can be read off from those in QCD Kelley:2010fn . However,
the bases in ref. Kelley:2010fn  are $T^{\prime}_{1}=t^{a}t^{b}$,
$T^{\prime}_{2}=t^{b}t^{a}$ and $T^{\prime}_{3}=\delta^{ab}$, while our bases
consist of $T_{1}=\delta^{ab}$, $T_{2}=if^{abc}t^{c}$ and
$T_{3}=d^{abc}t^{c}$. The change of basis can be obtained by noting the
relation

$t^{a}t^{b}=\frac{1}{2}\Bigl{[}\frac{1}{N}\delta^{ab}+(if^{abc}+d^{abc})t^{c}\Bigr{]}.$
(48)

Therefore the relation between $T^{\prime}_{i}$ and $T_{i}$ is given by

$(T^{\prime}_{1}\ T^{\prime}_{2}\ T^{\prime}_{3})=(T_{1}\ T_{2}\
T_{3})\mathbf{A},$ (49)

where the transformation matrix $\mathbf{A}$ is given by

$\mathbf{A}=\begin{pmatrix}1/2N&1/2N&1\\\ 1/2&\displaystyle-1/2&0\\\
1/2&1/2&0\end{pmatrix}.$ (50)

Then the sum of the operators $D^{\prime}_{I}O^{\prime}_{I}$ of the effective
Lagrangian in eq. (9) has the relation
$O^{\prime}_{I}D^{\prime}_{I}=O_{I}A_{IJ}D^{\prime}_{J}=O_{I}D_{I}$, from
which we obtain the relation $D_{I}=A_{IJ}D^{\prime}_{J}$. Since the hard
coefficients $H_{IJ}$ are proportional to $D_{I}D^{*}_{J}$,
$\mathbf{H}=\mathbf{A}\mathbf{H}^{\prime}\mathbf{A}^{\dagger}$.

In addition to the change of the basis, the appropriately adjusted color
factors for SU(2), and the fact that only the left-handed leptons contribute
to our hard functions should be implemented. With this in mind, the hard
function can be written as

$H_{IJ}(s,t,u)=4g^{4}H_{IJ}^{(0)}+8g^{4}\frac{\alpha}{4\pi}H_{IJ}^{(1)}+\cdots,$
(51)

where $H_{IJ}^{(0)}$ is the hard function at leading order (LO), and
$H_{IJ}^{(1)}$ at NLO. The Mandelstam variables $s$, $t$, $u$ are given by
$s=(p_{1}+p_{2})^{2}=(p_{3}+p_{4})^{2}$,
$t=(p_{1}-p_{3})^{2}=(p_{2}-p_{4})^{2}$,
$u=(p_{1}-p_{4})^{2}=(p_{2}-p_{3})^{2}$, where $p_{i}$ are the partonic
momenta. The variable $u$ is given by $u=-\omega_{1}\omega_{4}n_{1}\cdot
n_{4}/2=-\omega_{2}\omega_{3}n_{2}\cdot n_{3}/2$, and $s$, $t$ can be
expressed accordingly.

The LO hard coeffcients $H_{IJ,\mathrm{QCD}}^{\prime(0)}$ from ref.
Kelley:2010fn  are given as

$\mathbf{H}^{\prime(0)}=\frac{1}{s^{2}}\begin{pmatrix}\displaystyle\frac{u}{t}(t^{2}+u^{2})&\displaystyle
t^{2}+u^{2}&0\\\ \displaystyle
t^{2}+u^{2}&\displaystyle\frac{t}{u}(t^{2}+u^{2})&0\\\ 0&0&0\end{pmatrix},$
(52)

and in our basis it becomes

$\mathbf{H}_{W}^{(0)}=\frac{1}{4s^{2}tu}\begin{pmatrix}(t+u)^{2}(t^{2}+u^{2})/N^{2}&(u^{4}-t^{4})/N&(t+u)^{2}(t^{2}+u^{2})/N\\\
(u^{4}-t^{4})/N&(t-u)^{2}(t^{2}+u^{2})&u^{4}-t^{4}\\\
(t+u)^{2}(t^{2}+u^{2})/N&u^{4}-t^{4}&(t+u)^{2}(t^{2}+u^{2})\end{pmatrix}.$
(53)

The NLO hard functions $H^{\prime(1)}_{IJ}$ are given as

$\displaystyle H_{11}^{\prime(1)}$
$\displaystyle=\frac{u(t^{2}+u^{2})}{ts^{2}}\Bigl{[}-(C_{A}+C_{F})(L(s))^{2}+V_{1}(s,t,u)L(s)\Bigr{]}+\frac{tu}{s^{2}}W_{1}(s,t,u)+\frac{u^{3}}{ts^{2}}W_{2}(s,t,u),$
$\displaystyle H_{22}^{\prime(1)}$
$\displaystyle=\frac{t(t^{2}+u^{2})}{us^{2}}\Bigl{[}-(C_{A}+C_{F})(L(s))^{2}+V_{1}(s,u,t)L(s)\Bigr{]}+\frac{tu}{s^{2}}W_{1}(s,u,t)+\frac{t^{3}}{us^{2}}W_{2}(s,u,t),$
$\displaystyle H_{12}^{\prime(1)}$
$\displaystyle=\frac{t^{2}+u^{2}}{2s^{2}}\Bigl{[}-2(C_{A}+C_{F})(L(s))^{2}+V_{1}(s,t,u)L(s)+V_{1}(s,u,t)L(s)\Bigr{]}$
(54)
$\displaystyle+\frac{u^{2}}{2s^{2}}\Bigl{(}W_{1}(s,u,t)+W_{2}(s,t,u)\Bigr{)}+\frac{t^{2}}{2s^{2}}\Bigl{(}W_{1}(s,t,u)+W_{2}(s,u,t)\Bigr{)},$
$\displaystyle H_{13}^{\prime(1)}$
$\displaystyle=\frac{t}{2s}V_{2}(s,t,u)L(s)+\frac{u^{2}}{2st}L(s)+\frac{t}{2s}W_{4}(s,t,u)+\frac{u^{2}}{2st}W_{4}(s,u,t),$
$\displaystyle H_{23}^{\prime(1)}$
$\displaystyle=\frac{t^{2}}{2su}V_{2}(s,t,u)L(s)+\frac{u}{2s}V_{2}(s,u,t)L(s)+\frac{t^{2}}{2su}W_{4}(s,t,u)+\frac{u}{2s}W_{4}(s,u,t),\
\ H_{33}^{\prime(1)}=0,$

and $H^{\prime}_{JI}=H^{\prime*}_{IJ}$. The functions $V_{i}$ and $W_{i}$ are
given by

$\displaystyle W_{1}(s,t,u)$
$\displaystyle=(C_{A}-C_{F})\frac{s}{u}\Bigl{[}\Bigl{(}L(s)-L(t)\Bigr{)}^{2}+\pi^{2}\Bigr{]}+C_{A}-8C_{F}+(7C_{A}+C_{F})\frac{\pi^{2}}{6},$
(55) $\displaystyle W_{2}(s,t,u)$
$\displaystyle=\Bigl{(}-C_{F}\frac{s^{3}}{u^{3}}-C_{A}\frac{t^{3}+u^{3}-s^{3}}{2u^{3}}\Bigr{)}\Bigl{[}\Bigl{(}L(s)-L(t)\Bigr{)}^{2}+\pi^{2}\Bigr{]}$
$\displaystyle+\Bigl{(}2C_{A}\frac{ts}{u^{2}}+C_{F}\frac{s(2s-u)}{u^{2}}\Bigr{)}\Bigl{(}L(t)-L(s)\Bigr{)}+C_{F}\frac{t-7u}{u}-C_{A}\frac{t}{u}+(7C_{A}+C_{F})\frac{\pi^{2}}{6},$
$\displaystyle W_{3}(s,t,u)$
$\displaystyle=2C_{F}-2C_{A}-\frac{2t}{3s}(C_{A}-n_{f}),$ $\displaystyle
W_{4}(s,t,u)$
$\displaystyle=-\frac{3u}{4t}\Bigl{(}L(s)-L(u)\Bigr{)}^{2}-\Bigl{(}L(s)-L(t)\Bigr{)}\Bigl{(}L(s)-L(u)\Bigr{)}+\frac{3\pi^{2}}{2}\frac{u^{2}}{ts},$
$\displaystyle V_{1}(s,t,u)$
$\displaystyle=3C_{F}-2C_{A}\Bigl{(}L(t)-L(s)\Bigr{)}+\beta_{0}\ \
V_{2}(s,t,u)=\Bigl{(}L(s)-L(u)\Bigr{)}+\frac{t}{s}\Bigl{(}L(t)-L(u)\Bigr{)}.$

The function $L(x)$ as a function of the Mandelstam variables is given by

$L(t)=\ln\frac{-t}{\mu^{2}},\ L(u)=\ln\frac{-u}{\mu^{2}},\
L(s)=\ln\frac{s}{\mu^{2}}-i\pi.$ (56)

In order to obtain $H_{IJ}^{(1)}$ in our basis, we perform the transformation
$\mathbf{H}=\mathbf{A}\mathbf{H}^{\prime}\mathbf{A}^{\dagger}$.

##  7 Soft function

The soft function from $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ is
written in terms of $3\times 3$ matrices in the basis of the operators in eq.
(8). The computation for the virtual and the real contributions is the same as
that in
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, but
the complexity comes from the group theory factors in the soft Wilson lines.
[See eq. (3).] The number of the color factors to be computed is of the order
of 6000 to NLO for $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, while it is
of order 1000 for
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$. We
have constructed a Mathematica package to compute all the group theory
factors.

The relevant Feynman diagrams are shown in fig. 1. The virtual contributions
in fig. 1(a) are obtained by contracting different soft Wilson lines on the
same side of the unitarity cut. The contraction between the same Wilson lines
vanishes because $n_{i}^{2}=0$. The real contribution, shown in fig. 1(b) is
obtained by contracting the soft Wilson lines across the unitary cut.

Figure 1: Feynman diagrams for the emission of a soft gauge boson from the
soft Wilson lines $Y$ or $\mathcal{Y}$. The vertical lines are the final-state
cut. (a) The virtual contribution, (b) the real contribution. The gauge bosons
attached to the same $i$ do not contribute because $n_{i}^{2}=0$. All the
possible pairs of contractions with different $i$ and $j$ should be summed.

As in $\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$,
we decompose the soft function into the hemisphere and the non-hemisphere
parts. We implement the hemisphere function by the constraint function $F$,
with four independent labels, as Chay:2021arz

$F(k,\\{q_{i}\\})=F_{ij,\mathrm{hemi}}(k,\\{q_{i}\\})+F_{ij,ml}(k,\\{q_{i}\\})+F_{ij,lm}(k,\\{q_{i}\\})+(i\leftrightarrow
j),$ (57)

where the hemisphere measurement function for the full hemisphere
$q_{j}>q_{i}$ is given by

$F_{ij,\mathrm{hemi}}(k,\\{q_{i}\\})=\theta(q_{j}-q_{i})\delta(k-q_{i}).$ (58)

The indices $l$, $m$ refer to the remaining lightcone directions when we
contract the Wilson lines in the $i$ and $j$ directions. The non-hemisphere
functions are given as

$\displaystyle F_{ij,ml}(k,\\{q_{i}\\})$
$\displaystyle=\theta(q_{j}-q_{i})\theta(q_{i}-q_{m})\theta(q_{l}-q_{m})\Bigl{(}\delta(k-q_{m})-\delta(k-q_{i})\Bigr{)},$
$\displaystyle F_{ij,lm}(k,\\{q_{i}\\})$
$\displaystyle=\theta(q_{j}-q_{i})\theta(q_{i}-q_{l})\theta(q_{m}-q_{l})\Bigl{(}\delta(k-q_{l})-\delta(k-q_{i})\Bigr{)},$
(59)

which are the non-hemisphere measurement function for the regions $m$ and $l$
respectively. Note that the constraint function $F$ is constructed for the
gauge boson emitted from the soft Wilson lines $Y_{i}$ ($\mathcal{Y}_{i}$) and
$Y_{j}^{\dagger}$ ($\mathcal{Y}_{j}$). The hemisphere function for the $i$ and
$j$ jet directions contains the collinear and the soft divergences. It is
enough to focus on the hemisphere soft function, in which all the divergences
reside. We can obtain the relevant anomalous dimensions at NLL.

We cite the hemisphere virtual contributions $S_{ij,\mathrm{hemi}}^{V}$ and
the real contributions $S_{ij,\mathrm{hemi}}^{V}$ from ref. Chay:2021arz  as

$\displaystyle S_{ij,\mathrm{hemi}}^{V}$
$\displaystyle=\frac{\alpha}{2\pi}\delta(k)\Bigl{[}\frac{1}{\epsilon^{2}}-\frac{2}{\eta}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}-\frac{1}{\epsilon}\ln\frac{n_{ij}\nu^{2}}{\mu^{2}}+\frac{1}{2}\ln^{2}\frac{\mu^{2}}{M^{2}}-\ln\frac{n_{ij}\nu^{2}}{M^{2}}\ln\frac{\mu^{2}}{M^{2}}-\frac{\pi^{2}}{12}\Bigr{]},$
$\displaystyle S_{ij,\mathrm{hemi}}^{R}$
$\displaystyle=-\frac{\alpha}{2\pi}\Bigl{\\{}\delta(k)\Bigl{[}\frac{2}{\eta}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}-\frac{2}{\epsilon^{2}}+\frac{1}{\epsilon}\ln\frac{\nu^{2}}{\mu^{2}}+\ln\frac{\nu^{2}}{M^{2}}\ln\frac{\mu^{2}}{M^{2}}+\frac{\pi^{2}}{6}\Bigr{]}$
(60)
$\displaystyle+\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}\frac{2}{\mu}\mathcal{L}_{0}\Bigl{(}\frac{k}{\mu}\Bigr{)}\theta(k<\sqrt{n_{ij}}M)+\frac{2}{k}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{n_{ij}\mu^{2}}{k^{2}}\Bigr{)}\theta(k>\sqrt{n_{ij}}M)\Bigr{\\}},$

where $n_{ij}=n_{i}\cdot n_{j}/2$. In terms of these contributions, the
hemisphere soft function can be written as

$\mathbf{S}_{\mathrm{hemi}}(a_{1},a_{2},a_{3},a_{4})=\sum_{i\neq
j}\Bigl{[}\mathbf{S}_{ij}^{V}(a_{1},a_{2},a_{3},a_{4})S_{ij,\mathrm{hemi}}^{V}+\mathbf{S}_{ij}^{R}(a_{1},a_{2},a_{3},a_{4})S_{ij,\mathrm{hemi}}^{R}\Bigr{]}.$
(61)

The factors $\mathbf{S}_{ij}^{V,R}(a_{1},a_{2},a_{3},a_{4})$ in front of the
virtual and real contributions represent the corresponding color factors, and
they are represented in terms of the $3\times 3$ matrices for
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$. The indices $a_{i}$ denote the
color index in the presence of the nonsinglets from the originating $i$-th
collinear particle ($i=1,2$ for the incoming particles, and $i=3,4$ for the
outgoing particles in our convention). For example, the soft color matrix with
all the singlets is given by $\mathbf{S}(0,0,0,0)$ and the soft color matrix
with the nonsinglet contributions from 1 and 3 is denoted as
$\mathbf{S}(1,0,1,0)$, etc.. All the soft matrices at tree level are presented
in appendix B.

##  8 Anomalous dimensions and RG evolution

The factorized hard, collinear and soft parts contain the logarithms which
become small at their own characteristic scales. The typical hard scale is
$\mu_{H}\sim Q$, and the collinear scale is given by
$\mu_{C}\sim\sqrt{Q\mathcal{T}}\sim\sqrt{QM}$, while the soft scale is
$\mu_{S}\sim M$. However, if the common factorization scale $\mu_{F}$ is away
from these characteristic scales, the logarithms become so large that the
perturbation theory breaks down. We can resum large logarithms by evolving the
factorized parts from their characteristic scales to the common factorization
scale $\mu_{F}$ by solving the RG equations.

The nonsinglet contributions contain additional logarithms associated with the
rapidity divergence. The characteristic collinear and soft rapidity scales are
$\nu_{C}\sim Q$, $\nu_{S}\sim M$ respectively. We perform the double evolution
with respect to the rapidity scale $\nu$, as well as the renormalization scale
$\mu$. Because the order of the evolution is irrelevant Chiu:2012ir , we first
evolve with respect to $\nu$ from $\nu_{C}\sim Q$ for the collinear functions
and from $\nu_{S}\sim M$ for the soft functions to the common factorization
scale $\nu_{F}$. Then we evolve all the factorized functions with respect to
$\mu$. The path of the evolution is shown in fig. 2.

Figure 2: The evolutions of the hard, collinear and soft functions start from
their own characteristic scales $\mu_{H}$, ($\mu_{C}$, $\nu_{C}$), and
($\mu_{S}$, $\nu_{S}$) respectively. Since the order of the evolution is
irrelevant, we first evolve the collinear and soft functions with respect to
$\nu$, and then $\mu$. The path of the evolution is illustrated by the arrows.

###  8.1 Collinear functions

There are various collinear functions involved in the 2-jettiness. They are
the beam functions for the leptons and the gauge bosons, the corresponding
PDFs , and the semi-inclusive muon jet functions. Because of the crossing
symmetry, the anomalous dimensions of the beam functions and the corresponding
jet functions turn out to be the same.

It is convenient to express the $N$-jettiness in terms of the Laplace
transforms, in which the factorization consists of the products instead of the
convolution. The actual $N$-jettiness is obtained by taking the inverse
Laplace transform, which will be discussed more in the numerical estimates.
The Laplace transforms for the beam functions are written as

$\displaystyle\tilde{B}_{i}\Bigl{(}\ln\frac{\omega
Q_{L}}{\mu^{2}},z,M,\mu\Bigr{)}$
$\displaystyle=\int_{0}^{\infty}dk\,e^{-sk}B_{i}(\omega k,z,M,\mu),$
$\displaystyle B_{i}(\omega k,z,M,\mu)$ $\displaystyle=\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}ds\,e^{sk}\tilde{B}_{i}\Bigl{(}\ln\frac{1}{e^{\gamma_{\mathrm{E}}}s\mu},M,\mu\Bigr{)},$
(62)

with $s=1/(e^{\gamma_{\mathrm{E}}}Q_{L})$. Here $Q_{L}$ is the scale
introduced in performing the Laplace transform. The Laplace transforms of the
PDFs and the jet functions are defined similarly.

Let us denote $\tilde{f}_{s}$ and $\tilde{f}_{n}$ as the singlet and the
nonsinglet Laplace-transformed collinear functions, which can be either the
beam functions, or the semi-inclusive muon jet functions. The singlet
functions evolve with respect to the renormalization scale $\mu$ only, and the
RG equation is written as

$\frac{d}{d\ln\mu}\tilde{f}_{s}(\mu)=\tilde{\gamma}_{fs}\tilde{f}_{s}(\mu),$
(63)

and the nonsinglet functions $\tilde{f}_{n}$ satisfy the coupled RG equations,
which are given as

$\frac{d}{d\ln\mu}\tilde{f}_{n}(\mu,\nu)=\tilde{\gamma}_{fn}^{\mu}\tilde{f}_{n}(\mu,\nu),\
\
\frac{d}{d\ln\nu}\tilde{f}_{n}(\mu,\nu)=\tilde{\gamma}_{fn}^{\nu}\tilde{f}_{n}(\mu,\nu).$
(64)

Here $\tilde{\gamma}_{fs}$, $\tilde{\gamma}_{fn}^{\mu}$ and
$\tilde{\gamma}_{fn}^{\nu}$ are the anomalous dimensions of the singlets and
the nonsinglets. The $\mu$ and $\nu$ anomalous dimensions for the beam
functions of the gauge bosons in the Laplace transform are given as

$\displaystyle\tilde{\gamma}^{\mu}_{Ws}=2C_{A}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\omega}-2\gamma_{W},\
\ \tilde{\gamma}^{\nu}_{Ws}=0,$
$\displaystyle\tilde{\gamma}^{\mu}_{Wn}=\tilde{\gamma}^{\mu}_{Ws}-C_{A}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\nu},\
\ \tilde{\gamma}^{\nu}_{Wn}=C_{A}\Gamma_{c}\ln\frac{\mu}{M},$ (65)

where the non-cusp anomalous dimension at NLO is given as
$\gamma_{W}^{(1)}=-\beta_{0}\alpha/(4\pi)$. The cusp anomalous dimension
Korchemsky:1987wg ; Korchemskaya:1992je  $\Gamma_{c}(\alpha)$ can be expanded
as

$\Gamma_{c}(\alpha)=\frac{\alpha}{4\pi}\Gamma_{c}^{0}+\Bigl{(}\frac{\alpha}{4\pi}\Bigr{)}^{2}\Gamma_{c}^{1}+\cdots,$
(66)

with

$\Gamma_{c}^{0}=4,\
\Gamma_{c}^{1}=\Bigl{(}\frac{268}{9}-\frac{4}{3}\pi^{2}\Bigr{)}C_{A}-\frac{40n_{f}}{9}.$
(67)

To NLL accuracy, the cusp anomalous dimension to two loops is needed.

The anomalous dimensions for the beam function of the gauge bosons in eq.
(8.1), the beam function of the leptons and the semi-inclusive jet functions
Chay:2021arz  can be written in a general form as

$\displaystyle\tilde{\gamma}_{is}^{\mu}=2C_{i}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\omega}-2\gamma_{i},\
\
\tilde{\gamma}_{in}^{\mu}=\tilde{\gamma}_{is}^{\mu}-C_{A}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\nu},$
$\displaystyle\tilde{\gamma}_{is}^{\nu}=0,\ \
\tilde{\gamma}_{in}^{\nu}=C_{A}\Gamma_{c}\ln\frac{\mu}{M},$ (68)

with $i=\ell,W$, and the group theory factors are given as $C_{\ell}=C_{F}$,
$C_{W}=C_{A}$. The non-cusp anomalous dimensions at order $\alpha$ are given
by

$\gamma_{\ell}^{(1)}=-3C_{F}\frac{\alpha}{4\pi},\ \
\gamma_{W}^{(1)}=-\beta_{0}\frac{\alpha}{4\pi}.$ (69)

The evolution of the singlet collinear functions from eq. (63) can be written
as

$\tilde{f}_{is}(\mu_{F})=U_{is}(\mu_{F},\mu_{C})\tilde{f}_{is}(\mu_{C}),$ (70)

with $i=\ell,W$, where the evolution kernel is given as

$U_{is}(\mu_{F},\mu_{C})=\exp\Bigl{[}4C_{i}S_{\Gamma}(\mu_{F},\mu_{C})-2C_{i}a_{\Gamma}(\mu_{F},\mu_{C})\ln\frac{\omega
Q_{L}}{\mu_{C}^{2}}-2a_{\gamma_{i}}(\mu_{F},\mu_{C})\Bigr{]}.$ (71)

The various quantities in eq. (71) are defined as

$\int_{\mu_{C}}^{\mu_{F}}\frac{d\mu}{\mu}\Gamma_{c}(\alpha)\ln\frac{\mu^{2}}{\omega
Q_{L}}\equiv
2S_{\Gamma}(\mu_{F},\mu_{C})-a_{\Gamma}(\mu_{F},\mu_{C})\ln\frac{\omega
Q_{L}}{\mu_{C}^{2}}.$ (72)

where

$S_{\Gamma}(\mu,\mu_{i})=\int_{\alpha(\mu_{i})}^{\alpha(\mu)}d\alpha\frac{\Gamma_{c}(\alpha)}{\beta(\alpha)}\int_{\alpha(\mu_{i})}^{\alpha}\frac{d\alpha^{\prime}}{\beta(\alpha^{\prime})},\
\
a_{f}(\mu,\mu_{i})=\int_{\alpha(\mu_{i})}^{\alpha(\mu)}\frac{d\alpha}{\beta(\alpha)}f(\alpha),$
(73)

with $f(\alpha)=\Gamma_{c}(\alpha)$ or any other noncusp anomalous dimensions
such as $\gamma_{\ell}$ or $\gamma_{W}$.

On the other hand, the nonsinglet collinear functions satisfy the coupled RG
equations, eq. (64). And the evolution from both RG equations are written as

$\tilde{f}_{in}(\mu_{F},\nu_{F})=U_{in}(\mu_{F},\mu_{C};\nu_{F})V_{in}(\nu_{F},\nu_{C};\mu_{C})\tilde{f}_{in}(\mu_{C},\nu_{C}),$
(74)

where the evolution kernels are given as

$\displaystyle U_{in}(\mu_{F},\mu_{C};\nu_{F})$
$\displaystyle=U_{is}(\mu_{F},\mu_{C})\exp[-2C_{A}S_{\Gamma}(\mu_{F},\mu_{C})+C_{A}a_{\Gamma}(\mu_{F},\mu_{C})\ln\frac{\nu_{F}Q_{L}}{\mu_{C}^{2}}\Bigr{]},$
$\displaystyle V_{in}(\nu_{F},\nu_{C};\mu_{C})$
$\displaystyle=\exp\Bigl{[}C_{A}a_{\Gamma}(\mu_{C},M)\ln\frac{\nu_{F}}{\nu_{C}}\Bigr{]}.$
(75)

###  8.2 Hard function

The RG equation for the hard function $\mathbf{H}_{W}$ in the channel
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ is written as

$\frac{d}{d\ln\mu}\mathbf{H}_{W}=\mathbf{\Gamma}_{H}^{W}\mathbf{H}_{W}+\mathbf{H}_{W}\mathbf{\Gamma}_{H}^{W\dagger}.$
(76)

The relation between the anomalous dimensions $\bm{\Gamma}_{H}^{W}$ in our
basis, and the anomalous dimensions $\bm{\Gamma}_{H}^{\prime W}$ in ref.
Kelley:2010fn  is given by $\bm{\Gamma}_{H}^{W}=A\bm{\Gamma}_{H}^{\prime
W}A^{-1}$. Note that only the nondiagonal matrix is affected by this
transformation.

The anomalous dimension matrix Kelley:2010fn ; Chiu:2009mg  can be generalized
to include the channels $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ and
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$. The
anomalous dimensions $\bm{\Gamma}_{i}$ ($i=W,e$) of the hard functions
$\mathbf{H}_{i}$ can be written as

$\mathbf{\Gamma}_{H}^{i}(s,t,u)=\Bigl{[}\frac{c_{H}}{2}\Gamma_{c}(\alpha)L_{i}(\mu)+\gamma_{H}\Bigr{]}\mathbf{1}+\Gamma_{c}(\alpha)\mathbf{M}_{i}(s,t,u),$
(77)

where $c_{H}=n_{\ell}C_{F}+n_{W}C_{A}$ with $n_{\ell}$ fermions and $n_{W}$
gauge bosons involved in the channel. And
$\gamma_{H}=n_{\ell}\gamma_{\ell}+n_{W}\gamma_{W}$, with $L_{W}(\mu)=L(t)$,
$L_{e}(\mu)=L(s)$. The matrices $\mathbf{M}_{W}$ and $\mathbf{M}_{e}$, in the
two channels, are given by

$\displaystyle\mathbf{M}_{W}$
$\displaystyle=\begin{pmatrix}\displaystyle\frac{C_{A}+C_{F}}{2}\ln\frac{n_{12}n_{34}}{n_{13}n_{24}}&\displaystyle\frac{1}{2}\ln\frac{n_{13}n_{24}}{n_{14}n_{23}}&0\\\
\displaystyle\ln\frac{n_{13}n_{24}}{n_{14}n_{23}}&\displaystyle\frac{C_{A}}{4}\ln\frac{n_{14}n_{23}}{n_{13}n_{24}}+\frac{C_{F}}{2}\ln\frac{n_{12}n_{34}}{n_{13}n_{24}}&\displaystyle\frac{-3C_{A}+8C_{F}}{4}\ln\frac{n_{13}n_{24}}{n_{14}n_{23}}\\\
0&\displaystyle\frac{C_{A}}{4}\ln\frac{n_{13}n_{24}}{n_{14}n_{23}}&\displaystyle\frac{C_{A}}{4}\ln\frac{n_{14}n_{23}}{n_{13}n_{24}}+\frac{C_{F}}{2}\ln\frac{n_{12}n_{34}}{n_{13}n_{24}}\end{pmatrix},$
$\displaystyle\mathbf{M}_{e}$
$\displaystyle=\begin{pmatrix}\displaystyle\Bigl{(}2C_{F}-C_{A}/2\Bigr{)}\ln\frac{n_{13}n_{24}}{n_{14}n_{23}}-\frac{C_{A}}{2}\ln\frac{n_{12}n_{34}}{n_{14}n_{23}}&&\displaystyle\ln\frac{n_{13}n_{24}}{n_{14}n_{23}}\\\
\displaystyle-
C_{F}\Bigl{(}C_{F}-C_{A}/2\Bigr{)}\ln\frac{n_{13}n_{24}}{n_{14}n_{23}}&&0\end{pmatrix}.$
(78)

The evolution of the hard functions from the hard scale $\mu_{H}$ to the
factorization scale $\mu_{F}$ is written in the form

$\mathbf{H}_{i}(\mu_{F})=\Pi_{H}^{i}(\mu_{F},\mu_{H})\bm{\Pi}_{H}^{i}(\mu_{F},\mu_{H})\mathbf{H}_{i}(\mu_{H})\bm{\Pi}_{H}^{i\dagger}(\mu_{F},\mu_{H}),$
(79)

where $\Pi_{H}^{i}(\mu_{F},\mu_{H})$ describes the evolution from the
anomalous dimension, proportional to the idendity matrix, while
$\bm{\Pi}_{H}^{i}(\mu_{F},\mu_{H})$ is the evolution kernel from
$\mathbf{M}_{i}$. They are given as

$\displaystyle\Pi_{H}^{i}(\mu_{F},\mu_{H})$
$\displaystyle=\exp\Bigl{[}-2c_{H}S(\mu_{F},\mu_{H})+c_{H}a_{\Gamma}(\mu_{F},\mu_{H})L_{i}(\mu_{H})+2a_{\gamma_{H}}(\mu_{F},\mu_{H})\Bigr{]},$
$\displaystyle\bm{\Pi}_{H}^{i}(\mu_{F},\mu_{H})$
$\displaystyle=\exp\Bigl{[}a_{\Gamma}(\mu_{F},\mu_{H})\mathbf{M}_{i}\Bigr{]}.$
(80)

###  8.3 Soft function

The RG equation for the soft function $\tilde{\mathbf{S}}$ is written as

$\frac{d}{d\ln\mu}\tilde{\mathbf{S}}=\tilde{\bm{\Gamma}}_{S}^{\mu\dagger}\tilde{\mathbf{S}}+\tilde{\mathbf{S}}\tilde{\bm{\Gamma}}_{S}^{\mu},$
(81)

where $\tilde{\bm{\Gamma}}_{S}^{\mu}$ is the $\mu$-anomalous dimension matrix.
The relation of the soft functions between our basis $T_{I}$ and another basis
$T^{\prime}_{I}$ is obtained by requiring that
$\mathrm{tr}(\mathbf{H}^{\prime}\mathbf{S}^{\prime})=\mathrm{tr}(\mathbf{H}\mathbf{S})$,
from which we obtain that
$\mathbf{S}=(\mathbf{A}^{-1})^{\dagger}\mathbf{S}^{\prime}\mathbf{A}^{-1}$,
and
$\tilde{\bm{\Gamma}}_{S}^{\mu}=\mathbf{A}\tilde{\bm{\Gamma}}_{S}^{\prime\mu}\mathbf{A}^{-1}$.

At order $\alpha$, eq. (81) is written as

$\frac{d}{d\ln\mu}\tilde{\mathbf{S}}^{(1)}=\tilde{\bm{\Gamma}}_{S}^{\mu\dagger(1)}\tilde{\mathbf{S}}^{(0)}+\tilde{\mathbf{S}}^{(0)}\tilde{\bm{\Gamma}}_{S}^{\mu(1)},$
(82)

where $\tilde{\mathbf{S}}^{(1)}$ is the renormalized soft function. The
$\mu$-soft anomalous dimensions can be extracted from the requirement that the
sum of the anomalous dimensions from all the factorized parts should cancel.
The $\mu$-soft anomalous dimensions at NLO are written as

$\tilde{\bm{\Gamma}}_{iS}^{\mu(1)}=-\Bigl{(}\tilde{\bm{\Gamma}}_{iH}^{(1)}+(\tilde{\gamma}_{i}^{\mu(1)}+\tilde{\gamma}_{J}^{\mu(1)})\otimes\mathbf{1}\Bigr{)},$
(83)

where $\tilde{\gamma}_{i}^{\mu}$ are the anomalous dimensions of the beam
functions ($i=\ell,W$), and $\tilde{\gamma}_{J}^{\mu}$ is that of the jet
function. And the RG equation with respect to the rapidity scale $\nu$ is
written as

$\frac{d}{d\ln\nu}\tilde{\mathbf{S}}=\tilde{\Gamma}_{S}^{\nu}\tilde{\mathbf{S}}.$
(84)

The soft anomalous dimensions for
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ are given by

$\displaystyle\Bigl{(}\Gamma_{S}^{k\mu}\Bigr{)}_{W}$
$\displaystyle=\Gamma_{c}\Bigl{[}\frac{c_{H}}{4}\ln\frac{Q_{L}^{4}}{n_{13}n_{24}\mu^{4}}+\frac{1}{2}(C_{F}-C_{A})\ln\frac{n_{12}}{n_{34}}\Bigr{]}-\Gamma_{c}\mathbf{M}_{W}+\frac{kC_{A}}{2}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\nu},$
$\displaystyle\Bigl{(}\Gamma_{S}^{k\nu}\Bigr{)}_{W}$
$\displaystyle=-kC_{A}\Gamma_{c}\ln\frac{\mu}{M},$ (85)

where $k$ ($k=0,2,3,4$) is the number of the nonsinglets involved in the
process.333The case with $k=1$ vanishes due to the color conservation, and
only the even number of nonsinglets is allowed in SU(2). We can write the soft
anomalous dimensions in a compact form including both channels
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ and
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, respectively. The Laplace
transformed $\mu$ and $\nu$ soft anomalous dimensions can be written as

$\displaystyle\Bigl{(}\tilde{\Gamma}_{Sk}^{\mu}\Bigr{)}_{i}$
$\displaystyle=\Gamma_{c}\Bigl{[}\frac{c_{H}}{4}\ln\frac{Q_{L}^{4}}{\mu^{4}N_{i}}+B_{i}\Bigr{]}\mathbf{1}-\Gamma_{c}\mathbf{M}_{i}+k\frac{C_{A}}{2}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\nu}\mathbf{1},$
$\displaystyle\Bigl{(}\tilde{\Gamma}_{Sk}^{\nu}\Bigr{)}_{i}$
$\displaystyle=-kC_{A}\Gamma_{c}\ln\frac{\mu}{M},$ (86)

where $i=\ell$ or $W$. The quantities in eq. (8.3) are given as

$N_{\ell}=n_{12}n_{34},\ \ N_{W}=n_{13}n_{24},\ \ B_{\ell}=0,\ \
B_{W}=\frac{1}{2}(C_{F}-C_{A})\ln\frac{n_{12}}{n_{34}}.$ (87)

Note that $B_{W}=0$ for $n_{12}=n_{34}$ when the jets are back-to-back. The
list of the anomalous dimensions for all the factorized parts is shown in
table 1.

Functions | Type |  $\mu$-anomalous dimensions |  $\nu$-anomalous dimensions
---|---|---|---
Hard function | $-$ | $\bm{\Gamma}_{H}^{i}=\Bigl{[}\frac{c_{H}}{2}\Gamma_{c}L_{i}+\gamma_{H}\Bigr{]}\mathbf{1}+\Gamma_{c}\mathbf{M}_{i}$ | 0
Lepton Beam | singlet | $\tilde{\gamma}_{\ell s}^{\mu}=2C_{F}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\omega}-2\gamma_{\ell}$ | 0
nonsinglet | $\tilde{\gamma}_{\ell n}^{\mu}=\tilde{\gamma}_{\ell s}^{\mu}-C_{A}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\nu}$ | $\tilde{\gamma}_{\ell n}^{\nu}=C_{A}\Gamma_{c}\ln\frac{\mu}{M}$
Boson Beam | singlet | $\tilde{\gamma}_{Ws}^{\mu}=2C_{A}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\omega}-2\gamma_{W}$ | 0
nonsinglet | $\tilde{\gamma}_{Wn}^{\mu}=\tilde{\gamma}_{Ws}^{\mu}-C_{A}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\nu}$ | $\tilde{\gamma}_{bn}^{\nu}=C_{A}\Gamma_{c}\ln\frac{\mu}{M}$
Muon jet | singlet | $\tilde{\gamma}_{\ell s}^{\mu}=2C_{F}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\omega}-2\gamma_{\ell}$ | 0
nonsinglet | $\tilde{\gamma}_{\ell n}^{\mu}=\tilde{\gamma}_{\ell s}^{\mu}-C_{A}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\nu}$ | $\tilde{\gamma}_{\ell n}^{\nu}=C_{A}\Gamma_{c}\ln\frac{\mu}{M}$
Soft function | singlet | $(\tilde{\bm{\Gamma}}_{S}^{\mu})_{i}=\Gamma_{c}\Bigl{[}\frac{c_{H}}{4}\ln\frac{Q_{L}^{4}}{\mu^{4}N_{i}}+B_{i}\Bigr{]}\mathbf{1}-\Gamma_{c}\mathbf{M}_{i}$ | 0
nonsinglet | $(\tilde{\bm{\Gamma}}_{Sk}^{\mu})_{i}=(\bm{\Gamma}_{S}^{\mu})_{i}+k\frac{C_{A}}{2}\Gamma_{c}\ln\frac{\mu^{2}}{Q_{L}\nu}\mathbf{1}$ | $(\tilde{\Gamma}_{Sk}^{\nu})_{i}=-kC_{A}\Gamma_{c}\ln\frac{\mu}{M}$
Table 1: Anomalous dimensions of the factorized parts.

The evolution of the Laplace-transformed soft functions is written as

$\tilde{S}_{k}^{i}(\mu_{F},\nu_{F})=\Pi_{Sk}^{i}(\mu_{F},\mu_{S};\nu_{F},\nu_{S})\bm{\Pi}^{i\dagger}(\mu_{F},\mu_{S})\tilde{S}_{k}^{i}(\mu_{S},\nu_{S})\bm{\Pi}_{S}^{i}(\mu_{F},\mu_{S}),$
(88)

where $\bm{\Pi}^{i}$ is the evolution kernel from $\mathbf{M}_{i}$
($i=\ell,W$). And $\Pi_{Sk}^{i}$ is the kernel from the diagonal part, with
$k$ nonsinglets. $\Pi_{Sk}^{i}$ involves a double evolution with respect to
$\mu$ and $\nu$. Since the order of the evolution is irrelevant, we first
evolve with respect to $\nu$, and then $\mu$. As a result, $\Pi_{Sk}^{i}$ can
be written as

$\Pi_{Sk}^{i}(\mu_{F},\mu_{S},\nu_{F},\nu_{S})=U_{Sk}^{i}(\mu_{F},\mu_{S};\nu_{F})V_{Sk}^{i}(\nu_{F},\nu_{S};\mu_{S}),$
(89)

where

$\displaystyle V_{Sk}^{i}(\nu_{F},\nu_{S};\mu_{S})$
$\displaystyle=\exp\Bigl{[}-kC_{A}a_{\Gamma}(\mu_{S},M)\ln\frac{\nu_{F}}{\nu_{S}}\Bigr{]},$
(90) $\displaystyle U_{Sk}^{i}(\mu_{F},\mu_{S};\nu_{F})$
$\displaystyle=\exp\Bigl{[}-2c_{H}S_{\Gamma}(\mu_{F},\mu_{S})+\frac{c_{H}}{2}a_{\Gamma}(\mu_{F},\mu_{S})\ln\frac{Q_{L}^{4}}{\mu_{S}^{4}N_{i}}$
$\displaystyle+B_{i}a_{\Gamma}(\mu_{F},\mu_{S})+2kC_{A}S_{\Gamma}(\mu_{F},\mu_{S})+kC_{A}a_{\Gamma}(\mu_{F},\mu_{S})\ln\frac{\mu_{S}^{2}}{Q_{L}\nu_{F}}\Bigr{]}.$

##  9 Numerical estimates of the 2-jettiness in SU(2) near threshold

In order to estimate the singlet and the nonsinglet contributions to the
2-jettiness, we confine ourselves to the threshold limit in which the two jets
are emitted back-to-back near $\theta=\pi/2$ (perpendicular to the beam axis)
at NLL. It not only simplifies the kinematics, but also is available from the
current information obtained so far. We explain these points first, and
present the numerical analysis accordingly.

Firstly, the anomalous dimensions depend on the angles among the light-like
directions. [See, for example, eqs. (8.2) and (8.3).] Therefore the evolution
of the 2-jettiness in eq (3) involves the integration over these angles. The
numerical integration over the angles can be performed, but we would rather
consider the differential 2-jettiness $d\sigma/d\mathcal{T}_{2}dt$ at
$\theta=\pi/2$. It corresponds to the back-to-back jets, perpendicular to the
beam directions.

Secondly, the information on the PDFs is not available, in contrast to QCD.
That is, the information of the probability in finding a lepton with a certain
fraction of the longitudinal momentum is lacking. We naively expect that the
PDFs satisfy $f_{i}(x)\propto\delta(1-x)$ for $i=\ell,W$ at tree level, but
the detailed functional forms at higher orders are not known. These delta
functions will be adopted in our analysis. If we consider the normalized
2-jettiness
$(d\sigma_{i}/d\mathcal{T}_{2}dt)/(d\sigma_{i}^{(0)}/d\mathcal{T}_{2}dt)$,
where $d\sigma_{i}^{(0)}/d\mathcal{T}_{2}dt$ is the Born result, the
uncertainties on the PDF or the beam functions can be somewhat reduced in the
ratio.

Finally, we have computed the hemisphere soft functions only at NLO, from
which the anomalous dimensions can be obtained. The non-hemisphere part of the
soft function is finite, and is needed for the computation beyond NLL. Here we
confine to NLL accuracy such that only the hemisphere soft function is needed.

In summary, we consider the 2-jettiness near threshold and at NLL accuracy,
and the jets are produced perpendicular to the beam directions. In this case,
$\omega_{1}=\omega_{2}=\omega_{3}=\omega_{4}=\omega$,
$t=(p_{1}-p_{3})^{2}=-\omega_{1}\omega_{3}n_{1}\cdot
n_{3}/2=-\omega^{2}(1-\cos\theta)/2$ by choosing $n_{1}^{\mu}=(1,0,0,1)$,
$n_{2}^{\mu}=(1,0,0,-1)$, $n_{3}^{\mu}=(1,\sin\theta,0,\cos\theta)$ and
$n_{4}^{\mu}=(1,-\sin\theta,0,-\cos\theta)$.

With the SU(2) gauge interaction, the nonsinglet beam functions for the gauge
bosons do not contribute ($d^{abc}=0$), and the beam functions for the leptons
and the semi-inclusive muon jet functions are employed at NLL. The product of
the hard and soft functions for $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$
is written as

$\displaystyle\sum_{IJ}H_{IJ}(\mu_{F})\tilde{S}_{JI}^{00gh}(\mu_{F},\nu_{F})$
$\displaystyle=\exp\Bigl{[}-4(C_{A}+C_{F})S_{\gamma}(\mu_{F},\mu_{H})-2(C_{A}+C_{F})a_{\Gamma}(\mu_{F},\mu_{H})\ln\frac{\mu_{H}^{2}}{\omega^{2}n_{13}}$
$\displaystyle+4\Bigl{(}a_{\gamma_{W}}(\mu_{F},\mu_{H})+a_{\gamma_{\mu}}(\mu_{F},\mu_{H})\Bigr{)}\Bigr{]}$
$\displaystyle\times\exp\Bigl{[}-4(C_{A}+C_{F})S_{\Gamma}(\mu_{F},\mu_{S})+2(C_{A}+C_{F})a_{\Gamma}(\mu_{F},\mu_{S})\ln\frac{Q_{L}^{2}}{\mu_{S}^{2}n_{13}}\Bigr{]}$
$\displaystyle\times\exp\Bigl{[}k\Bigl{(}2C_{A}S_{\Gamma}(\mu_{F},\mu_{S})+C_{A}a_{\Gamma}(\mu_{F},\mu_{S})\ln\frac{\mu_{S}^{2}}{Q_{L}\nu_{F}}-C_{A}a_{\Gamma}(\mu_{S},M)\ln\frac{\nu_{F}}{\nu_{S}}\Bigl{)}\Bigr{]}$
$\displaystyle\times\mathrm{Tr}\Bigl{(}H^{(0)}\exp\Bigl{[}(a_{\Gamma}(\mu_{S},\mu_{H})\mathbf{M}^{\dagger}\Bigr{]}\mathbf{S}^{(0)k}\exp\Bigl{[}a_{\Gamma}(\mu_{S},\mu_{H})\mathbf{M}\Bigr{]}\Bigr{)}.$
(91)

Note that there should be an even number of nonsinglets in the jet functions.

Taking the inverse Laplace transform of the collinear and soft functions by
noting that

$\mathcal{L}^{-1}[Q_{L}^{a}]=\mathcal{L}^{-1}\Bigl{[}\Bigl{(}\frac{1}{s\exp(\gamma_{\mathrm{E}})}\Bigr{)}^{a}\Bigr{]}=\frac{1}{\mathcal{T}_{2}\Gamma(a)}\Bigl{(}\frac{\mathcal{T}_{2}}{\exp(\gamma_{\mathrm{E}})}\Bigr{)}^{a}=\frac{\exp\Bigl{(}a\ln\mathcal{T}_{2}\Bigr{)}}{\mathcal{T}_{2}\Gamma(a)\exp(a\gamma_{\mathrm{E}})},$
(92)

the differential 2-jettiness from $WW\rightarrow\mu^{-}\mu^{+}$ is written as

$\displaystyle\frac{d\sigma_{W}}{d\mathcal{T}_{2}dt}$
$\displaystyle=\frac{4\pi\alpha^{2}}{\omega^{4}}\sum_{k=0,2}\Bigl{(}\frac{1}{6}\Bigr{)}^{2}\Bigl{(}\frac{1}{4}\Bigr{)}^{2-k}\Bigl{(}\frac{1}{2}\Bigr{)}^{k}\frac{1}{\mathcal{T}_{2}\Gamma(a_{k})\exp[a_{k}\gamma_{\mathrm{E}}]}$
$\displaystyle\times\exp\Bigl{[}-4(C_{A}+C_{F})S_{\Gamma}(\mu_{F},\mu_{H})-2(C_{A}+C_{F})a_{\Gamma}(\mu_{F},\mu_{H})\ln\frac{\mu_{H}^{2}}{\omega^{2}n_{13}}$
$\displaystyle+4\Bigl{(}a_{\gamma_{W}}(\mu_{F},\mu_{H})+a_{\gamma_{\ell}}(\mu_{F},\mu_{H})\Bigr{)}\Bigr{]}$
$\displaystyle\times\exp\Bigl{[}8(C_{A}+C_{F})S_{\Gamma}(\mu_{F},\mu_{C})-4(C_{A}+C_{F})a_{\Gamma}(\mu_{F},\mu_{C})\ln\frac{\omega\mathcal{T}_{2}}{\mu_{C}^{2}}$
$\displaystyle-4\Bigl{(}a_{\gamma_{W}}(\mu_{F},\mu_{C})+a_{\gamma_{\ell}}(\mu_{F},\mu_{C})\Bigr{)}$
$\displaystyle+kC_{A}\Bigl{(}-2S_{\Gamma}(\mu_{F},\mu_{C})+a_{\Gamma}(\mu_{F},\mu_{C})\ln\frac{\nu_{F}\mathcal{T}_{2}}{\mu_{C}^{2}}+a_{\Gamma}(\mu_{F},M)\ln\frac{\nu_{F}}{\nu_{C}}\Bigr{)}\Bigr{]}$
$\displaystyle\times\exp\Bigl{[}-4(C_{A}+C_{F})S_{\Gamma}(\mu_{F},\mu_{S})+2(C_{A}+C_{F})a_{\Gamma}(\mu_{F},\mu_{S})\ln\frac{\mathcal{T}_{2}^{2}}{\mu_{S}^{2}n_{13}}$
$\displaystyle+kC_{A}\Bigl{(}2S_{\Gamma}(\mu_{F},\mu_{S})+a_{\Gamma}(\mu_{F},\mu_{S})\ln\frac{\mu_{S}^{2}}{\nu_{F}\mathcal{T}_{2}}-a_{\Gamma}(\mu_{S},M)\ln\frac{\nu_{F}}{\nu_{S}}\Bigr{)}\Bigr{]}$
$\displaystyle\times\mathrm{Tr}\Bigl{\\{}\mathbf{H}_{W}^{(0)}\exp\Bigl{[}a_{\Gamma}(\mu_{S},\mu_{H})\mathbf{M}_{W}^{\dagger}\Bigr{]}\mathbf{S}_{W}^{(0)k}\exp\Bigl{[}a_{\Gamma}(\mu_{S},\mu_{H})\mathbf{M}_{W}^{\dagger}\Bigr{]}\Bigr{\\}},$
(93)

where $a_{k}$ is given by

$a_{k}=[4(C_{A}+C_{F})-kC_{A}]a_{\Gamma}(\mu_{C},\mu_{S}),$ (94)

and $n_{13}=(1-\cos\theta)/2$. And the numerical factors in front come from
the product of the color factors $k_{i}^{a}$ in eq. (3) and the muon
projections in eq. (46).

The matrix $\mathbf{M}_{W}$ near threshold is given as

$\mathbf{M}_{W}^{\mathrm{th}}=\begin{pmatrix}-(C_{F}+C_{A})\ln
n_{13}&\displaystyle\ln\frac{n_{13}}{1-n_{13}}&0\\\ \displaystyle
2\ln\frac{n_{13}}{1-n_{13}}&\displaystyle\frac{C_{A}}{2}\ln\frac{1-n_{13}}{n_{13}}-C_{F}\ln
n_{13}&\displaystyle\frac{1}{2}(-3C_{A}+8C_{F})\ln\frac{n_{13}}{1-n_{13}}\\\
0&\displaystyle\frac{C_{A}}{2}\ln\frac{n_{13}}{1-n_{13}}&\displaystyle\frac{C_{A}}{2}\ln\frac{1-n_{13}}{n_{13}}+C_{F}\ln
n_{13}\end{pmatrix}.$ (95)

The hard functions and the soft functions at LO are given as

$\displaystyle\mathbf{H}_{W}^{(0)}$
$\displaystyle=\frac{1}{n_{13}(1-n_{13})}\begin{pmatrix}\displaystyle\frac{1-2n_{13}(1-n_{13})}{N^{2}}&\displaystyle\frac{(1-n_{13})^{4}-n_{13}^{4}}{N}&\displaystyle\frac{1-2n_{13}(1-n_{13})}{N}\\\
\displaystyle\frac{(1-n_{13})^{4}-n_{13}^{4}}{N}&(1-2n_{13})^{2}\Bigl{(}1-2n_{13}(1-n_{13})\Bigr{)}&(1-n_{13})^{4}-n_{13}^{4}\\\
\displaystyle\frac{1-2n_{13}(1-n_{13})}{N}&(1-n_{13})^{4}-n_{13}^{4}&1-2n_{13}(1-n_{13})\end{pmatrix},$
$\displaystyle\mathbf{S}_{W}^{(0)0}$
$\displaystyle=S^{(0)}(0,0,0,0)=\begin{pmatrix}2C_{F}C_{A}^{2}&0&0\\\
0&C_{F}C_{A}^{2}&0\\\ 0&0&\displaystyle 0\end{pmatrix},$
$\displaystyle\mathbf{S}_{W}^{(0)2}$
$\displaystyle=S^{(0)}(0,0,1,1)=\begin{pmatrix}2C_{F}C_{A}&0&0\\\
0&\displaystyle C_{A}(C_{F}-C_{A}/2)&0\\\ 0&0&0\end{pmatrix}(34).$ (96)

Here $(a_{1}a_{2}\cdots a_{n})$ denotes $\mathrm{Tr}(t^{a_{1}}t^{a_{2}}\cdots
t^{a_{n}})$. All the other soft matrices at treel level are presented in
appendix B.

Note that the hard and soft anomalous dimensions, hence the corresponding
evolution kernels, depend on $n_{13}=(1-\cos\theta)/2$, as well as the hard
functions themselves. And we fix $\theta=\pi/2$ in the numerical estimation.
It corresponds to the back-to-back jets perpendicular to the beam axis.

Similarly, we can write the differential 2-jettiness for the channel
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$. In
contrast to the channel $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, there
are nonsinglet contributions from the electron beam functions. The only
constraint is that there should be even number of nonsinglet contributions as
a whole. It is given as

$\displaystyle\frac{d\sigma_{e}}{d\mathcal{T}_{2}dt}$
$\displaystyle=\frac{64\pi\alpha^{2}}{\omega^{4}}\sum_{k=0,2,4}\Bigl{(}\frac{1}{4}\Bigr{)}^{4-k}\Bigl{(}\frac{1}{2}\Bigr{)}^{k}\frac{1}{\mathcal{T}_{2}\Gamma(a_{k})\exp(a_{k}\gamma_{\mathrm{E}})}$
$\displaystyle\times\exp\Bigl{[}-8C_{F}S_{\Gamma}(\mu_{F},\mu_{H})-4C_{F}a_{\Gamma}(\mu_{F},\mu_{H})\ln\frac{\mu_{H}^{2}}{\omega^{2}}+8a_{\gamma_{\ell}}(\mu_{F},\mu_{H})\Bigr{]}$
$\displaystyle\times\exp\Bigl{[}16C_{F}S_{\Gamma}(\mu_{F},\mu_{C})-8C_{F}a_{\Gamma}(\mu_{F},\mu_{C})\ln\frac{\omega\mathcal{T}_{2}}{\mu_{C}^{2}}-8a_{\gamma_{\ell}}(\mu_{F},\mu_{C})$
$\displaystyle+kC_{A}\Bigl{(}-2S_{\Gamma}(\mu_{F},\mu_{C})+a_{\Gamma}(\mu_{F},\mu_{C})\ln\frac{\nu_{F}\mathcal{T}_{2}}{\mu_{C}^{2}}+a_{\Gamma}(\mu_{F},M)\ln\frac{\nu_{F}}{\nu_{C}}\Bigl{)}\Bigr{]}$
$\displaystyle\times\exp\Bigl{[}-8C_{F}S_{\Gamma}(\mu_{F},\mu_{S})+4C_{F}a_{\Gamma}(\mu_{F},\mu_{S})\ln\frac{\mathcal{T}_{2}^{2}}{\mu_{S}^{2}}$
$\displaystyle+kC_{A}\Bigl{(}2S_{\Gamma}(\mu_{F},\mu_{S})+a_{\Gamma}(\mu_{F},\mu_{S})\ln\frac{\mu_{S}^{2}}{\nu_{F}\mathcal{T}_{2}}-a_{\Gamma}(\mu_{F},M)\ln\frac{\nu_{F}}{\nu_{S}}\Bigr{)}\Bigr{]}$
$\displaystyle\times\mathrm{Tr}\Bigl{\\{}\mathbf{H}_{e}^{(0)}\exp\Bigl{[}a_{\Gamma}(\mu_{S},\mu_{H})\mathbf{M}_{e}^{\dagger}\Bigr{]}\mathbf{S}_{e}^{(0)k}\exp\Bigl{[}a_{\Gamma}(\mu_{S},\mu_{H})\mathbf{M}_{e}\Bigr{]}\Bigr{\\}},$
(97)

where we choose an even number of nonsinglets out of $T^{e}$, $T^{f}$, $T_{g}$
and $T^{h}$, which we denote as $k$. In this case, $a_{k}$ is given as

$a_{k}=(8C_{F}-kC_{A})a_{\Gamma}(\mu_{C},\mu_{S}).$ (98)

The numerical factor in front comes from the product of the color factors
$k_{i}^{a}$ in eq. (3) and the muon projections in eq. (46).

Here all the matrices $\mathbf{M}_{e}$, $\mathbf{H}_{e}^{(0)}$ and
$\mathbf{S}^{(0)}$ are $2\times 2$ matrices. The matrix $\mathbf{M}_{e}$ at
threshold is written as

$\mathbf{M}_{e}=\begin{pmatrix}\displaystyle\Bigl{(}4C_{F}-C_{A}\Bigr{)}\ln\frac{n_{13}}{1-n_{13}}+C_{A}\ln(1-n_{13})&&\displaystyle
2\ln\frac{n_{13}}{1-n_{13}}\\\ \displaystyle-
C_{F}\Bigl{(}2C_{F}-C_{A}\Bigr{)}\ln\frac{n_{13}}{1-n_{13}}&&0\end{pmatrix}.$
(99)

The hard function $\mathbf{H}_{e}^{(0)}$ is given as

$\mathbf{H}_{e}^{(0)}=\begin{pmatrix}(1-n_{13})^{2}&0\\\ 0&0\end{pmatrix},$
(100)

and the soft functions can be found in ref. Chay:2021arz . The tree-level soft
factors are given as follows: For the singlet with $k=0$,

$\mathbf{S}_{e}^{(0)}(0,0,0,0)=\begin{pmatrix}C_{A}C_{F}/2&0\\\
0&C_{A}^{2}\end{pmatrix},$ (101)

and there are six possible nonsinglet soft functions with $k=2$ as

$\displaystyle\mathbf{S}_{e}^{(0)}(1,1,0,0)$
$\displaystyle=\frac{1}{2}\begin{pmatrix}\displaystyle
C_{F}-\frac{C_{A}}{2}&0\\\ 0&2C_{A}\end{pmatrix}(12),\
\mathbf{S}_{e}^{(0)}(0,0,1,1)=\frac{1}{2}\begin{pmatrix}\displaystyle
C_{F}-\frac{C_{A}}{2}&0\\\ 0&2C_{A}\end{pmatrix}(34),$
$\displaystyle\mathbf{S}_{e}^{(0)}(1,0,1,0)$
$\displaystyle=\frac{1}{2}\begin{pmatrix}\displaystyle
2C_{F}-\frac{C_{A}}{2}&1\\\ 1&0\end{pmatrix}(13),\
\mathbf{S}_{e}^{(0)}(0,1,0,1)=\frac{1}{2}\begin{pmatrix}\displaystyle
2C_{F}-\frac{C_{A}}{2}&1\\\ 1&0\end{pmatrix}(24),$
$\displaystyle\mathbf{S}_{e}^{(0)}(1,0,0,1)$
$\displaystyle=\frac{1}{2}\begin{pmatrix}2C_{F}-C_{A}&1\\\
1&0\end{pmatrix}(14),\
\mathbf{S}_{e}^{(0)}(0,1,1,0)=\begin{pmatrix}\displaystyle
C_{F}-\frac{C_{A}}{2}&\displaystyle\frac{1}{2}\\\
\displaystyle\frac{1}{2}&0\end{pmatrix}(23).$ (102)

Finally the soft function with $k=4$ is given by

$\displaystyle\mathbf{S}_{e}^{(0)}(1,1,1,1)$
$\displaystyle=\frac{1}{2}\begin{pmatrix}\displaystyle
C_{F}-\frac{C_{A}}{2}&0\\\
1&0\end{pmatrix}(1243)+\frac{1}{2}\begin{pmatrix}\displaystyle
C_{F}-\frac{C_{A}}{2}&1\\\ 0&0\end{pmatrix}(1342)$
$\displaystyle+\begin{pmatrix}\displaystyle(C_{F}-\frac{C_{A}}{2})^{2}&\displaystyle
C_{F}-\frac{C_{A}}{2}\\\ \displaystyle
C_{F}-\frac{C_{A}}{2}&1\end{pmatrix}(12)(34)+\begin{pmatrix}\displaystyle\frac{1}{4}&0\\\
0&0\end{pmatrix}(13)(24).$ (103)

Figure 3: The contributions to the 2-jettiness at NLL from the singlet (blue),
the nonsinglet (red), and the total (green) contributions in (a)
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, (b)
$\ell_{e}\overline{\ell}_{e}\rightarrow\mu^{-}\mu^{+}$ near threshold. The
bands show the theoretical uncertainties.

The plots of the 2-jettiness near threshold at $\theta=\pi/2$ are shown in
fig. 3. The singlet, the nonsinglet, and the total contributions are shown
respectively with the red, blue, and green curves with the bands showing the
theoretical uncertainties at NLL. By normalizing with the tree-level
quantities, the uncertainties from the beam functions for the gauge bosons and
the electrons in the initial state may be alleviated because there is a
partial cancellation at higher orders.

As we can see in fig. 3, the nonsinglet contribution in
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ is comparable to that of the
singlet contribution, while the nonsinglet contribution in
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$
amounts to about 20% of the singlet contribution. This result is puzzling at
first sight because the nonsinglet beam functions for the gauge bosons do not
contribute to the 2-jettiness in
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, resulting in fewer nonsinglet
channels, compared to
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$.
However, the nonsinglet contribution is suppressed compared to the singlet
contribution in
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$
because there is a partial cancellation in summing over the possible
nonsinglet contributions. This can be seen in eqs. (9) and (9), the upper left
matrix elements are proportional to $C_{F}-C_{A}/2=3/4-1<0$, while there is no
such partial cancellation in $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$.

Since we do not have information on the beam functions beyond the required
accuracy, we do not try to estimate the total contribution to the 2-jettiness
by adding the contributions from both channels. However, it is clear that the
nonsinglet contributions contribute to the 2-jettiness in the weak interaction
significantly in both channels, compared to QCD, in which there is a singlet
contribution only. Especially, the nonsinglet contribution is significant in
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$.

##  10 Conclusions

We have analyzed the 2-jettiness in high-energy electroweak scattering in
which muon and anti-muon dijets are observed. The underlying processes consist
of $\ell_{e}\overline{\ell}_{e},WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$.
The central issue in considering the 2-jettiness in weak interaction stems
from the main difference between the electroweak and the QCD processes. That
is, there exist gauge nonsinglets (electrons, muons or neutrinos) in weak
interaction, while only the color singlets, or the hadrons are observed in
strong interaction. This difference results in many interesting aspects, which
do not appear in QCD, and exhibits more intricate structure of the
factorization.

In inclusive quantities in QCD, the Sudakov logarithm from virtual corrections
is cancelled by that from real contributions. For this reason, the PDF is free
of the Sudakov logarithm and it satisfies the DGLAP evolution equation.
However, if we consider exclusive processes, the phase space in the virtual
contribution does not coincide with that in the real contribution. Therefore
the residual Sudakov logarithm survives. As a result, the beam function with
the small lightcone momentum induces the Sudakov logarithm and obeys a
different evolution equation. This type of the Sudakov logarithm appears both
in QCD and in weak interaction.

In weak interaction, additional Sudakov logarithms appear even in inclusive
quantities because the observed particles include gauge nonsinglets. In this
case, the group theory factors in the virtual contributions and in the real
contributions are different. Therefore even for inclusive quantities, the sum
of the virtual and real contributions does not cancel due to the different
group theory factors and it brings out the Sudakov logarithm. It is known as
the Block-Nordsieck violation in electroweak processes. As a result, the PDF
for nonsinglets contain the Sudakov logarithm and does not satisfy the DGLAP
equation, but a different evolution equation. Furthermore, the factorization
in weak interaction may be violated due to the Block-Nordsieck violation
Glauber exchange between spectator partons may violate factorization when the
weak charges of the final states are specified Baumgart:2018ntv . The possible
breakdown of the factorization may start at order $\alpha^{4}$ of the
magnitude $\sim\alpha^{4}\ln^{4}(M^{2}/Q^{2})$, and it is due to the fact that
the group-theory factors for the exchange of two Glauber gauge bosons in
different configurations across the unitarity cuts are different and the
overall effects do not cancel. This should be considered seriously in
ascertaining the factorization in electroweak interaction, but it is beyond
the scope of this paper, and has not been considered here.

The nonsinglet beam functions contribute to the 2-jettiness, as well as the
singlet contributions. Due to the Sudakov logarithm and the existence of the
rapidity divergence, they obey coupled RG equations, and the evolution is
distinct from that of the singlet beam functions. It is also interesting to
note that the additional contribution to the nonsinglet contributions is
proportional to $C_{A}$, whether it is the beam function (or the PDF) for the
gauge bosons or for the lepton. And the matching coefficients which relate the
PDF and the beam functions for the singlets and the nonsinglets are
proportional to each other, both for the leptons and the gauge bosons. The
nonsinglet soft functions get additional contributions, which are also
proportional to $C_{A}$, and they also satisfy the coupled RG equations.

We have considered the singlet and the nonsinglet contributions to the
2-jettiness in both channels
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$ and
$WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$. Their contributions are
comparable in $WW\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, while the
nonsinglet contribution is relatively small in
$\ell_{e}\overline{\ell}_{e}\rightarrow\ell_{\mu}\overline{\ell}_{\mu}$, but
appreciable. The nonsinglet contribution plays a significant role in the
2-jettiness, and possibly in other jet shape observables.

In conclusion, we have established a consistent factorization for the
2-jettiness including the nonsinglet contributions, and their contributions
are numerically significant. This result is in contrast to QCD, in which there
are only singlet contributions. If we have more information on the beam
function (or PDF) for the weak gauge bosons and the electrons, the prediction
can be more precise.

We admit that the 2-jettiness analyzed here is far from realistic experimental
comparison because we confine the gauge group to SU(2). But the main focus
here is to present clearly the difference between QCD and the weak
interaction. And it turns out that the nonsinglet contributions in electroweak
high-energy scattering is appreciable, in comparison to QCD. Of course, we
have to choose the electroweak gauge group SU(2)$\times$U(1), to be realistic.
The phenomenology with the electroweak gauge group will be studied in the near
future.

##  Appendix A NLO calculation of the gauge-boson beam function and PDF

###  A.1 Beam functions for the gauge bosons

The Feynman rules for the gauge boson $\mathcal{B}_{n\perp}^{\mu c}$ get
complicated due to the insertion of the rapidity regulator. The Feynman rules
are shown in fig. 4. At NLO, the Feynman diagrams, contributing to the beam
function, are presented in fig. 5. Figure 5 (f) represents the mixing between
the fermion and the gauge boson, but it does not possess any UV divergence,
hence the beam function does not have any mixing. In obtaining the final
result, we take the limit of small mass $M$, and keep only the logarithmic
terms in $M$. We also refer to fig. 4 for the Feynman diagrams of the PDFs
with the caveat that $t$ is not measured for the PDFs.

Figure 4: Feynman rules for $\mathcal{B}_{n\perp}^{\mu c}$ with one and two
gauge bosons.

Let us first consider the singlet matrix elements which are proportional to
$(G^{0})^{ab}$. Only the group theory factors differ in the nonsinglet matrix
elements. The naive contribution from fig. 5(a) is written as

$\displaystyle\tilde{M}_{(a)}(G^{0})^{ab}=-\mu^{2\epsilon}_{\overline{\mathrm{MS}}}\int\frac{d^{D}{\ell}}{(2\pi)^{D}}\epsilon_{\alpha}^{a}\epsilon_{\beta}^{b*}\frac{1}{(\ell^{2}-M^{2})^{2}}V_{3\beta\rho\sigma}^{bde}(-p,\ell,p-\ell)V_{3\alpha\sigma\nu}^{aec}(p,\ell,-\\-\ell)$
(104) $\displaystyle\times
B_{0}^{\mu\rho}(-\ell)B_{0\mu\nu}(\ell)(2\pi)\delta\Bigl{(}(p-\ell)^{2}-M^{2}\Bigr{)}\delta\Bigl{(}\frac{t}{\omega}-(p^{+}-\ell^{+})\Bigr{)}\delta(\ell^{-}-\omega)\theta(\omega)\theta(t)\theta(p^{-}-\ell^{-}),$

where we average over the possible $(D-2)$ polarization, but not over the weak
charges because we consider the beam function with a fixed weak charge.

$\frac{1}{D-2}\sum_{\mathrm{pol.}}\epsilon_{\alpha}\epsilon_{\beta}^{*}=-\frac{g_{\perp\alpha\beta}}{D-2}.$
(105)

And $V_{3\mu\nu\rho}^{abc}(p_{1},p_{2},p_{3})$ is the three-boson vertex,
which is given as

$V_{3\mu\nu\rho}^{abc}(p_{1},p_{2},p_{3})=gf^{abc}\Bigl{[}g^{\mu\nu}(p_{1}-p_{2})^{\rho}+g^{\nu\beta}(p_{2}-p_{3})^{\mu}+g^{\rho\mu}(p_{3}-p_{1})^{\nu}\Bigr{]}.$
(106)

Figure 5: Feynman diagrams for the beam functions and the PDFs. The wavy lines
with solid lines denote collinear gauge bosons, and the solid lines are
leptons. The mirror images of (b), (d) and (e) are omitted. For the beam
functions, the virtuality $t$ and the longitudinal momentum fraction $z$ are
measured, but only $z$ is measured for the PDFs.

Since $\tilde{M}_{(a)}$ is finite and we can put $\epsilon=\eta=0$ and perform
the integration. After some algebra, $\tilde{M}_{(a)}$ is given as

$\displaystyle\tilde{M}_{(a)}$ $\displaystyle=-\frac{\alpha
C_{A}}{4\pi}\frac{1}{x(t-x^{2}M^{2})^{2}}\Bigl{[}2t(-2+2x-3x^{2}+2x^{3})+M^{2}x(4+7x^{2}+2x^{3})\Bigr{]}$
$\displaystyle\longrightarrow\frac{\alpha
C_{A}}{2\pi}\Bigl{[}\Bigl{(}2\frac{1-x}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}\Bigr{)}\delta(t)$
$\displaystyle+2\Bigl{(}\frac{1-x}{x}+x(1-x)+\frac{x}{2}\Bigr{)}\Bigl{(}\delta(t)\ln\frac{(1-x)\mu^{2}}{x(1-x+x^{2})M^{2}}+\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\Bigr{)}\Bigr{]}.$
(107)

In the small mass limit, it is regarded as the distribution in $t$. The
coefficient of $\delta(t)$ is obtained by integrating $\tilde{M}_{(a)}$ over
$t$. The remainder is determined by the fact that $\tilde{M}_{(a)}$ is
independent of $\mu$. The zero-bin contribution is suppressed, hence neglected
at leading order in SCET. The functions $\mathcal{L}_{n}(x)$ are defined as

$\mathcal{L}_{n}(x)\equiv\Bigl{[}\frac{\theta(x)\ln^{n}x}{x}\Bigr{]}_{+}=\lim_{\beta\rightarrow
0}\Bigl{[}\frac{\theta(x-\beta)\ln^{n}x}{x}+\delta(x-\beta)\frac{\ln^{n+1}\beta}{n+1}\Bigr{]}.$
(108)

The naive contribution $\tilde{M}_{(b)}$ from fig. 5(b) is given as

$\displaystyle\tilde{M}_{(b)}(G^{0})^{ab}=\mu^{2\epsilon}_{\overline{\mathrm{MS}}}\int\frac{d^{D}{\ell}}{(2\pi)^{D}}\frac{g_{\perp\alpha\beta}}{D-2}B_{1,bcd}^{\mu\rho\beta}(p-\ell,-p)\delta(\ell^{-}-\omega)\delta(t+\omega(p^{+}-\ell^{+}))$
$\displaystyle\times(2\pi)\delta((p-\ell)^{2}-M^{2})V_{3\alpha\rho\nu}^{adc}(p,\ell-p,-\ell)\frac{1}{\ell^{2}-M^{2}}B_{0}^{\mu\nu}(\ell)\theta(\omega)\theta(p^{-}-\ell^{-})\theta(t)$
$\displaystyle=\frac{\alpha
C_{A}}{4\pi}(G^{0})^{ab}\frac{(\mu^{2}e^{\gamma_{\mathrm{E}}})^{\epsilon}}{\Gamma(1-\epsilon)}\int
d\ell^{+}d\ell^{-}d\bm{\ell}_{\perp}^{2}(\bm{\ell}_{\perp}^{2})^{-\epsilon}\frac{\ell^{-}+p^{-}}{\ell^{-}-p^{-}}\frac{1}{\ell^{2}-M^{2}}\Bigl{(}\frac{\nu^{2}}{\ell^{-}(p^{-}-\ell^{-})}\Bigr{)}^{\eta}$
$\displaystyle\times\delta(\ell^{-}-\omega)\delta(\ell^{+}-(p^{+}-b^{+}))\delta\Bigl{(}\bm{\ell}_{\perp}^{2}-\Bigl{(}\frac{(1-x)t}{x}-M^{2}\Bigr{)}\Bigr{)}\theta(\omega)\theta(p^{-}-\omega)$
$\displaystyle=\frac{\alpha
C_{A}}{4\pi}(G^{0})^{ab}\frac{(\mu^{2}e^{\gamma_{\mathrm{E}}})^{\epsilon}}{\Gamma(1-\epsilon)}\frac{1+x}{1-x}\Bigl{(}\frac{(1-x)t}{x}-M^{2}\Bigr{)}^{-\epsilon}\frac{x}{t-x^{2}M^{2}}\Bigl{(}\frac{\nu^{2}}{x(1-x)(p^{-})^{2}}\Bigr{)}^{\eta}\theta(t).$
(109)

In computing the group theory factor, we have used the fact that
$f^{adc}f^{bcd}=-\delta^{ab}/N$.

Since it does not possess any divergence, we can put $\epsilon=\eta=0$, and it
is written as

$\tilde{M}_{(b)}=\frac{\alpha
C_{A}}{4\pi}\frac{x(1+x)}{1-x}\frac{1}{t-x^{2}M^{2}}\theta\Bigl{(}\frac{(1-x)t}{x}-M^{2}\Bigr{)}\theta(x)\theta(1-x).$
(110)

In the limit of small $M$, it can be written as a double distribution in terms
of $t$ and $x$. In order to extract the coefficient for $\delta(t)$, we first
look at the integral

$\frac{x(1+x)}{1-x}\int_{xM^{2}/(1-x)}^{\mu^{2}}\frac{dt}{t-x^{2}M^{2}}=\frac{x(1+x)}{1-x}\ln\frac{(1-x)(\mu^{2}-M^{2}x^{2})}{x(1-x+x^{2})M^{2}},$
(111)

where $\mu$ is chosen to be an arbitrary scale. The right-hand side of eq.
(111) can be expressed in terms of the distributions by considering the
dependence on $x$ as

$\displaystyle\frac{x(1+x)}{1-x}\ln\frac{(1-x)(\mu^{2}-M^{2}x^{2})}{x(1-x+x^{2})M^{2}}$
$\displaystyle=A\delta(1-x)+x(1+x)\mathcal{L}_{0}(1-x)\ln\frac{\mu^{2}}{x(1-x+x^{2})M^{2}}$
$\displaystyle+x(1+x)\mathcal{L}_{1}(1-x),$ (112)

where $A$ is a constant to be determined. Integrating the left-hand side of
eq. (A.1) yields

$\int_{0}^{k}dx\frac{x(1+x)}{1-x}\ln\frac{(1-x)(\mu^{2}-M^{2}x^{2})}{x(1-x+x^{2})M^{2}}=-\frac{9}{2}+\frac{5\sqrt{3}}{6}+\frac{4}{9}\pi^{2}-\frac{5}{2}\ln\frac{\mu^{2}}{M^{2}}+\ln^{2}\frac{\mu^{2}}{M^{2}},$
(113)

with $k=\mu^{2}/(\mu^{2}+M^{2})$, and integrating the right-hand side yields

$\displaystyle\int_{0}^{1}dx\Bigl{[}A\delta(1-x)+x(1+x)\mathcal{L}_{0}(1-x)\ln\frac{\mu^{2}}{x(1-x+x^{2})M^{2}}+x(1+x)\mathcal{L}_{1}(1-x)\Bigr{]}$
$\displaystyle=A-\frac{9}{2}+\frac{5\sqrt{3}}{6}+\frac{4}{9}\pi^{2}-\frac{5}{2}\ln\frac{\mu^{2}}{M^{2}}.$
(114)

By comparing two quantities, we obtain $A=\ln^{2}(\mu^{2}/M^{2})$.

As a result, we can write $\tilde{M}_{(b)}$ as

$\displaystyle\tilde{M}_{(b)}=\frac{\alpha
C_{A}}{4\pi}\Bigl{\\{}\delta(t)\Bigl{[}\delta(1-x)\ln^{2}\frac{\mu^{2}}{M^{2}}+x(1+x)\mathcal{L}_{0}(1-x)\ln\frac{\mu^{2}}{x(1-x+x^{2})M^{2}}$
$\displaystyle+x(1+x)\mathcal{L}_{1}(1-x)\Bigr{]}+f_{b}(x,\mu)\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}+g_{b}(x,\mu)\frac{1}{\mu^{2}}\mathcal{L}_{1}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\Bigr{\\}}.$
(115)

The unknown functions $f_{b}(x,\mu)$ and $g_{b}(x,\mu)$ are determined by the
requirement that $\tilde{M}_{(b)}$ is independent of $\mu^{2}$, that is,
$d\tilde{M}_{(b)}/d\ln\mu^{2}=0$. They are given as

$f_{b}(x,\mu)=2\delta(1-x)\ln\frac{\mu^{2}}{M^{2}}+x(1+x)\mathcal{L}_{0}(1-x),\
\ g_{b}(x)=2\delta(1-x).$ (116)

The final result for $\tilde{M}_{(b)}$ is written as

$\displaystyle\tilde{M}_{(b)}$ $\displaystyle=\frac{\alpha
C_{A}}{4\pi}\Bigl{\\{}\delta(t)\Bigl{[}\delta(1-x)\ln^{2}\frac{\mu^{2}}{M^{2}}+x(1+x)\mathcal{L}_{0}(1-x)\ln\frac{\mu^{2}}{x(1-x+x^{2})M^{2}}$
$\displaystyle+x(1+x)\mathcal{L}_{1}(1-x)\Bigr{]}+\Bigl{(}2\ln\frac{\mu^{2}}{M^{2}}\delta(1-x)+x(1+x)\mathcal{L}_{0}(1-x)\Bigr{)}\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}$
$\displaystyle+2\delta(1-x)\frac{1}{\mu^{2}}\mathcal{L}_{1}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\Bigr{\\}}.$
(117)

From eq. (A.1), the zero-bin contribution $M_{(b)}^{\varnothing}$ can show up
when $p-\ell$ or $\ell$ become soft. However, the contribution for the soft
$\ell$ is suppressed and only the case with soft $p-\ell$ contributes. It is
given, with $b^{+}=t/\omega$, as

$\displaystyle M_{(b)}^{\varnothing}$ $\displaystyle=\frac{\alpha
C_{A}}{4\pi}\frac{(\mu^{2}e^{\gamma_{\mathrm{E}}})^{\epsilon}}{\Gamma(1-\epsilon)}\delta(1-x)\int_{M^{2}/b^{+}}^{\infty}d\ell^{-}\frac{2}{b^{+}p^{-}\ell^{-}}(\ell^{-}b^{+}-M^{2})^{-\epsilon}\Bigl{(}\frac{\nu}{\ell^{-}}\Bigr{)}^{\eta}$
$\displaystyle=\frac{\alpha
C_{A}}{2\pi}\delta(1-x)\Bigl{(}\frac{\mu^{2}e^{\gamma_{\mathrm{E}}}}{M^{2}}\Bigr{)}^{\epsilon}\frac{\Gamma(\epsilon+\eta)}{\Gamma(1+\eta)}\Bigl{(}\frac{\nu^{3}}{\omega
M^{2}}\Bigr{)}^{\eta}\frac{\nu^{-2\eta}}{t^{1-\eta}}$
$\displaystyle=\frac{\alpha
C_{A}}{2\pi}\delta(1-x)\Bigl{\\{}\delta(t)\Bigl{[}\Bigl{(}\frac{1}{\eta}+\ln\frac{\nu}{p^{-}}\Bigr{)}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}-\frac{1}{\epsilon^{2}}+\frac{1}{2}\ln^{2}\frac{\mu^{2}}{M^{2}}+\frac{\pi^{2}}{12}\Bigr{]}$
$\displaystyle+\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}\frac{1}{\mu^{2}}\mathcal{L}_{0}\Bigl{(}\frac{t}{\mu^{2}}\Bigr{)}\Bigr{\\}}.$
(118)

The naive contribution $\tilde{M}_{(d)}$ is given as

$\displaystyle\tilde{M}_{(d)}(G^{0})^{ab}=-\mu^{2\epsilon}_{\overline{\mathrm{MS}}}\frac{g_{\perp\alpha\beta}}{2(D-2)}\int\frac{d^{D}{\ell}}{(2\pi)^{D}}B_{0}^{\mu\beta}(-p)B_{1,\mu\sigma\rho}^{bcd}(\ell,p-\ell)V_{3adc}^{\alpha\rho\sigma}(p,\ell-p,-\ell)$
$\displaystyle\times\frac{1}{\ell^{2}-M^{2}}\frac{1}{(p-\ell)^{2}-M^{2}}\delta(p^{-}-\omega)\delta(t/\omega)$
$\displaystyle=ig^{2}\frac{(G^{0})^{ab}C_{A}}{2}\mu^{2\epsilon}_{\overline{\mathrm{MS}}}\delta(1-x)\delta(t)\int\frac{d^{D}{\ell}}{(2\pi)^{D}}\frac{1}{\ell^{2}-M^{2}}\frac{1}{(p-\ell)^{2}-M^{2}}$
$\displaystyle\times\Bigl{[}\frac{\ell^{-}+p^{-}}{p^{-}-\ell^{-}}\Bigl{(}\frac{\nu}{p^{-}-\ell^{-}}\Bigr{)}^{\eta}+\frac{2p^{-}-\ell^{-}}{\ell^{-}}\Bigl{(}\frac{\nu}{\ell^{-}}\Bigr{)}^{\eta}\Bigr{]}$
$\displaystyle=ig^{2}(G^{0})^{ab}C_{A}\mu^{2\epsilon}_{\overline{\mathrm{MS}}}\delta(1-x)\delta(t)\int\frac{d^{D}{\ell}}{(2\pi)^{D}}\frac{1}{\ell^{2}-M^{2}}\frac{1}{(p-\ell)^{2}-M^{2}}\frac{\ell^{-}+p^{-}}{p^{-}-\ell^{-}}\Bigl{(}\frac{\nu}{p^{-}-\ell^{-}}\Bigr{)}^{\eta}$
$\displaystyle=-\frac{\alpha
C_{A}}{4\pi}(G^{0})^{ab}\delta(1-x)\delta(t)\Gamma(\epsilon)\Bigl{(}\frac{\mu^{2}}{M^{2}}\Bigr{)}^{\epsilon}\Bigl{(}\frac{\nu}{p^{-}}\Bigr{)}^{\eta}\int_{0}^{1}dy\frac{1+y}{(1-y)^{1+\eta}}(1-y+y^{2})^{-\epsilon}.$
(119)

The zero-bin contribution $M_{(d)}^{\varnothing}$ comes from the region in
which $p^{\mu}-\ell^{\mu}$ becomes soft, while the contribution from the
region in which $\ell^{\mu}$ becomes soft is suppressed. It is given as

$\displaystyle M_{(d)}^{\varnothing}$ $\displaystyle=\frac{\alpha
C_{A}}{4\pi}\delta(1-x)\delta(t)e^{\gamma_{\mathrm{E}}\epsilon}\Bigl{(}\frac{\mu^{2}}{M^{2}}\Bigr{)}^{\epsilon}\Bigl{(}\frac{\nu}{M}\Bigr{)}^{\eta}\frac{\Gamma(\epsilon+\eta/2)}{\Gamma(1+\eta/2)}\int_{0}^{\infty}y^{-1+\eta}.$
(120)

The net contribution $M_{(d)}=\tilde{M}_{(d)}-M_{(d)}^{\varnothing}$ is
computed and the result is expanded near $\eta=0$ first, and then near
$\epsilon=0$, and the result is given as

$M_{(d)}=\frac{\alpha
C_{A}}{4\pi}\delta(1-x)\delta(t)\Bigl{[}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}\Bigl{(}\frac{2}{\eta}+2\ln\frac{\nu}{p^{-}}+1\Bigr{)}+2-\frac{\pi}{\sqrt{3}}-\frac{\pi^{2}}{9}\Bigr{]}.$
(121)

In fig. 5(e), it involves the self-energy of the gauge boson and the blob
contains the loops of the gauge bosons, the ghosts, the fermions and the
scalar particles. The Feynman diagrams are shown in figs. 6 and 7. Contrary to
QCD, we compute the self-energy corrections with the gauge boson mass $M$ with
the left-handed fermions only. And the additional contributions come from the
interacton of the complex scalar multiplets.

Figure 6: Feynman diagrams for the self-energy of the gauge bosons from the
loops of (a, b) the gauge bosons, (c) the ghost loops, and (d) the fermions.

The self-energy of the gauge boson can be written as

$\Pi(p^{2})\delta^{ab}\Bigl{(}g^{\mu\nu}-\frac{p^{\mu}p^{\nu}}{p^{2}}\Bigr{)},$
(122)

and the field-strength renormalization $Z$ is given as

$Z=\Bigl{(}1-\frac{d\Pi(p^{2})}{dp^{2}}\Big{|}_{p^{2}=M^{2}}\Bigr{)}^{-1}.$
(123)

The field-strength renormalization $Z_{W}$ from the gauge bosons and the ghost
particles, common to all the SU($N$) gauge interactions, comes from fig. 6 (a)
to (c). At order $\alpha$, it is given as

$Z_{W}^{(1)}=\frac{\alpha
C_{A}}{4\pi}\Bigl{[}\frac{5}{3}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}-\frac{5}{9}+\frac{\pi}{\sqrt{3}}\Bigr{]}.$
(124)

The contribution $Z_{f}$ from the fermion loop in fig. 6 (d) at order $\alpha$
is given by

$Z_{f}^{(1)}=\frac{\alpha}{4\pi}\frac{1}{2}n_{f}T_{F}\Bigl{[}-\frac{4}{3}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}-\frac{8}{9}\Bigr{]}.$
(125)

The factor 1/2 is to remind that only the left-handed fermions contribute.

Figure 7: Feynman diagrams for the self-energy of the gauge bosons from the
scalar particles. (a) contains a divergence, (b) is finite, and (c) does not
contribute to the field-strength renormalization.

The remaining contribution comes from the scalar particles. For the SU(2)
gauge group, we follow the standard scalar fields in the Standard Model.
However, for general SU($N$) gauge theories, the symmetry breaking pattern may
depend on the model. We just assume that all the gauge bosons have a common
mass $M$, as in SU(2). The Feynman diagrams from the scalar particles are
shown in fig. 7. Fig. 7(a) contains an UV divergence, and its contribution is
denoted as $Z_{sa}$ below. Fig. 7(b) depends on the structure of the complex
scalar multiplets, to yield the same masses to all the gauge bosons in the
SU($N$) gauge theory. However, since it is finite and does not contribute to
the field-strength renormalization, we neglect this. Note that this
contribution cancels in obtaining the matching coefficients in eq. (5). And
fig. 7(c) vanishes.

The field-strength renormalization $Z$ at order $\alpha$ is given as

$\displaystyle Z^{(1)}$ $\displaystyle=Z_{W}^{(1)}+Z_{f}^{(1)}+Z_{sa}^{(1)}$
(126)
$\displaystyle=\frac{\alpha}{4\pi}\Bigl{[}(\beta_{0}-2C_{A})\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}+C_{A}(-\frac{5}{9}+\frac{\pi}{\sqrt{3}}\Bigr{)}-\frac{4}{9}n_{f}T_{F}+n_{s}T_{F}\Bigl{(}-\frac{17}{9}+\frac{\pi}{\sqrt{3}}\Bigr{)}\Bigr{]},$

where $\beta_{0}=11C_{A}/3-2n_{f}T_{F}/3-n_{s}T_{F}/3$ is the beta function at
leading order.

Adding all the contributions, the singlet beam function $B_{Ws}^{(1)}$ for the
gauge boson at NLO is written as

$\displaystyle
B_{Ws}^{(1)}(t,x,\mu)=M_{(a)}+2(M_{(b)}-M_{(b)}^{\varnothing}+M_{(d)}-M_{(d)}^{\varnothing})+Z^{(1)}\delta(1-x)\delta(t)$
$\displaystyle=\frac{\alpha}{2\pi}\Bigl{\\{}\delta(1-x)\Bigl{[}C_{A}\Bigl{(}\frac{2}{\epsilon^{2}}\delta(t)-\frac{2}{\epsilon}\frac{1}{\mu^{2}}\mathcal{L}_{0}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}\Bigr{)}+\frac{\beta_{0}}{2}\frac{1}{\epsilon}\Bigr{]}$
$\displaystyle+\delta(t)\Bigl{(}\frac{\beta_{0}}{2}\delta(1-x)+C_{A}P_{WW}(x)\Bigr{)}\ln\frac{\mu^{2}}{M^{2}}+C_{A}\delta(t)\Bigl{(}P_{WW}(x)\ln\frac{1-x}{x}-\frac{\pi^{2}}{6}\delta(1-x)\Bigr{)}$
$\displaystyle+C_{A}\Bigl{[}\delta(1-x)\frac{2}{\mu^{2}}\mathcal{L}_{1}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}+P_{WW}(x)\frac{1}{\mu^{2}}\mathcal{L}_{0}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}\Bigr{]}$
$\displaystyle+C_{A}\delta(t)\Bigl{[}\delta(1-x)\Bigl{(}\frac{31}{18}-\frac{\pi^{2}}{9}-\frac{\pi}{2\sqrt{3}}\Bigr{)}-\Bigl{(}\frac{2(1-x)}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}\Bigr{)}$
$\displaystyle-
P_{WW}(x)\ln(1-x+x^{2})\Bigr{]}+\delta(1-x)\delta(t)\Bigl{(}-\frac{2}{9}n_{f}T_{F}+\frac{1}{2}n_{s}T_{F}\Bigl{(}-\frac{17}{9}+\frac{\pi}{\sqrt{3}}\Bigr{)}\Bigr{\\}}.$
(127)

For the nonsinglet beam functions, only the group theory factors are
different. Extracting the nonsinglet matrix elements, proportional to
$d^{abc}$, the matrix elements $\tilde{M}_{(a)}$, $\tilde{M}_{(b)}$ and
$M_{(b)}^{\varnothing}$ in the singlet calculation are replaced by
$\tilde{M}_{(a)}/2$, $\tilde{M}_{(b)}/2$ and $M_{(b)}^{\varnothing}/2$
respectively, while the remaining matrix elements are the same. Therefore, the
nonsinglet beam functions at NLO are given as

$\displaystyle
B_{Wn}^{(1)}(t,x,\mu,\nu)=\frac{1}{2}\tilde{M}_{(a)}+\tilde{M}_{(b)}-M_{(b)}^{\varnothing}+2(\tilde{M}_{(d)}-M_{(d)}^{\varnothing})+Z^{(1)}\delta(1-x)\delta(t)$
$\displaystyle=B_{Ws}^{(1)}(t,x,\mu)-\frac{1}{2}\Bigl{(}\tilde{M}_{(a)}+2(\tilde{M}_{(b)}-M_{(b)}^{\varnothing})\Bigr{)}$
$\displaystyle=B_{Ws}^{(1)}(t,x,\mu)$ $\displaystyle+\frac{\alpha
C_{A}}{4\pi}\Bigl{\\{}2\delta(1-x)\Bigl{[}\delta(t)\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}\Bigl{(}\frac{1}{\eta}+\ln\frac{\nu}{p^{-}}\Bigr{)}-\frac{1}{\epsilon^{2}}\delta(t)+\frac{1}{\epsilon}\frac{1}{\mu^{2}}\mathcal{L}_{0}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}-\frac{1}{\mu^{2}}\mathcal{L}_{1}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}\Bigr{]}$
$\displaystyle-
P_{WW}(x)\Bigl{(}\delta(t)\ln\frac{\mu^{2}}{M^{2}}+\frac{1}{\mu^{2}}\mathcal{L}_{0}\bigl{(}\frac{t}{\mu^{2}}\bigr{)}\Bigr{)}+\delta(t)\Bigl{(}P_{WW}(x)\ln\frac{x(1-x+x^{2})}{1-x}+\frac{\pi^{2}}{6}\delta(1-x)$
$\displaystyle+2\frac{1-x}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}\Bigr{)}\Bigr{\\}}.$
(128)

###  A.2 PDF for the gauge bosons

We can compute the PDF either by computing the matrix elements in eq. (36), or
by integrating the results of the beam functions with respect to the variable
$t$. This can be seen by comparing the definitions of the beam functions in
eq. (5) and those of the PDFs in eq. (36). We adopt the second approach here.

We first consider the matrix elements for the singlet PDF. The naive
contribution from fig. 5(a) can be obtained from eq. (A.1) by integrating with
respect to $t$, and it is given as

$\displaystyle\tilde{M}_{(a)}$ $\displaystyle=\frac{\alpha
C_{A}}{2\pi}\Bigl{[}2\Bigl{(}\frac{1-x}{x}+\frac{3x}{2}-x^{2}\Bigr{)}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{(1-x+x^{2})M^{2}}\Bigr{)}$
$\displaystyle+\frac{2(1-x)}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}\Bigr{]}.$
(129)

Integrating eq. (A.1) with respect to $t$ yields $\tilde{M}_{b}$ for the PDF.
However, the correct rapidity divergence is obtained only after the zero-bin
subtraction is performed. We present the true collinear contribution
$M_{b}=\tilde{M}_{b}-M^{\varnothing}_{b}$ as

$\displaystyle M_{(b)}$ $\displaystyle=\frac{\alpha
C_{A}}{4\pi}\Bigl{[}-2\Bigl{(}\frac{1}{\eta}+\ln\frac{\nu}{p^{-}}\Bigr{)}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}+x(1+x)\mathcal{L}_{0}(1-x)\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}(1-x+x^{2})}\Bigr{)}\Bigr{]}.$
(130)

In a similar approach to computing $M_{(b)}$, $M_{(d)}$ can be obtained as

$M_{(d)}=\frac{\alpha
C_{A}}{4\pi}\delta(1-x)\Bigl{[}\Bigl{(}\frac{2}{\eta}+2\ln\frac{\nu}{p^{-}}+1\Bigr{)}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}+2-\frac{\pi}{\sqrt{3}}-\frac{\pi^{2}}{9}\Bigr{]}.$
(131)

The wave function renormalization $M_{(e)}$ is given by

$\displaystyle Z^{(1)}$
$\displaystyle=\frac{\alpha}{4\pi}\Bigl{[}(\beta_{0}-2C_{A})\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}+C_{A}\Bigl{(}-\frac{5}{9}+\frac{\pi}{\sqrt{3}}\Bigr{)}-\frac{4}{9}n_{f}T_{F}+n_{s}T_{F}\Bigl{(}-\frac{17}{9}+\frac{\pi}{\sqrt{3}}\Bigr{)}\Bigr{]}.$
(132)

Combining all the matrix elements, the singlet PDF at NLO is given as

$\displaystyle
f_{Ws}^{(1)}(x,\mu,\nu)=\tilde{M}_{a}+2(M_{b}+M_{d})+Z^{(1)}\delta(1-x)$
$\displaystyle=\frac{\alpha}{2\pi}\Bigl{\\{}C_{A}\Bigl{(}P_{WW}(x)+\frac{\beta_{0}}{2}\delta(1-x)\Bigr{)}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}$
$\displaystyle-
C_{A}\Bigl{(}P_{WW}(x)\ln(1-x+x^{2})+2\frac{1-x}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}\Bigr{)}$
$\displaystyle+\delta(1-x)\Bigl{[}C_{A}\Bigl{(}\frac{31}{18}-\frac{\pi}{2\sqrt{3}}-\frac{\pi^{2}}{9}\Bigr{)}-\frac{2}{9}n_{f}T_{F}+n_{s}T_{F}\Bigl{(}-\frac{17}{9}+\frac{\pi}{\sqrt{3}}\Bigr{)}\Bigr{\\}}.$
(133)

The matrix elements for the nonsinglet PDFs have different group theory
factors, and the relation between the singlets and the nonsinglets is the same
as that in the beam functions. The nonsinglet PDFs at NLO are written as

$\displaystyle
f_{Wn}^{(1)}(x,\mu,\nu)=\frac{1}{2}M_{a}+\tilde{M}_{(b)}-M_{(b)}^{\varnothing}+2(\tilde{M}_{(d)}-M_{(d)}^{\varnothing})+Z^{(1)}\delta(1-x)$
$\displaystyle=f_{Ws}^{(1)}(x,\mu)-\frac{1}{2}\Bigl{(}\tilde{M}_{(a)}+2(\tilde{M}_{(b)}-M_{(b)}^{\varnothing})\Bigr{)}$
$\displaystyle=f_{Ws}^{(1)}(x,\mu)+\frac{\alpha}{2\pi}\frac{C_{A}}{2}\Bigl{[}\delta(1-x)\Bigl{(}\frac{2}{\eta}+2\ln\frac{\nu}{p^{-}}\Bigr{)}\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}-\Bigl{(}\frac{1}{\epsilon}+\ln\frac{\mu^{2}}{M^{2}}\Bigr{)}P_{WW}(x)$
$\displaystyle+2\frac{1-x}{x}+2x(1-x)+\frac{3}{2}\frac{x(1-x)}{1-x+x^{2}}+P_{WW}(x)\ln(1-x+x^{2})\Bigr{]}.$
(134)

##  Appendix B Soft matrix elements at tree level

When there are singlets only, the soft matrix is given by

$S^{(0)}(0,0,0,0)=\begin{pmatrix}N(N^{2}-1)&0&0\\\
0&\frac{1}{2}N(N^{2}-1)&0\\\ 0&0&\frac{(N-2)(N+2)(N^{2}-1)}{2N}\end{pmatrix}.$
(135)

When there are 2 nonsinglets, the soft matrices are given as

$\displaystyle S^{(0)}(1,1,0,0)$
$\displaystyle=(N-2)\begin{pmatrix}2(N+2)&0&0\\\ 0&\frac{1}{2}(N+2)&0\\\
0&0&\frac{(N+2)(N^{2}-12)}{2N^{2}}\end{pmatrix}(12),$ $\displaystyle
S^{(0)}(1,0,1,0)$ $\displaystyle=(N-2)\begin{pmatrix}0&0&\frac{N+2}{N}\\\
0&\frac{1}{4}(N+2)&\frac{1}{4}(N+2)\\\
\frac{N+2}{N}&\frac{1}{4}(N+2)&\frac{(N+2)(N^{2}-12)}{4N^{2}}\end{pmatrix}(13),$
$\displaystyle S^{(0)}(1,0,0,1)$
$\displaystyle=(N-2)\begin{pmatrix}0&0&\frac{N+2}{N}\\\
0&\frac{1}{4}(N+2)&-\frac{1}{4}(N+2)\\\
\frac{N+2}{N}&-\frac{1}{4}(N+2)&\frac{(N+2)(N^{2}-12)}{4N^{2}}\end{pmatrix}(14),$
$\displaystyle S^{(0)}(0,1,1,0)$
$\displaystyle=(N-2)\begin{pmatrix}0&0&\frac{N+2}{N}\\\
0&\frac{1}{4}(N+2)&-\frac{1}{4}(N+2)\\\
\frac{N+2}{N}&-\frac{1}{4}(N+2)&\frac{(N+2)(N^{2}-12)}{4N^{2}}\end{pmatrix}(23),$
$\displaystyle S^{(0)}(0,1,0,1)$
$\displaystyle=(N-2)\begin{pmatrix}0&0&\frac{N+2}{N}\\\
0&\frac{1}{4}(N+2)&\frac{1}{4}(N+2)\\\
\frac{N+2}{N}&\frac{1}{4}(N+2)&\frac{(N+2)(N^{2}-12)}{4N^{2}}\end{pmatrix}(24),$
$\displaystyle S^{(0)}(0,0,1,1)$
$\displaystyle=\begin{pmatrix}(N^{2}-1)&0&0\\\ 0&-\frac{1}{2}&0\\\
0&0&-\frac{(N-2)(N+2)}{2N^{2}}\end{pmatrix}(34).$ (136)

Note that, with $N=2$, only the last term $S^{(0)}(0,0,1,1)$ survives, which
correspond to the case with no nonsinglet beam functions and 2 jet functions.

For the case of 3 nonsinglets, they are given as

$\displaystyle S^{(0)}(1,1,1,0)$
$\displaystyle=(N-2)\begin{pmatrix}0&\frac{N+2}{2N}\big{(}(123)-(132)\big{)}&\frac{N^{2}-12}{2N}\frac{(123)+(132)}{N-2}\\\
-\frac{N+2}{2N}\big{(}(123)-(132)\big{)}&-\frac{(123)+(132)}{N-2}&\frac{N+2}{N^{2}}\big{(}(123)-(132)\big{)}\\\
\frac{N^{2}-12}{2N}\frac{(123)+(132)}{N-2}&-\frac{N+2}{N^{2}}\big{(}(123)-(132)\big{)}&-\frac{2(N^{2}-10)}{N^{2}}\frac{(123)+(132)}{N-2}\end{pmatrix},$
$\displaystyle S^{(0)}(1,1,0,1)$
$\displaystyle=(N-2)\begin{pmatrix}0&\frac{N+2}{2N}\big{(}(124)-(142)\big{)}&\frac{N^{2}-12}{2N}\frac{(124)+(142)}{N-2}\\\
-\frac{N+2}{2N}\big{(}(124)-(142)\big{)}&-\frac{(124)+(142)}{N-2}&\frac{N+2}{N^{2}}\big{(}(124)-(142)\big{)}\\\
\frac{N^{2}-12}{2N}\frac{(124)+(142)}{N-2}&-\frac{N+2}{N^{2}}\big{(}(124)-(142)\big{)}&-\frac{2(N^{2}-10)}{N^{2}}\frac{(124)+(142)}{N-2}\end{pmatrix},$
$\displaystyle S^{(0)}(1,0,1,1)$
$\displaystyle=(N-2)\begin{pmatrix}0&0&\frac{N+2}{N}(143)\\\
0&-\frac{1}{2}\frac{(134)+(143)}{N-2}&0\\\
\frac{N+2}{N}(134)&0&-\frac{N^{2}-12}{2N^{2}}\frac{(134)+(143)}{N-2}\end{pmatrix},$
$\displaystyle S^{(0)}(0,1,1,1)$
$\displaystyle=(N-2)\begin{pmatrix}0&0&\frac{N+2}{N}(243)\\\
0&-\frac{1}{2}\frac{(234)+(243)}{N-2}&0\\\
\frac{N+2}{N}(234)&0&-\frac{N^{2}-12}{2N^{2}}\frac{(234)+(243)}{N-2}\end{pmatrix}.$
(137)

For SU(2), all the soft matrices with an odd number of nonsinglets vanish.

Finally, when there are 4 nonsinglets, they are given as

$\displaystyle S^{(0)}(1,1,1,1)$
$\displaystyle=(N-2)\begin{pmatrix}\frac{2(N+2)}{N}(12)(34)&\frac{N+2}{2N}\big{(}(1243)-(2143)\big{)}&S_{13}\\\
-\frac{N+2}{2N}\big{(}(1234)-(2134)\big{)}&S_{22}&S_{23}\\\
S_{31}&S_{32}&S_{33}\end{pmatrix},$

where $S_{ij}$ are given by

$\displaystyle S_{13}$
$\displaystyle=\frac{N^{2}-12}{2N^{2}}\frac{N\big{(}(1243)+(2143)\big{)}-2(12)(34)}{N-2},$
$\displaystyle S_{22}$
$\displaystyle=\frac{N}{4(N-2)}\Bigl{(}(14)(23)+(13)(24)-2(1324)-2(1423)-(1234)-(1243)-(1342)-(1432)\Bigr{)},$
$\displaystyle S_{23}$ $\displaystyle=\frac{1}{4N^{2}(N-2)}\
\Bigl{[}16\Bigl{(}(1342)-(1234)\Bigr{)}+N(N^{2}-8)\Bigl{(}(13)(24)-(14)(23)\Bigr{)}$
$\displaystyle+N^{2}\Bigl{(}(1234)-(1243)-(1342)+(1432)\Bigr{)}\Bigr{]},$
$\displaystyle S_{31}$
$\displaystyle=\frac{N^{2}-12}{2N^{2}(N-2)}\Bigl{[}N\Bigl{(}(1234)+(2134)\Bigr{)}-2(12)(34)\Bigr{]},$
$\displaystyle S_{32}$
$\displaystyle=\frac{1}{4N^{2}(N-2)}\Bigl{[}16\Bigl{(}(1243)-(1432)\Bigr{)}+N(N^{2}-8)\Bigl{(}(13)(24)-(14)(23)\Bigr{)}$
$\displaystyle+N^{2}\Bigl{(}(1234)-(1243)-(1342)+(1432)\Bigr{)}\Bigr{]},$
$\displaystyle S_{33}$
$\displaystyle=\frac{1}{4N^{3}(N-2)}\Bigl{[}-64(12)(34)+16N\Bigl{(}2(1234)+2(1243)+2(1342)+2(1432))$
$\displaystyle+(1324)+(1423)\Bigr{)}+N^{2}(N^{2}-16)\Bigl{(}(14)(23)+(13)(24)\Bigr{)}$
$\displaystyle-N^{3}\Bigl{(}(1234)+(1243)+(1342)+(1432)-2(1423)-2(1324)\Bigr{)}\Bigl{]}.$
(138)

###### Acknowledgements.

This work is supported by Basic Science Research Program through the National
Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant
No. NRF-2019R1F1A1060396). T. K. is supported by Basic Science Research
Program through the National Research Foundation of Korea (NRF) funded by the
Ministry of Science and ICT (Grants No. NRF-2021R1A2C1008906).

## References

  * (1) C. W. Bauer, S. Fleming, and M. E. Luke, Summing Sudakov logarithms in $B\to X_{s}\gamma$ in effective field theory, Phys. Rev. D63 (2000) 014006, [hep-ph/0005275].
  * (2) C. W. Bauer, S. Fleming, D. Pirjol, and I. W. Stewart, An Effective field theory for collinear and soft gluons: Heavy to light decays, Phys. Rev. D63 (2001) 114020, [hep-ph/0011336].
  * (3) C. W. Bauer and I. W. Stewart, Invariant operators in collinear effective theory, Phys. Lett. B 516 (2001) 134–142, [hep-ph/0107001].
  * (4) C. W. Bauer, D. Pirjol, and I. W. Stewart, Soft collinear factorization in effective field theory, Phys. Rev. D65 (2002) 054022, [hep-ph/0109045].
  * (5) D. d’Enterria, Physics at the FCC-ee, in 17th Lomonosov Conference on Elementary Particle Physics, 2, 2016.  arXiv:1602.05043.
  * (6) ILC Collaboration, G. Aarons et al., International Linear Collider Reference Design Report Volume 2: Physics at the ILC, arXiv:0709.1893.
  * (7) FCC Collaboration, A. Abada et al., FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2, Eur. Phys. J. ST 228 (2019), no. 2 261–623.
  * (8) CLICdp, CLIC Collaboration, T. K. Charles et al., The Compact Linear Collider (CLIC) - 2018 Summary Report, arXiv:1812.06018.
  * (9) M. Ciafaloni, P. Ciafaloni, and D. Comelli, Bloch-Nordsieck violating electroweak corrections to inclusive TeV scale hard processes, Phys. Rev. Lett. 84 (2000) 4810–4813, [hep-ph/0001142].
  * (10) M. Ciafaloni, P. Ciafaloni, and D. Comelli, Bloch-Nordsieck violation in spontaneously broken Abelian theories, Phys. Rev. Lett. 87 (2001) 211802, [hep-ph/0103315].
  * (11) P. Ciafaloni and D. Comelli, The Importance of weak bosons emission at LHC, JHEP 09 (2006) 055, [hep-ph/0604070].
  * (12) A. Manohar, B. Shotwell, C. Bauer, and S. Turczyk, Non-cancellation of electroweak logarithms in high-energy scattering, Phys. Lett. B 740 (2015) 179–187, [arXiv:1409.1918].
  * (13) L. N. Lipatov, The parton model and perturbation theory, Yad. Fiz. 20 (1974) 181–198.
  * (14) V. N. Gribov and L. N. Lipatov, Deep inelastic e p scattering in perturbation theory, Sov. J. Nucl. Phys. 15 (1972) 438–450.
  * (15) G. Altarelli and G. Parisi, Asymptotic Freedom in Parton Language, Nucl. Phys. B 126 (1977) 298–318.
  * (16) Y. L. Dokshitzer, Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics., Sov. Phys. JETP 46 (1977) 641–653.
  * (17) J.-y. Chiu, A. Jain, D. Neill, and I. Z. Rothstein, The Rapidity Renormalization Group, Phys. Rev. Lett. 108 (2012) 151601, [arXiv:1104.0881].
  * (18) A. V. Manohar and W. J. Waalewijn, Electroweak Logarithms in Inclusive Cross Sections, JHEP 08 (2018) 137, [arXiv:1802.08687].
  * (19) J. Chay, T. Ha, and T. Kwon, N-jettiness in electroweak high-energy processes, JHEP 01 (2022) 007, [arXiv:2105.06672].
  * (20) R. Kelley and M. D. Schwartz, 1-loop matching and NNLL resummation for all partonic 2 to 2 processes in QCD, Phys. Rev. D 83 (2011) 045022, [arXiv:1008.2759].
  * (21) I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn, N-Jettiness: An Inclusive Event Shape to Veto Jets, Phys. Rev. Lett. 105 (2010) 092002, [arXiv:1004.2489].
  * (22) T. T. Jouttenus, I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn, The Soft Function for Exclusive N-Jet Production at Hadron Colliders, Phys. Rev. D 83 (2011) 114030, [arXiv:1102.4344].
  * (23) C. W. Bauer, A. Hornig, and F. J. Tackmann, Factorization for generic jet production, Phys. Rev. D 79 (2009) 114013, [arXiv:0808.2191].
  * (24) J. Collins, Foundations of Perturbative QCD.  Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology. Cambridge University Press, 2011.
  * (25) A. Idilbi and T. Mehen, On the equivalence of soft and zero-bin subtractions, Phys. Rev. D 75 (2007) 114017, [hep-ph/0702022].
  * (26) A. Idilbi and T. Mehen, Demonstration of the equivalence of soft and zero-bin subtractions, Phys. Rev. D 76 (2007) 094015, [arXiv:0707.1101].
  * (27) T. Becher and G. Bell, Analytic Regularization in Soft-Collinear Effective Theory, Phys. Lett. B 713 (2012) 41–46, [arXiv:1112.3907].
  * (28) J.-Y. Chiu, A. Jain, D. Neill, and I. Z. Rothstein, A Formalism for the Systematic Treatment of Rapidity Logarithms in Quantum Field Theory, JHEP 05 (2012) 084, [arXiv:1202.0814].
  * (29) Y. Li, D. Neill, and H. X. Zhu, An exponential regulator for rapidity divergences, Nucl. Phys. B 960 (2020) 115193, [arXiv:1604.00392].
  * (30) M. A. Ebert, I. Moult, I. W. Stewart, F. J. Tackmann, G. Vita, and H. X. Zhu, Subleading power rapidity divergences and power corrections for qT, JHEP 04 (2019) 123, [arXiv:1812.08189].
  * (31) J. Chay and C. Kim, Consistent treatment of rapidity divergence in soft-collinear effective theory, JHEP 03 (2021) 300, [arXiv:2008.00617].
  * (32) C. F. Berger, C. Marcantonini, I. W. Stewart, F. J. Tackmann, and W. J. Waalewijn, Higgs Production with a Central Jet Veto at NNLL+NNLO, JHEP 04 (2011) 092, [arXiv:1012.4480].
  * (33) G. P. Korchemsky and A. V. Radyushkin, Renormalization of the Wilson Loops Beyond the Leading Order, Nucl. Phys. B 283 (1987) 342–364.
  * (34) I. A. Korchemskaya and G. P. Korchemsky, On lightlike Wilson loops, Phys. Lett. B 287 (1992) 169–175.
  * (35) J.-y. Chiu, A. Fuhrer, R. Kelley, and A. V. Manohar, Factorization Structure of Gauge Theory Amplitudes and Application to Hard Scattering Processes at the LHC, Phys. Rev. D 80 (2009) 094013, [arXiv:0909.0012].
  * (36) M. Baumgart, O. Erdoğan, I. Z. Rothstein, and V. Vaidya, Breakdown of the naive parton model in super-weak scale collisions, Phys. Rev. D 100 (2019), no. 9 096008, [arXiv:1811.04120].