aainstitutetext: Max-Planck-Institut für Physik, Werner-Heisenberg-Institut,
D-80805 München, Germanybbinstitutetext: Dipartimento di Fisica and Arnold-
Regge Center, Università di Torino, and INFN, Sezione di Torino, Via P. Giuria
1, I-10125 Torino, Italyccinstitutetext: INPAC, Shanghai Key Laboratory for
Particle Physics and Cosmology, School of Physics and Astronomy, Shanghai Jiao
Tong University, Shanghai 200240, Chinaddinstitutetext: Key Laboratory for
Particle Astrophysics and Cosmology (MOE), Shanghai 200240, China

# A new method for calculating the soft anomalous dimension matrix for massive
particle scattering

Johannes Henn b Calum Milloy c,d Kai Yan<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The general structure of infrared divergences in the scattering of massive
particles is captured by the soft anomalous dimension matrix. The latter can
be computed from a correlation function of multiple Wilson lines. The state-
of-the-art two-loop result has a tantalizingly simple structure that is not
manifest in the calculations. We argue that the complexity in intermediate
steps of the known calculations comes from spurious, regulator-dependent
terms. Based on this insight we propose a different infrared regulator that is
associated to only one of the Wilson lines. We demonstrate that this
streamlines obtaining the two-loop result: computing the required Feynman
integrals via the differential equations method, only multiple polylogarithmic
functions appear (to all orders in the dimensional regulator), as opposed to
elliptic polylogarithms. We show that the new method is promising for higher-
loop applications by computing a three-loop diagram of genuine complexity, and
provide the answer in terms of multiple polylogarithms. The relatively simple
symbol alphabet we obtain may be of interest for bootstrap approaches.

††preprint: MPP-2023-242

##  1 Introduction

Being able to predict the structure of infrared divergences in scattering
processes is relevant for collider physics. For massless particles, the
general structure of divergences is known to three loops. The infrared
divergences factorize, similar to ultraviolet divergences. Due to color
dependence, the infrared factorization takes the form of a matrix in color
space. The essential part of the matrix is computed by the soft anomalous
dimension matrix, which is determined from a correlation function of Wilson
lines meeting at the origin, and pointing in the direction of the scattered
particles’ momenta Korchemskaya:1994qp .

Up to two loops the soft anomalous dimension matrix consists only of so-called
dipole terms, which depend on pairs of scattered particles only Dixon:2008gr ;
Becher:2009cu . Starting from three loops, correlations between multiple lines
(in particular, four lines) are also present. Interestingly, the four-line
correlations depend on a function of two scale-invariant ratios, which has
been computed at three loops Almelid:2015jia , and its presence has been
confirmed in four-particle scattering amplitudes Henn:2016jdu ; Caola:2021rqz
.

The structure of infrared divergences of massive particles is more intricate.
In the two-line case, the state of the art for the cusp anomalous dimension is
three-loop order in QCD Grozin:2014hna ; Grozin:2015kna  and four-loop order
in QED Bruser:2020bsh . Much less is known beyond the dipole terms. The two-
loop soft anomalous dimension matrix involving three massive particles was
computed by several groups Mitov:2009sv ; Ferroglia:2009ep ; Mitov:2010xw ;
Chien:2011wz . The result includes a color structure beyond the dipole terms,
namely $f^{abc}{\bf T}_{i}^{a}{\bf T}_{j}^{b}{\bf T}_{k}^{c}$, where
$\mathbf{T}_{i}^{a}$ is the colour matrix in the representation of particle
$i$ Catani:1996vz ; Bassetto:1983mvz . (At higher loop orders, further color
structure may appear.) The accompanying function depends on the scale
invariant combination of the scalar products between the velocities of the
particles,

$\displaystyle\beta_{ij}=\frac{v_{i}\cdot
v_{j}}{\sqrt{v_{i}^{2}}\sqrt{v_{j}^{2}}}\,.$ (1)

The function multiplying this color structure depends on three variables of
the type (1), since at two loops Feynman diagrams connect up to three Wilson
lines, cf. Fig. 2. Remarkably, the result is much simpler than a genuine
three-variable function: it takes the form of a sum over products of
logarithms and dilogarithms that individually depend on one variable only.

This simplicity is not transparent in the state-of-the-art calculations. It is
interesting to note that the calculations in the literature were done by
different means, via Mellin-Barnes methods Ferroglia:2009ep , by a position-
space calculation Mitov:2009sv ; Mitov:2010xw , and by a special gauge choice
Chien:2011wz . While all produce the correct result, it is not clear whether
the methods can be easily generalized to the next loop order. Moreover, for
the first two methods, the answer takes a simple form appears at the end of
the calculation only. This motivates us to understand why this is the case,
and to look for other approaches.

In recent years, the differential equation method Kotikov:1990kg ; Bern:1993kr
; Gehrmann:1999as  has proven extremely useful for the calculation of Feynman
integrals, especially thanks to new ideas about the transcendental structure
of Feynman integrals Henn:2013pwa . It is therefore interesting to see how the
two-loop calculation would look like with this method. It turns out
HennSimmonsDuffinUnbuplished ; Milloy:2020hzi  that the most complicated two-
loop Feynman diagram, shown in Fig. 2(b), involves elliptic functions. This
can be seen by a cut analysis, or by inspecting the differential equations
matrix. What this means is that the $\mathcal{O}(1/\epsilon)$ part of the
calculation, which contributes to the soft anomalous dimension matrix, is
simple, but the higher order terms in the $\epsilon$-expansion are much more
complicated. Although extracting the relevant information from the complicated
differential equations is possible Milloy:2020hzi , by carefully expanding
them in $\epsilon$, this is rather cumbersome. So, the question arises, is
there a simpler approach?

It is important to realize that the soft anomalous dimension matrix
corresponds to the leading divergence in $\epsilon$ of the Feynman integrals
(or, linear combinations of Feynman diagrams organized in so-called webs
Gatheral:1983cz ; Frenkel:1984pz ; Sterman:1981jc ). In the calculations, a
regularization choice is usually made to extract this divergence. In the
language of position-space Wilson line correlators, the divergence we are
interested in is a short-distance one, originating from webs of Feynman
diagrams that homogeneously scale to zero towards the origin. This divergence
is regulated by dimensional regularization, with $d=4-2\epsilon$, and
$\epsilon>0$. The soft anomalous dimension matrix corresponds to the
coefficient of the $1/\epsilon$ pole. However, the Wilson lines extend all the
way to infinity, and this introduces a long-distance divergence that we are
not interested in. It is common in the literature to regulate this by an
exponential suppression factor. It is important to realize that due to the
cutoff, only the leading $1/\epsilon$ pole is gauge-invariant and meaningful,
while e.g. the finite term is not. So, it is no surprise that the functions
contributing to the finite term can be very complicated - they do not have a
physical meaning, and depend on the regularization procedure!

This motivates us to search for a different way of regulating the long-
distance divergences of the webs. Our idea is the following: to make the webs
well-defined, it is sufficient to regulate one of the Wilson lines with an
exponential suppression factor. In this paper, we provide an argument why this
procedure leads to the correct soft anomalous dimension matrix. We introduce
the new method in section 2. We then show the effectiveness of the method in
section 3, by easily reproducing the known two-loop result, and by computing
the contribution of a previously unknown three-loop web function for the
three-lines soft anomalous dimension matrix, in section 4. We discuss
implications for prospects of future investigations in section 5.

##  2 Definitions and discussion of the method

###  2.1 Soft anomalous dimension matrix

Scattering amplitudes of massive particles contain infrared divergences due to
virtual exchange of small quanta between heavy particles. In particular, an
$n$-point gauge theory amplitude with massive external legs obeys a
factorization formula which characterizes its all-loop infrared structure,

$\displaystyle\mathcal{M}_{n}(p_{i},m_{i},\mu;\alpha_{s}(\mu);\epsilon)=S_{n}(v_{i};\alpha_{s}(\mu);\epsilon)H_{n}(p_{i},m_{i},\mu;\alpha_{s}(\mu))\,,$
(2)

where $S_{n}$ is the soft function encoding all the infrared singularities,
and $H_{n}$ collects the finite remainder of the amplitude ${\cal M}_{n}$. All
quantities are renormalized with ultraviolet (UV) poles removed. The running
coupling $\alpha_{s}(\mu)$ is evaluated at the scale $\mu$ in $\overline{\rm
MS}$ scheme. It satisfies the following equation in $d=4-2\epsilon$
dimensions,

$\displaystyle\frac{d\ln\alpha_{s}}{d\ln\mu}=-2\epsilon-2\sum_{k}\left(\frac{\alpha_{s}}{4\pi}\right)^{k}b_{k}\,.$
(3)

The $b_{k}$ coefficients are coefficients of the beta function
$\beta(\alpha_{s})$.

The $S_{n}$ is a universal function depending on the velocities $v_{i}\equiv
p_{i}/m_{i}$ of the scattering particles. It governs the singularity of the
amplitudes both in the infrared region and in the infinite mass limit of the
external particles, which can be studied in the context of Heavy Quark
Effective Theory (HQET) (see, e.g. manohar_wise_2000 ). In the limit
$m_{Q}\rightarrow\infty$, the heavy quark traveling with velocity $v_{i}$
behaves like a classical source, whose interactions with soft gluons are
represented by a Wilson line

$\displaystyle
Y_{i}(x)\equiv{\mathbf{P}}\,\text{exp}\left[ig\int_{0}^{\infty}ds\,v_{i}\cdot
A^{a}(x+v_{i}s)T^{a}\right]\,.$ (4)

Hence $S_{n}$ admits a formal definition as the correlation function of $n$
semi-infinite Wilson lines emanating from the origin, whose renormalization
group equation defines the corresponding soft anomalous dimension,

$\displaystyle
S_{n}=\langle{0}|\prod_{i=1}^{n}\,Y_{i}(0)|{0}\rangle_{\rm{ren.}},\quad\;\frac{d}{d\ln\mu}S_{n}=-S_{n}\,\mathbf{\Gamma}_{n}\,.$
(5)

The anomalous dimension $\mathbf{\Gamma}_{n}$ is a matrix operator in color
space acting on the multi-point amplitudes. It can be expanded as a sum over
dipole, tripole and quadrupole functions, etc,

$\displaystyle\mathbf{\Gamma}_{n}(v_{i};\alpha_{s})=\sum_{(i,j)}\mathbf{\Gamma}_{ij}(x_{ij};\alpha_{s})+\sum_{(i,j,k)}\mathbf{\Gamma}_{ijk}(x_{ij},x_{jk},x_{ik};\alpha_{s})+\cdots,$
(6)

which are parametrized by the variables

$\displaystyle v_{i}\cdot
v_{j}\equiv\cosh\phi_{ij}\equiv-\frac{1}{2}(x_{ij}+1/x_{ij})\,,$ (7)

where we assume for simplicity $v_{i}^{2}=1$. The operators
$\mathbf{\Gamma}_{ij\cdots}$ come from multi-Wilson-line diagrams describing
the interactions between Wilson lines labeled by $i,j,\ldots$, and hence
depend on the associated color charges, as well as on the angles between the
respective lines. The corresponding light-like soft anomalous dimension matrix
(which is reached via $i\phi_{ij}{\text{'}s}\rightarrow\infty$) is known up to
the three-loop order Almelid:2015jia . In this case, beyond-dipole terms
vanish at two loops and begin at three loops. This is in contrary to the
massive case, where beyond-dipole terms are already present at two loops.

The soft function obeys the exponentiation theorem for Wilson-line matrix
elements, whose logarithm defines a web function free from sub-divergences in
the infrared regime. More explicitly, it can be constructed through its
$\alpha_{s}-$expansion,

$\displaystyle W(\alpha_{s}(\mu))\equiv$ $\displaystyle\log\,\langle
Y_{1}\cdots Y_{n}\rangle_{\text{ren.}}$ (8) $\displaystyle=$
$\displaystyle\log\bigg{\\{}1+\sum_{k=1}^{\infty}\left(\frac{\alpha_{s}(\mu)}{4\pi}\right)^{k}S_{n}^{(k)}\bigg{\\}}$
(9) $\displaystyle=$
$\displaystyle\sum_{k=1}^{\infty}\left(\frac{\alpha_{s}(\mu)}{4\pi}\right)^{k}W^{(k)}\,.$
(10)

In the abelian theory, e.g. QED, the $L-$loop web function determines the
$L-$loop soft anomalous dimension through its evolution equation
$\frac{d}{d\ln\mu}W_{\rm QED}^{(L)}=-\Gamma^{(L)}_{\rm QED}$. In particular,
in the cases where all fermions are treated as massive, the abelian web
function is one-loop exact. In the non-abelian theory like QCD, higher-loop
web functions contain maximally non-abelian color structures, which, at a
given loop order, vanish in the abelian theory and cannot be decomposed onto
an abelian product of lower-loop color structures. More intuitively, the
$L-$loop color factors can be associated with web diagrams containing $L-1$
interation vertices among the virtual soft particles. These non-abelian
contributions give rise to multi-line color correlations in the web function.

We are interested in the tripole function $\mathbf{\Gamma}_{ijk}$ in the soft
anomalous dimension. Up to the three-loop order it has the following
structure,

$\displaystyle\mathbf{\Gamma}_{123}=\,$ $\displaystyle
if^{abc}\mathbf{T}_{1}^{a}\mathbf{T}_{2}^{b}\mathbf{T}_{3}^{c}\,F_{1}(x_{12},x_{13},x_{23};\alpha_{s})\,$
$\displaystyle+\Big{\\{}[\mathbf{T}_{1}^{c},\mathbf{T}_{1}^{a}]\,[\mathbf{T}_{2}^{b},\mathbf{T}_{2}^{c}]\,\\{\mathbf{T}_{3}^{a},\mathbf{T}_{3}^{b}\\}F_{2}(x_{12},x_{13},x_{23};\alpha_{s})\,+$
$\displaystyle\hskip
85.35826pt+(1,2,3)\rightarrow(2,3,1)+(1,2,3)\rightarrow(3,1,2)\Big{\\}}\,.$
(11)

Here the $\alpha_{s}$-expansion of $F_{1}$ starts at the two-loop order, i.e.

$\displaystyle
F_{1}=\sum_{k=2}^{\infty}\left(\frac{\alpha_{s}}{4\pi}\right)^{k}F_{1}^{(k)}\,,$
(12)

and that of $F_{2}$ starts at three-loop order, i.e.

$\displaystyle
F_{2}=\sum_{k=3}^{\infty}\left(\frac{\alpha_{s}}{4\pi}\right)^{k}F_{2}^{(k)}\,.$
(13)

In the light-like limit, the all-order $F_{1}$ vanishes due to Bose symmetry,
and $F_{2}^{(3)}$ reduces to a constant Almelid:2015jia . In this work we
introduce a new method to compute the full angle-dependent functions
$F_{1}^{(3)}$ and $F_{2}^{(3)}$.

###  2.2 Proposal of a new infrared regularization procedure

In dimensional regularization, Wilson-line correlators do not depend on a
physical scale. The soft function can be regularized by introducing
offshellness regulators $\delta_{i}$ on each Wilson line, which cut off the
infrared (IR) divergences,

$\displaystyle\mathcal{S}_{\delta_{1},\cdots,\delta_{n}}\equiv\langle
0|\prod_{i}Y_{i,\delta_{i}}|0\rangle,\quad
Y_{i,\delta_{i}}\equiv{\mathbf{P}}_{\delta_{i}}\,\text{exp}\left[ig\int_{0}^{\infty}ds\,v_{i}\cdot
A^{a}(sv_{i})T^{a}\right]\,,$ (14)

where the modified path-ordered exponential reads

$\displaystyle{\mathbf{P}}_{\delta}\,\text{exp}\left[\int ds\right]\equiv 1$
$\displaystyle+\int ds\,\theta(s)\,e^{i\delta s}+\frac{1}{2!}\iint
ds_{1}ds_{2}\,\theta(s_{1})\theta(s_{2}-s_{1})\,e^{i\delta s_{2}}$
$\displaystyle+\frac{1}{2!}\iint
ds_{1}ds_{2}\,\theta(s_{2})\theta(s_{1}-s_{2})\,e^{i\delta s_{1}}+\cdots\,.$
(15)

The modified soft function contains only UV $\epsilon-$poles, which can be
collected into a multiplicative renormalization $\mathcal{Z}-$factor for the
Wilson-line operator such that

$\displaystyle{\cal
S}_{\delta_{1},\dots,\delta_{n}}(v_{i},\delta_{i};\alpha_{s};\epsilon)={\cal
S}^{\text{fin.}}_{\delta_{1},\cdots,\delta_{n}}(v_{i},\delta_{i}/\delta_{n},\mu/\delta_{n};\alpha_{s}(\mu);\epsilon)\,\mathcal{Z}_{n}^{-1}(v_{i};\alpha_{s}(\mu);\epsilon)\,,$
(16)

where ${\cal S}^{\text{fin.}}$ is finite at $\epsilon=0$. The factor
$\mathcal{Z}_{n}$ is a matrix in color space satisfying the evolution equation

$\displaystyle\frac{d}{d\ln\mu}\mathcal{Z}_{n}=-\mathcal{Z}_{n}\,\mathbf{\Gamma}_{n}\,.$
(17)

In this way the soft anomalous dimension matrix can be determined order-by-
order in $\alpha_{s}$ by extracting the UV $\epsilon$-poles of the modified
corrector ${\cal S}_{\delta_{1},\cdots,\delta_{n}}$.

The matrix $\mathcal{Z}_{n}$ encodes the UV divergences and is therefore
independent of the choice of IR regulators. We have the freedom to set the
parameters $\delta_{i}$ to be any nonzero value. The standard choice is to
consider ${\cal S}_{1,\cdots,1}$ where all $\delta_{i}$’s are set to 1, see
e.g. Grozin:2015kna . In this work we will adopt a new approach by setting all
but one of them to be zero, namely to only introduce the regulator on the last
external leg. Hence we have the advantage to compute simpler HQET integrals.

In the following we will explain how to apply the new method to determine the
soft anomalous dimension matrix. For simplicity we focus on the three-line
case where we consider ${\cal S}_{0,0,1}$ in the $\overline{\rm MS}$ scheme.

Let us argue that ${\cal S}_{0,0,1}$ has the same $\epsilon-$pole structures
as the ratio between the standard three-line and two-line correlators, i.e.
$[{\cal S}_{1,1}]^{-1}{\cal S}_{1,1,1}$, which implies that the finite
remainder $\mathcal{S}_{0,0,1}^{\text{fin.}}(v_{i},\mu;\alpha_{s}(\mu))$ obeys

$\displaystyle\left[\mu\frac{\partial}{\partial\mu}+\beta(\alpha_{s})\frac{\partial}{\partial\alpha_{s}}\right]\mathcal{S}_{0,0,1}^{\text{fin.}}={\bf\Gamma_{2}}\,\mathcal{S}_{0,0,1}^{\text{fin.}}-\mathcal{S}_{0,0,1}^{\text{fin.}}\,{\bf{\Gamma_{3}}}\,.$
(18)

To see why the statement is true, it is helpful to introduce modified
correlators ${\cal S}_{\delta,\delta}$ and ${\cal S}_{\delta,\delta,1}$ where
$\delta\ll 1$. Let us consider the ratio $[{\cal
S}_{\delta,\delta}(v_{1},v_{2})]^{-1}S_{\delta,\delta,1}(v_{1},v_{2},v_{3})$.
Since the parameter $\delta$ is an IR cutoff, which characterizes the
virtuality of the heavy particle $p^{2}-m^{2}\simeq m\delta$, the infrared
singularity in ${\cal S}_{\delta,\delta,1}$, as $\delta$ goes to zero,
factorizes and cancels with the two-line correlator. Therefore the ratio
$[{\cal S}_{\delta,\delta}]^{-1}S_{\delta,\delta,1}$ is regular in the
$\delta\rightarrow 0$ limit at fixed order in $\alpha_{s}$. Formally we can
take the limit $\delta\rightarrow 0$ in this ratio, yielding

$\displaystyle\lim_{\delta\rightarrow
0}\;[\mathcal{S}_{\delta,\delta}]^{-1}\mathcal{S}_{\delta,\delta,1}=[\mathcal{S}_{0,0}]^{-1}\mathcal{S}_{0,0,1}=\mathcal{S}_{0,0,1}\,,$
(19)

where we set the scaleless function ${\cal S}_{0,0}=1$. Likewise,

$\displaystyle\lim_{\delta\rightarrow
0}\;[S^{\rm{fin.}}_{\delta,\delta}]^{-1}S_{\delta,\delta,1}^{\rm{fin.}}=\mathcal{S}_{0,0,1}^{\rm{fin.}}\,.$
(20)

Meanwhile, given the universal pole structure eq. (16), we have

$\displaystyle[\mathcal{Z}_{2}]^{-1}[S_{\delta,\delta}]^{-1}{\cal
S}_{\delta,\delta,1}\mathcal{Z}_{3}=[S^{\rm{fin.}}_{\delta,\delta}]^{-1}S_{\delta,\delta,1}^{\rm{fin.}}\,.$
(21)

Taking the $\delta\rightarrow 0$ limit on both sides we find

$\displaystyle\mathcal{Z}_{2}^{-1}\mathcal{S}_{0,0,1}\mathcal{Z}_{3}={\cal
S}_{0,0,1}^{\rm fin.}\,,$ (22)

which is equivalent to the statement in eq. (18). Using $\Gamma_{2}$, we can
determine $\Gamma_{3}$ from $\mathcal{S}_{0,0,1}$ by demanding that the left-
hand side of (22) is finite.

The procedure works with any number of Wilson lines at any loop order, as we
demonstrate presently. Indeed, writing
$\Gamma_{n}=\sum_{i=1}^{\infty}\alpha_{s}^{i}\Gamma_{n}^{(i)}$, and solving
the renormalization group equation (17) gives

$\displaystyle\log\mathcal{Z}_{n}=$
$\displaystyle\,\alpha_{s}\frac{\Gamma_{n}^{(1)}}{2\epsilon}+\alpha_{s}^{2}\left(\frac{\Gamma_{n}^{(2)}}{4\epsilon}-\frac{b_{0}}{4\epsilon^{2}}\Gamma_{n}^{(1)}\right)$
(23)
$\displaystyle+\alpha_{s}^{3}\left(\frac{\Gamma_{n}^{(3)}}{6\epsilon}+\frac{1}{48\epsilon^{2}}\left[\Gamma_{n}^{(1)},\Gamma_{n}^{(2)}\right]-\frac{1}{6\epsilon^{2}}\left(b_{0}\Gamma_{n}^{(2)}+b_{1}\Gamma_{n}^{(1)}\right)+\frac{b_{0}^{2}}{6\epsilon^{3}}\Gamma_{n}^{(1)}\right)$
$\displaystyle+\alpha_{s}^{4}\bigg{(}\frac{\Gamma_{n}^{(4)}}{8\epsilon}+\frac{1}{48\epsilon^{2}}\left[\Gamma_{n}^{(1)},\Gamma_{n}^{(3)}\right]-\frac{b_{0}}{8\epsilon^{2}}\Gamma_{n}^{(3)}+\frac{\Gamma_{n}^{(2)}}{8\epsilon^{2}}\left(\frac{b_{0}^{2}}{\epsilon}-b_{1}\right)$
$\displaystyle\hskip
56.9055pt-\frac{\Gamma_{n}^{(1)}}{8\epsilon^{2}}\left(\frac{b_{0}^{3}}{\epsilon^{2}}-\frac{2b_{0}b_{1}}{\epsilon}+b_{2}\right)-\frac{b_{0}}{48\epsilon^{3}}\left[\Gamma_{n}^{(1)},\Gamma_{n}^{(2)}\right]\bigg{)}.$

Next, we write

$\displaystyle{\cal
S}_{0,\dots,0,\delta}=\exp\left(\sum_{i}\alpha_{s}^{i}w^{(i)}_{n}\right)=\exp\left(\sum_{i,j}\alpha_{s}^{i}\epsilon^{\,j}w^{(i,j)}_{n}\right)\,,$
(24)

where $w_{n}^{(i)}$ is the coefficient of $\alpha_{s}^{i}$ in the exponent and
$w_{n}^{(i,j)}$ is the coefficient of $\alpha_{s}^{i}\epsilon^{\,j}$. We have
also reinstated an arbitrary regulator $\delta$. We then impose
$\mathcal{Z}_{n-1}^{-1}{\cal S}_{n}\mathcal{Z}_{n}=\text{finite}$, at each
order in the coupling constant. At the first order this implies the relation

$\displaystyle\Gamma_{n}^{(1)}=$
$\displaystyle\,\Gamma_{n-1}^{(1)}-2w_{n}^{(1,-1)}\,.$ (25)

Hence we obtain $\Gamma_{n}^{(1)}$ from $\Gamma_{n-1}^{(1)}$, supplemented by
$w_{n}^{(1,-1)}$. Proceeding to the next orders, we have

$\displaystyle\Gamma_{n}^{(2)}=$
$\displaystyle\,\Gamma_{n-1}^{(2)}-4w_{n}^{(2,-1)}+2\left[\Gamma_{n-1}^{(1)},w_{n}^{(1,0)}\right]-2\left[w_{n}^{(1,-1)},w_{n}^{(1,0)}\right]\,,$
(26) $\displaystyle\Gamma_{n}^{(3)}=$
$\displaystyle\Gamma_{n-1}^{(3)}-6w^{(3,-1)}_{n}+\frac{3}{2}b_{0}\left[w^{(1,-1)}_{n},w^{(1,1)}_{n}\right]+3\left[w^{(1,0)}_{n},w^{(2,-1)}_{n}\right]$
$\displaystyle+3\left[w^{(2,0)}_{n},w^{(1,-1)}_{n}\right]+\left[w^{(1,0)}_{n},\left[w^{(1,-1)}_{n},w^{(1,0)}_{n}\right]\right]-\left[w^{(1,-1)}_{n},\left[w^{(1,-1)}_{n},w^{(1,1)}_{n}\right]\right]$
$\displaystyle+\frac{3}{2}\left[\Gamma_{n-1}^{(2)},w^{(1,0)}_{n}\right]-\frac{3}{2}b_{0}\left[\Gamma_{n-1}^{(1)},w^{(1,1)}_{n}\right]$
$\displaystyle+\frac{3}{4}\left[\Gamma_{n-1}^{(1)},\left[w^{(1,-1)}_{n},w^{(1,1)}_{n}\right]\right]+\frac{3}{4}\left[w^{(1,-1)}_{n},\left[\Gamma_{n-1}^{(1)},w^{(1,1)}_{n}\right]\right]+3\left[\Gamma_{n-1}^{(1)},w^{(2,0)}_{n}\right]$
$\displaystyle-\frac{3}{4}\left[\Gamma_{n-1}^{(1)},\left[\Gamma_{n-1}^{(1)},w^{(1,1)}_{n}\right]\right]+\frac{3}{2}\left[w^{(1,0)}_{n},\left[w^{(1,0)}_{n},\Gamma_{n-1}^{(1)}\right]\right]\,.$
(27)

Eqs. (25), (26) and (27) specify all the ingredients needed to determine the
soft anomalous dimension matrix up to three-loop order. Requiring higher-order
poles to cancel implies further relations. For instance we have

$\displaystyle
w^{(2,-2)}_{n}=\frac{1}{4}\left[\Gamma_{n-1}^{(1)},w^{(1,-1)}_{n}\right]\,.$
(28)

This can be used as a consistency check when calculating the two-loop web
function $w^{(2)}_{n}$.

We will demonstrate in the following section that eq. (26) reproduces the
correct two-loop result for ${\bf{\Gamma}}_{3}^{(2)}$.

##  3 The two-loop soft anomalous dimension matrix via the new method

###  3.1 One loop

The one-loop soft anomalous dimension can be found from the recurrence
equation in eq. (25). Solving eq. (25) we find

$\displaystyle\Gamma_{n}^{(1)}=-2\sum_{i=1}^{n}w^{(1,-1)}_{i},$ (29)

where $w^{(1,-1)}_{i}$ is the single pole of the soft function that involves
Wilson lines $1$ to $i$, in eq. (24). We recall that only the $n$-th Wilson
line is IR-regulated.

There are two contributions to $w^{(1)}_{n}$: firstly, a gluon exchange
between the regulated $n$-th Wilson line and an unregulated $i$-th line, with
$i=\\{0,\dots,n-1\\}$, and secondly, the self-energy (SE) graph involving the
$n$-th Wilson line. The former is shown in Figure 1(a), whereas the latter in
Figure 1(b). Other potential contributions, due to gluon exchanges involving
unregulated lines only, vanish as the integrals are scaleless. For example,
the self-energy on an unregulated line is shown in Figure 1(c). Note that this
involves a cancellation between the IR and the UV.

It is convenient to compute the one-loop exchanges in configuration space
where the Feynman-gauge gluon propagator is given by

$\displaystyle\int\frac{d^{d}k}{(2\pi)^{d}}\frac{e^{-ik\cdot
x}(-ig_{\mu\nu})}{k^{2}+i\varepsilon}=-\frac{\Gamma(1-\epsilon)}{4\pi^{2-\epsilon}}\frac{g_{\mu\nu}}{(-x^{2}+i\varepsilon)^{1-\epsilon}}\,,$
(30)

and where $x$ is the spacetime distance associated to the gluon propagator. In
the diagram shown in Figure 1(a), it is given by $x=sv_{n}-tv_{i}$. Here $s,t$
are line parameters that are integrated along the semi-infinite Wilson line
segments, i.e. $s,t\in[0,\infty]$. In Figure 1(b) we have $t\in[0,s]$ and
$s\in[0,\infty]$.

(a) One-loop gluon exchange. More text.

(b) Self-energy correction on the regulated line.

(c) Self-energy correction on the unregulated line.

Figure 1: The three one-loop diagrams.

We write the total contribution to $w_{n}^{(1)}$ as

$\displaystyle\alpha_{s}w_{n}^{(1)}=$
$\displaystyle\,g_{s}^{2}\mu^{2\epsilon}\frac{\Gamma(1-\epsilon)}{4\pi^{2-\epsilon}}\bigg{[}\sum_{i=1}^{n-1}\mathbf{T}_{i}\cdot\mathbf{T}_{n}\,v_{i}\cdot
v_{n}\int_{0}^{\infty}dsdt\frac{e^{i\delta
s}}{(-(sv_{n}-tv_{i})^{2}+i\varepsilon)^{1-\epsilon}}$
$\displaystyle+\frac{C_{n}}{2}\,\int_{0}^{\infty}dse^{i\delta
s}\int_{0}^{s}dt\frac{1}{(-(sv_{n}-tv_{n})^{2}+i\varepsilon)^{1-\epsilon}}\bigg{]}\,,$
(31)

where $C_{i}=\mathbf{T}_{i}^{a}\mathbf{T}_{i}^{a}$ is the quadratic Casimir
associated to line $i$. The integrals in eq. (31) can be performed to all
orders in $\epsilon$, resulting in a ${}_{2}F_{1}$ function. However, we will
only require the expansion in $\epsilon$. In order to display the results, it
is convenient to define the integral $f(x,\epsilon)$ and its expansion in
$\epsilon$ as follows

$\displaystyle f(x,\epsilon)=$
$\displaystyle\,-\frac{1}{2}\,\left(x+\frac{1}{x}\right)\int_{0}^{\infty}dt\,\frac{1}{(t+\frac{1}{x})^{1-\epsilon}(t+x)^{1-\epsilon}}=\frac{r(x)}{2}\sum_{n=0}^{\infty}\frac{\epsilon^{n}}{n!}M_{n00}(x),$
(32a) $\displaystyle r(x)=$ $\displaystyle\,\frac{1+x^{2}}{1-x^{2}},$ (32b)
$\displaystyle M_{000}(x)=$ $\displaystyle\,2\log(x),$ (32c) $\displaystyle
M_{100}(x)=$
$\displaystyle\,2\text{Li}_{2}(x^{2})+4\log(x)\log(1-x^{2})-2\log^{2}(x)-2\zeta_{2}.$
(32d)

The functions $M_{n00}$ are a subset of a larger class of functions $M_{klm}$
which form a basis of functions for multiple-gluon exchange webs
Falcioni:2014pka , where three- and four-gluon vertices are absent. We only
display the explicit results up to $M_{100}$ as this is all we require to two
loops. The one-loop web contribution in eq. (31) can be written as

$\displaystyle\alpha_{s}w_{n}^{(1)}=$
$\displaystyle\,-\frac{\alpha_{s}}{\pi}\left(\frac{\mu^{2}}{\delta^{2}}\right)^{\epsilon}\frac{\Gamma(1-\epsilon)}{\pi^{-\epsilon}}\bigg{[}\sum_{i=1}^{n-1}\mathbf{T}_{i}\cdot\mathbf{T}_{n}\,f(x_{in},\epsilon)\Gamma(2\epsilon)+\frac{C_{n}}{2}\Gamma(-1+2\epsilon)\bigg{]}$
(33) $\displaystyle=$
$\displaystyle\,-\frac{\alpha_{s}}{\pi}\left(\frac{\bar{\mu}^{2}}{\overline{\delta}^{2}}\right)^{\epsilon}e^{\epsilon\gamma_{E}}\Gamma(1-\epsilon)\bigg{[}\sum_{i=1}^{n-1}\mathbf{T}_{i}\cdot\mathbf{T}_{n}\,f(x_{in},\epsilon)\Gamma(2\epsilon)+\frac{C_{n}}{2}\Gamma(-1+2\epsilon)\bigg{]},$
(34)

where $\bar{\mu}^{2}=4\pi e^{-\gamma_{E}}\mu^{2}$ is the
$\overline{\text{MS}}$ scale, $\overline{\delta}=2\delta$ is a modified
regulator and we have used the convenient $x_{ij}$ variables defined in eq.
(7). The coefficients as an expansion in $\epsilon$ are then

$\displaystyle\alpha_{s}w_{n}^{(1,-1)}=$
$\displaystyle-\frac{\alpha_{s}}{4\pi}\left[\sum_{i=1}^{n-1}\mathbf{T}_{i}\cdot\mathbf{T}_{n}r(x_{in})M_{000}(x_{in})-C_{n}\right],$
(35) $\displaystyle\alpha_{s}w_{n}^{(1,0)}=$
$\displaystyle-\frac{\alpha_{s}}{4\pi}\left[\sum_{i=1}^{n-1}\mathbf{T}_{i}\cdot\mathbf{T}_{n}r(x_{in})M_{100}(x_{in})-2\,C_{n}\right]$
$\displaystyle-\frac{\alpha_{s}}{4\pi}\log\left(\frac{\bar{\mu}^{2}}{\overline{\delta}^{2}}\right)\left[\sum_{i=1}^{n-1}\mathbf{T}_{i}\cdot\mathbf{T}_{n}r(x_{in})M_{000}(x_{in})-C_{n}\right].$
(36)

Then using eq. (35) in eq. (29), we have the result for the one-loop anomalous
dimension

$\displaystyle\alpha_{s}\,\Gamma_{n}^{(1)}=\frac{\alpha_{s}}{2\pi}\sum_{i=1}^{n}\left[\sum_{j=1}^{i-1}\mathbf{T}_{i}\cdot\mathbf{T}_{j}\,r(x_{ij})M_{000}(x_{ij})-\,C_{i}\right].$
(37)

Eq. (37) is in agreement with the classic result of Korchemsky and Radyushkin
Korchemsky:1987wg . In our computation of $\Gamma_{n}^{(1)}$, we have had to
handle conceptual issues around extra infrared poles, as correlations between
unregulated Wilson lines vanish. These are then accounted for in eq. (29) by
effectively adding them back. We will see in the remaining sections that using
our regularisation scheme is combinatorically more complex than when all the
lines are regulated due to a lack of symmetry. However, the integrals become
simpler and can be computed using standard methods, as only one line is
regulated.

(a)

(b)

Figure 2:  Contributions to the two-loop soft anomalous dimension matrix.

###  3.2 Two loops

The two-loop contribution to the soft anomalous dimension is described by eq.
(26). In this section we compute it for the case of three Wilson lines,
focusing on structures connecting all three lines, cf. eq. (11). In this case,
we do not have the first term of eq. (26), $\Gamma_{2}^{(2)}$, and we do not
consider any two-particle colour structures. Thus we write

$\Gamma_{3}^{(2)}\bigg{\rvert}_{{\text{connected}}}=\left\\{-4w_{3}^{(2,-1)}+2\left[\Gamma_{2}^{(1)},w_{3}^{(1,0)}\right]-2\left[w_{3}^{(1,-1)},w_{3}^{(1,0)}\right]\right\\}\bigg{\rvert}_{\text{connected}}.$
(38)

At two loops, there are two types of graphs connecting all three lines. One of
them is the gluon exchange diagram shown in Figure 2(a), which we label as
$[1,2,1]$. (This notation keeps track of the number of gluon attachments to
each line.) The other graph is the three-gluon-vertex diagram shown in Figure
2(b). The calculation of these graphs is instructive. The calculation of the
gluon exchange diagram shows us how to deal with a subtlety related to
infrared divergences, while that of the three-gluon-vertex diagram shows us
the computational advantages of regulating only one Wilson line.

####  3.2.1 Gluon exchange diagram

For the case of the $[1,2,1]$, the contributions to $w_{3}^{(2,-1)}$ in eq.
(38) are shown in Figure 3. The diagrams correspond to each permutation of
gluon attachments to the Wilson lines. The thick Wilson line, $v_{3}$ is the
regulated one. Diagrams 3(f) and 3(e) can be obtained from diagrams 3(a) and
3(b), respectively, by interchanging $v_{1}$ and $v_{2}$. Therefore, we only
need to compute diagrams 3(a) - (d).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3:  Two-loop gluon exchange diagrams contributing to three-line color
structures.

Diagrams 3(a) and 3(b) have the kinematic factor

$F_{\ref{fig:121comp}(a)}-F_{\ref{fig:121comp}(b)}=(v_{1}\cdot
v_{2})(v_{2}\cdot v_{3})\int_{0}^{\infty}dudt_{1}dt_{2}ds\,e^{i\delta
s}\frac{\theta(t_{1}-t_{2})-\theta(t_{2}-t_{1})}{(-(sv_{3}-t_{1}v_{2})^{2})^{1-\epsilon}(-(uv_{1}-t_{2}v_{2})^{2})^{1-\epsilon}}\,,$
(39)

whereas diagrams 3(c) and 3(d) are

$F_{\ref{fig:121comp}(c)}-F_{\ref{fig:121comp}(d)}=(v_{1}\cdot
v_{3})(v_{2}\cdot v_{3})\int_{0}^{\infty}dudtds_{1}ds_{2}\frac{e^{i\delta
s_{1}}\theta(s_{1}-s_{2})-e^{i\delta
s_{2}}\theta(s_{2}-s_{1})}{(-(s_{1}v_{3}-uv_{1})^{2})^{1-\epsilon}(-(s_{2}v_{3}-tv_{2})^{2})^{1-\epsilon}},$
(40)

where the relative minus signs comes from non-Abelian exponentiation. It is
reminded that the exponential damping factor $e^{i\delta s}$ is only added to
the outer most attachment on $v_{3}$. The contribution to $w_{3}^{(2)}$ is

$w^{(2)}_{3,\,[1,2,1]}=\frac{\mathcal{T}_{123}}{2\pi^{2}}\left[(F_{\ref{fig:121comp}(a)}-F_{\ref{fig:121comp}(b)})-(F_{\ref{fig:121comp}(a)}-F_{\ref{fig:121comp}(b)})(v_{1}\leftrightarrow
v_{2})+F_{\ref{fig:121comp}(c)}-F_{\ref{fig:121comp}(d)}\right].$ (41)

Here we denoted the overall tripole colour factor by $\mathcal{T}_{123}\equiv
if^{abc}\mathbf{T}_{1}^{a}\mathbf{T}_{2}^{b}\mathbf{T}_{3}^{c}$. The
combination $F_{\ref{fig:121comp}(a)}-F_{\ref{fig:121comp}(b)}$ can be
computed to give

$\displaystyle F_{\ref{fig:121comp}(a)}-F_{\ref{fig:121comp}(b)}=$
$\displaystyle\,e^{-4\epsilon\gamma_{E}}(\delta^{2})^{-2\epsilon}\frac{r(x_{12})r(x_{23})}{16}\bigg{[}\frac{1}{\epsilon^{2}}M_{000}(x_{12})M_{000}(x_{23})$
(42) $\displaystyle\hskip
17.07182pt+\frac{1}{\epsilon}\left(M_{000}(x_{12})M_{100}(x_{23})+M_{100}(x_{12})M_{000}(x_{23})\right)\bigg{]},$

whereas

$F_{\ref{fig:121comp}(c)}-F_{\ref{fig:121comp}(d)}=0,,$ (43)

due to symmetry. We define a modified $M_{100}$ function to absorb regulator-
dependent terms as follows,

$M_{100}^{(\kappa)}(x)=M_{100}(x)+\log(e^{-\kappa\gamma_{E}}\delta^{-2\kappa})M_{000}(x).$
(44)

Using $w^{(1,-1)}$ from eq. (35), $w^{(1,0)}$ from eq. (3.1) and
$\Gamma_{2}^{(1)}$ from eq. (37), the three contributions to eq. (38) are then
written as

$\displaystyle-4\,w_{3,[1,2,1]}^{(2,-1)}=$
$\displaystyle\,\frac{\mathcal{T}_{123}}{8\pi^{2}}r(x_{12})\bigg{\\{}r(x_{13})\left[M_{100}(x_{12})M_{000}(x_{13})+M_{000}(x_{12})M_{100}^{(2)}(x_{13})\right]$
$\displaystyle\hskip
56.9055pt-r(x_{23})\left[M_{100}(x_{12})M_{000}(x_{23})+M_{000}(x_{12})M_{100}^{(2)}(x_{23})\right]\bigg{\\}},$
$\displaystyle 2\left[\Gamma_{2}^{(1)},w_{3}^{(1,0)}\right]=$
$\displaystyle\,\frac{\mathcal{T}_{123}}{4\pi^{2}}r(x_{12})M_{000}(x_{12})\left[r(x_{23})M_{100}^{(1)}(x_{23})-r(x_{13})M_{100}^{(1)}(x_{13})\right],$
$\displaystyle-2\left[w_{3}^{(1,-1)},w_{3}^{(1,0)}\right]=$
$\displaystyle\,\frac{\mathcal{T}_{123}}{8\pi^{2}}r(x_{13})r(x_{23})\left[M_{100}^{(1)}(x_{13})M_{000}(x_{23})-M_{000}(x_{13})M_{100}^{(1)}(x_{23})\right].$
(45)

The sum of the above gives

$\Gamma_{3,[1,2,1]}^{(2)}=\mathcal{T}_{123}\left(\frac{1}{8\pi^{2}}\right)\sum_{i,j,k}\varepsilon_{ijk}\,r(x_{ij})\,r(x_{ik})\,M_{100}(x_{ij})\,M_{000}(x_{ik}),$
(46)

where all the dependence on the regulator $\delta$ cancels, as expected.

####  3.2.2 The three-gluon-vertex diagram

To complete the calculation of the three-line soft anomalous dimension at two
loops we require the coefficient of the single pole of Figure 2(b). Using our
new method that places a regulator only on one Wilson line, we find that only
four basis integrals are required in the top sector, as opposed to six in the
traditional approach HennSimmonsDuffinUnbuplished ; Milloy:2020hzi . (See also
Waelkens:2017zwt  for a calculation using unitarity cuts.)

The contribution to the exponent of the soft function is written as

$\displaystyle w^{(2)}_{3,\,3gv}=$
$\displaystyle\,\frac{\mathcal{T}_{123}}{16\pi^{2}}\left(\bar{\mu}^{2}\right)^{2\epsilon}e^{2\epsilon\gamma_{E}}$
$\displaystyle\times\int\frac{d^{d}k_{1}}{i\pi^{d/2}}\frac{d^{d}k_{2}}{i\pi^{d/2}}\int
d^{d}k_{3}\frac{\delta^{(d)}(k_{1}+k_{2}+k_{3})\sum_{i,j,k}\epsilon_{ijk}\left(v_{i}\cdot
v_{j}\right)\left(k_{i}\cdot
v_{k}\right)}{k_{1}^{2}k_{2}^{2}(k_{1}+k_{2})^{2}(k_{1}\cdot v_{1})(k_{2}\cdot
v_{2})(k_{3}\cdot v_{3}-\delta)}.$ (47)

We will conveniently set $\delta=1$ as the dependence on $\delta$ can be
easily recovered by counting mass dimensions. We define the integral family as

$\displaystyle I_{a_{1}a_{2}a_{3}a_{4}a_{5}a_{6}a_{7}a_{8}a_{9}}^{[3gv]}$
$\displaystyle=e^{2\epsilon\gamma_{E}}\int\frac{d^{d}k_{1}}{i\pi^{d/2}}\frac{d^{d}k_{2}}{i\pi^{d/2}}$
(48) $\displaystyle\times\frac{(k_{1}\cdot v_{2})^{-a_{7}}(k_{2}\cdot
v_{3})^{-a_{8}}(k_{2}\cdot
v_{1})^{-a_{9}}}{(k_{1}^{2})^{a_{1}}(k_{2}^{2})^{a_{2}}((k_{1}+k_{2})^{2})^{a_{3}}(k_{1}\cdot
v_{1})^{a_{4}}(k_{2}\cdot v_{2})^{a_{5}}(-(k_{1}+k_{2})\cdot
v_{3}-1)^{a_{6}}}.$

Using Kira Klappert:2020nbg  we identify 16 basis integrals, 4 of which are in
the top sector $I_{111111000}^{[3gv]}$. Differential equations are derived
with respect to the three $x_{ij}$ variables. The rotation to canonical form
Henn:2013pwa  is found using a combination of techniques. For the lower
sectors, using the automated package CANONICA Meyer:2017joq , is enough. For
the top sector we use DlogBasis Henn:2020lye  to reveal uniform weight
integrals. In this way we find the following three integrals,

$\displaystyle g_{13}^{[3gv]}$
$\displaystyle=\epsilon^{4}\frac{1-x_{13}^{2}}{x_{13}}\,\,I_{111111-100}^{[3gv]}\,,$
(49a) $\displaystyle g_{14}^{[3gv]}$
$\displaystyle=\epsilon^{4}\frac{1-x_{12}^{2}}{x_{12}}\,\,I_{1111110-10}^{[3gv]}\,,$
(49b) $\displaystyle g_{15}^{[3gv]}$
$\displaystyle=\epsilon^{4}\frac{1-x_{23}^{2}}{x_{23}}\,\,I_{11111100-1}^{[3gv]}\,.$
(49c)

The last required integral is found using INITIAL Dlapa:2020cwj

$\displaystyle g_{16}^{[3gv]}$
$\displaystyle=\epsilon^{3}\,\left(2I_{111112-10-1}^{[3gv]}+(v_{1}\cdot
v_{3})\,I_{111112-100}^{[3gv]}+(v_{2}\cdot
v_{3})\,I_{11111200-1}^{[3gv]}\right)+\text{lower sectors}.$ (50)

Using this integral basis, we find a differential equation in the canonical
form Henn:2013pwa

$\text{d}\,\mathbf{g}^{[3gv]}=\epsilon\sum_{\ell\in\mathcal{A}^{[3gv]}}c_{\ell}^{[3gv]}\,\text{d}\log(\ell)\,\mathbf{g}^{[3gv]}\,,$
(51)

where the set of $\ell$, the alphabet $\mathcal{A}$, is

$\displaystyle\mathcal{A}^{[3gv]}$
$\displaystyle=\\{1-x_{12},\,x_{12},\,1+x_{12},\,1-x_{13},\,x_{13},\,1+x_{13},\,1-x_{23},\,x_{23},\,1+x_{23},$
(52) $\displaystyle\hskip
71.13188ptx_{12}x_{13}+x_{23},\,x_{12}x_{23}+x_{13},\,x_{12}+x_{13}x_{23},\,1+x_{12}x_{13}x_{23}\\}\,.$

To solve the system in eq. (51) we need a boundary condition. We choose
$x_{ij}=1$ for all $i$ and $j$. At this boundary point, the lower sector
integrals are just iterated eikonal bubbles which evaluate to

$\displaystyle\int\frac{d^{d}k}{i\pi^{d/2}}\frac{1}{(k^{2})^{\alpha}(k\cdot
v-1)^{\beta}}$
$\displaystyle=(-1)^{\alpha+\beta}\frac{2^{d-2\alpha}\Gamma(d/2-\alpha)\Gamma(2\alpha+\beta-d)}{\Gamma(\alpha)\Gamma(\beta)},$
(53)

or, because of the differing regulators on different lines,

$\displaystyle\hskip-1.42262pt\int\frac{d^{d}k}{i\pi^{d/2}}$
$\displaystyle\frac{1}{(k^{2})^{\alpha}(k\cdot v-1)^{\beta}(k\cdot
v)^{\gamma}}$ (54)
$\displaystyle=(-1)^{\alpha+\beta+\gamma}\frac{2^{d-2\alpha}\Gamma(d/2-\alpha)\Gamma(d-2\alpha-\gamma)\Gamma(2\alpha+\beta+\gamma-d)}{\Gamma(\alpha)\Gamma(\beta)\Gamma(d-2\alpha)}.$

The top sector integrals $g_{13}^{[3gv]}$ to $g_{15}^{[3gv]}$ in eqs.
(49a)-(49c) clearly vanish at such a boundary. However, $g_{16}^{[3gv]}$ in
eq. (50) does not and depends on integrals more complicated than the bubble
integrals in eqs. (53) and (54). Instead of attempting to compute them or
determining from physical consistency conditions obtained from the
differential equations (see e.g. Henn:2020lye ), we will fix the remaining
boundary value from knowledge of the lightlike limit.

Solving the differential equation gives the following expansion in $\epsilon$
for $w_{3,\,3gv}^{(2)}$ in eq. (3.2.2)

$\displaystyle w^{(2,-3)}_{3,\,3gv}$
$\displaystyle=\mathcal{T}_{123}\frac{c_{0}}{16\pi^{2}}\left(r(x_{13})M_{000}(x_{13})-r(x_{23})M_{000}(x_{23})\right)$
(55a) $\displaystyle w^{(2,-2)}_{3,\,3gv}$
$\displaystyle=\mathcal{T}_{123}\frac{c_{1}}{16\pi^{2}}\left(r(x_{13})M_{000}(x_{13})-r(x_{23})M_{000}(x_{23})\right)$
(55b) $\displaystyle w^{(2,-1)}_{3,\,3gv}$
$\displaystyle=\frac{\mathcal{T}_{123}}{32\pi^{2}}\bigg{\\{}\frac{1}{2}\sum_{ijk}\epsilon_{ijk}r(x_{ij})M_{000}(x_{ij})M_{000}(x_{ik})^{2}$
(55c) $\displaystyle\hskip 85.35826pt-
c_{2}\left[r(x_{13})M_{000}(x_{13})-r(x_{23})M_{000}(x_{23})\right]\bigg{\\}}\,,$

where the $c_{i}$ are constants coming from the remaining integral and the
subscript denotes its transcendental weight. The web only has a single pole
divergence, thus $c_{0}=c_{1}=0$. Together with eq. (46) the complete three-
line two-loop soft anomalous dimension is

$\displaystyle\Gamma_{3}^{(2)}=-\mathcal{T}_{123}\frac{1}{8\pi^{2}}$
$\displaystyle\Bigg{\\{}\sum_{ijk}\epsilon_{ijk}r(x_{ij})M_{000}(x_{ij})\left(r(x_{ik})M_{100}(x_{ik})+\frac{1}{2}M_{000}(x_{ik})^{2}\right)$
$\displaystyle\hskip 48.36958pt-
c_{2}\bigg{(}r(x_{13})M_{000}(x_{13})-r(x_{23})M_{000}(x_{23})\bigg{)}\Bigg{\\}}.$
(56)

In the lightlike limit, when all $x_{ij}\to 0$ simultaneously, we have

$\displaystyle\lim_{x_{ij}\to
0}\Gamma_{3}^{(2)}=\mathcal{T}_{123}\frac{c_{2}}{4\pi^{2}}\log\frac{x_{13}}{x_{23}}.$
(57)

It is known that there is no tripole colour structure when all Wilson lines
are lightlike, since we cannot construct conformally invariant cross ratios
Dixon:2009ur . Therefore we conclude that $c_{2}=0$, giving

$\displaystyle\Gamma_{3}^{(2)}=-\mathcal{T}_{123}\frac{1}{8\pi^{2}}\sum_{ijk}\epsilon_{ijk}r(x_{ij})M_{000}(x_{ij})\left(r(x_{ik})M_{100}(x_{ik})+\frac{1}{2}M_{000}(x_{ik})^{2}\right),$
(58)

in agreement with previous computations Mitov:2009sv ; Ferroglia:2009ep ;
Mitov:2010xw ; Chien:2011wz .

In summary, we see that the calculation of the two-loop soft anomalous
dimension matrix is greatly simplified by regulating only one of the multiple
Wilson lines. In the next section, we show that the new method is promising
for higher-loop applications as well.

##  4 Three-loop calculation of a three-line web function

Figure 4: The three-loop $w_{122}^{(3)}$ diagram

To illustrate the potential of the new method, we evaluate in this section a
three-loop web diagram. This is a contribution to the soft anomalous dimension
matrix of three Wilson lines at the three-loop order. Naturally, many further
diagrams are required for the full calculation. Here, our intention is to
provide a proof-of-concept by calculating a diagram, shown in Fig. 4, of
genuine complexity.

The kinematic part of the integral shown in Fig. 4 can be conveniently written
as

$\displaystyle\mathcal{F}_{122}^{(3)}=$
$\displaystyle\int\frac{d^{d}k_{1}}{i\pi^{d/2}}\frac{d^{d}k_{2}}{i\pi^{d/2}}\frac{d^{d}k_{3}}{i\pi^{d/2}}$
(59) $\displaystyle\times$ $\displaystyle\int
d^{d}k_{4}\frac{\delta^{(d)}(k_{1}+k_{2}+k_{3})\sum_{i,j,k}\epsilon_{ijk}(v_{i}\cdot
v_{j})(v_{k}\cdot k_{i})}{k_{1}^{2}k_{2}^{2}k_{3}^{2}k_{4}^{2}(k_{1}\cdot
v_{1})((k_{1}-k_{4})\cdot v_{1})(k_{4}\cdot v_{2})((k_{2}+k_{4})\cdot
v_{2})(k_{3}\cdot v_{3}-1)}.$

Using Kira Klappert:2020nbg  we find that there are 66 basis integrals,
belonging to 22 sectors. The dimension of the sector with the most number of
integrals is 8. We then find the differential equation satisfied by these
integrals.

For the transformation of the differential equation to canonical form, we
found it convenient to transform sector-by-sector. The diagonal block
transformations were found using either CANONICA Meyer:2017joq  or using
DlogBasis Henn:2020lye  to identify one uniform weight integral and then using
INITIAL Dlapa:2020cwj  for the transformation.

Once all the diagonal blocks are in canonical form, we then need to transform
all the off-diagonal blocks. For this we used a combination of CANONICA and
Libra Lee:2020zfb . We write our equation as

$\text{d}\,\mathbf{g}^{[122]}=\epsilon\sum_{\ell\in\mathcal{A}^{[122]}}c_{\ell}^{[122]}\,\text{d}\log(\ell)\,\mathbf{g}^{[122]}\,,$
(60)

where $c_{\ell}^{[122]}$ are $66\times 66$ constant matrices. The canonical
differential equation matrices $c^{[122]}$, as well as the definition of the
canonical basis $\mathbf{g}^{[122]}$, are included in the ancillary files. The
alphabet appearing in the differential equation involves two additional
letters compared to the two-loop case of eq. (52), namely

$\displaystyle\\{x_{13}+x_{12}x_{23}+x_{12}x_{13}^{2}x_{23}+x_{13}x_{23}^{2},x_{12}x_{13}+x_{23}+x_{13}^{2}x_{23}+x_{12}x_{13}x_{23}^{2}\\}\,.$
(61)

In order to fully solve the differential equation we need a boundary
condition. In practice, we consider kinematic points
$v_{i}=v,v_{j}=v,v_{k}=-v$, where 36 of the 66 integrals can be found
analytically. Either they vanish or reduce to real-valued iterated integrals
of the form given in eqs. (53) and (54). This does not completely fix the
required constants for the single pole of $\mathcal{F}_{122}^{(3)}$ which we
label as $\mathcal{F}_{122}^{(3,-1)}$. We can find the remaining constants by
imposing finiteness of the differential equation on single-variable slices.
Consider the single-variable regime $(x_{12},x_{13},x_{23})=(1,x,x)$. Then we
impose the condition that eq. (60) remains finite as $x\to 1$ and $x\to\pm i$.
Doing this also for the other permutations we find we can fix all the boundary
conditions.

Thus we arrive at the result

$\displaystyle\mathcal{F}_{122}^{(3,-1)}=r(x_{12})^{2}$ $\displaystyle
f_{A}(x_{12},x_{13},x_{23})+r(x_{12})r(x_{13})f_{B}(x_{12},x_{13},x_{23})$
$\displaystyle+r(x_{12})r(x_{23})f_{C}(x_{12},x_{13},x_{23}),$ (62)

where the functions $f_{A}$, $f_{B}$ and $f_{C}$ are of pure weight five. They
can be written in terms of classical polylogarithms $\text{Li}_{n}(-s)$ up to
weight $n=5$, with the following arguments,

$\displaystyle
s\in\\{x_{12},-x_{12},x_{13},-x_{13},\frac{x_{12}}{x_{13}x_{23}},\frac{x_{12}x_{13}}{x_{23}},\frac{x_{12}x_{23}}{x_{13}},x_{12}x_{13}x_{23}\\}.$
(63)

Their symbols are drawn from a smaller set of alphabets given in eq. (52). The
expressions for $f_{A,B,C}$ are given explicitly in the ancillary files.

We would like to present the total contribution of web function
$w^{(3)}_{122}$ after symmetrization under permutations of Wilson lines. Hence
let us introduce a basis $\\{g_{A}^{\pm},g_{B}^{\pm},g^{\pm}_{C}\\}$ defined
as follows

$\displaystyle g_{A}^{\pm}(3,2,1)$ $\displaystyle=\pm
g_{A}^{\pm}(3,1,2)\equiv\frac{1}{2}\left[f_{A}(x_{12},x_{13},x_{23})\pm(x_{13}\leftrightarrow
x_{23})\right],$ $\displaystyle g_{B}^{\pm}(3,2,1)$ $\displaystyle=\pm
g_{B}^{\pm}(2,3,1)\equiv\frac{1}{2}\left[f_{B}(x_{12},x_{13},x_{23})\pm(x_{12}\leftrightarrow
x_{13})\right],$ $\displaystyle g_{C}^{\pm}(3,2,1)$ $\displaystyle=\pm
g_{C}^{\pm}(1,2,3)\equiv\frac{1}{2}\left[f_{C}(x_{12},x_{13},x_{23})\pm(x_{12}\leftrightarrow
x_{23})\right]\,.$ (64)

where we adopt a short-hand notation
$g(i,j,k)=g(\alpha_{i},\alpha_{j},\alpha_{k})$ where

$\displaystyle\alpha_{3}\equiv x_{12}\,,\quad\alpha_{2}\equiv
x_{13}\,,\quad\alpha_{1}\equiv x_{23}\,.$ (65)

In this notation the permutation among Wilson lines $v_{i,j,k}$ is equivalent
to the permutation among arguments of the function $g(i,j,k)$.

Summing over all six permutations of the web diagram ${\cal F}^{(3)}_{122}$,
we obtain

$\displaystyle\mathcal{C}(1,2,3)\,{\cal
F}_{122}^{(3)}(3,2,1)+{\rm{perm.}}=\frac{1}{\epsilon}\times$ (66)
$\displaystyle\Bigg{\\{}r(1)^{2}\,\left([\mathcal{C}(3,2,1)+\mathcal{C}(2,3,1)]\,g_{A}^{+}(1,2,3)+[\mathcal{C}(3,2,1)-\mathcal{C}(2,3,1)]\,g_{A}^{-}(1,2,3)\right)$
$\displaystyle+r(3)r(2)\,\left([\mathcal{C}(1,2,3)+\mathcal{C}(1,3,2)]\,g_{B}^{+}(3,2,1)+[\mathcal{C}(1,2,3)-\mathcal{C}(1,3,2)]\,g_{B}^{-}(3,2,1)\right.$
$\displaystyle\qquad\qquad\;\;+\left.[\mathcal{C}(2,1,3)+\mathcal{C}(3,1,2)]\,g_{C}^{+}(3,1,2)+[\mathcal{C}(2,1,3)-\mathcal{C}(3,1,2)]\,g_{C}^{-}(3,1,2)\right)$
$\displaystyle+{\rm cyc.}\Bigg{\\}}.$

The expression in the braces is manifestly symmetric under exchange
$2\leftrightarrow 3\,$, and cyc. represents the sum of two cyclic
permutations. The color factor
$\mathcal{C}(1,2,3)=\mathbf{T}_{1}^{c}\,[\mathbf{T}_{3}\cdot\mathbf{T}_{1},\mathbf{T}_{3}\cdot\mathbf{T}_{2}]\,\mathbf{T}_{2}^{c}$
is associated with the diagram we computed, whose maximally non-abelian
component can be decomposed onto the following basis,

$\displaystyle\mathbf{T}_{1}^{c}\,[\mathbf{T}_{3}\cdot\mathbf{T}_{1},\mathbf{T}_{3}\cdot\mathbf{T}_{2}]\,\mathbf{T}_{2}^{c}-\frac{1}{2}\Big{\\{}\mathbf{T}_{1}\cdot\mathbf{T}_{2},[\mathbf{T}_{3}\cdot\mathbf{T}_{1},\mathbf{T}_{3}\cdot\mathbf{T}_{2}]\Big{\\}}=$
(67)
$\displaystyle\frac{1}{4}\,[\mathbf{T}_{3}^{a},\mathbf{T}_{3}^{b}]\,[\mathbf{T}_{1}^{c},\mathbf{T}_{1}^{a}]\,\\{\mathbf{T}_{2}^{b},\mathbf{T}_{2}^{c}\\}+\frac{1}{4}\,[\mathbf{T}_{3}^{a},\mathbf{T}_{3}^{b}]\,\\{\mathbf{T}_{1}^{c},\mathbf{T}_{1}^{a}\\}\,[\mathbf{T}_{2}^{b},\mathbf{T}_{2}^{c}]+\frac{1}{2}\,[\mathbf{T}_{3}^{a},\mathbf{T}_{3}^{b}]\,[\mathbf{T}_{1}^{c},\mathbf{T}_{1}^{a}]\,[\mathbf{T}_{2}^{b},\mathbf{T}_{2}^{c}]\,.$

We would like to isolate the contributions to $F_{1}-$ and $F_{2}-$term in the
soft anomalous dimension. For this purpose it is convenient to regroup terms
associated with $r(3)r(2)$ in the following way

$\displaystyle
r(3)r(2)\,\big{(}[\mathcal{C}(1,2,3)+\mathcal{C}(2,1,3)]\,g_{BB}+[\mathcal{C}(1,3,2)+\mathcal{C}(3,1,2)]\,g_{CC}$
$\displaystyle\qquad\qquad+\left.[\mathcal{C}(1,2,3)-\mathcal{C}(2,1,3)+\mathcal{C}(1,3,2)-\mathcal{C}(3,1,2)]\,g_{+-}\right.$
$\displaystyle\qquad\qquad+[\mathcal{C}(1,2,3)-\mathcal{C}(2,1,3)-\mathcal{C}(1,3,2)+\mathcal{C}(3,1,2)]\,g_{--}\big{)},$
(68)

where

$\displaystyle
g_{BB}=g_{B}^{+}(3,2,1)+g_{C}^{+}(3,1,2)+g_{B}^{-}(3,2,1)+g_{C}^{-}(3,1,2),\quad
g_{CC}=g_{BB}|_{2\leftrightarrow 3}$ $\displaystyle
g_{+-}=g_{B}^{+}(3,2,1)-g_{C}^{+}(3,1,2),\quad
g_{--}=g_{B}^{-}(3,2,1)-g_{C}^{-}(3,1,2).$ (69)

Investigating the structure of color factors, we observe that only
$\\{g_{A}^{-},g_{--}\\}$ contribute to the $F_{1}-$term in the soft anomalous
dimension , whereas only $\\{g_{A}^{+},g_{BB},g_{CC}\\}$ contribute to the
$F_{2}-$term. The expressions for
$\\{g_{A}^{-}(1,2,3),g_{--},g_{A}^{+}(1,2,3),g_{BB}\\}$ are given in the
following, which involve products of logarithms with classical polylogarithms
up to weight four.

$\displaystyle g_{A}^{-}$ $\displaystyle=-\frac{8}{9}\,l_{+}l_{-}h_{3}(1),$
(70) $\displaystyle g_{--}$
$\displaystyle=-\frac{4}{9}\,\bigg{[}-6\,l_{-}{\cal L}_{4}^{c}+6\,l_{+}{\cal
L}_{4}^{d}+4\,l_{3}h^{b}_{4}(2)-4\,l_{2}h^{b}_{4}(3)$ (71)
$\displaystyle-3\,l_{-}l_{+}{\cal L}_{3}^{b}+6\,l_{1}l_{-}{\cal
L}_{3}^{c}-6\,l_{1}l_{+}{\cal L}_{3}^{d}$
$\displaystyle+3\,l_{1}l_{-}l_{+}{\cal
L}_{2}^{b}-l_{-}(3\,l_{1}^{2}+\pi^{2}){\cal
L}_{2}^{c}+l_{+}(3\,l_{1}^{2}+\pi^{2}){\cal L}_{2}^{d}$
$\displaystyle+l_{2}l_{3}l_{-}l_{+}{\cal
L}_{1}^{a}+\frac{1}{2}\,l_{-}l_{+}(-3\,l_{1}^{2}+l_{2}^{2}+l_{3}^{2}+\pi^{2}){\cal
L}_{1}^{b}+l_{1}l_{-}(l_{1}^{2}+\pi^{2}){\cal
L}_{1}^{c}-l_{1}l_{+}(l_{1}^{2}+\pi^{2}){\cal L}_{1}^{d}\bigg{]},$
$\displaystyle g_{A}^{+}$ $\displaystyle=\frac{4}{9}\,l_{1}\bigg{[}-6\,{\cal
L}_{4}^{a}-6\,l_{+}{\cal L}_{3}^{c}-6\,l_{-}{\cal L}_{3}^{d}-8\,h^{a}_{4}(1)$
(72) $\displaystyle-\,(-3\,l_{1}^{2}+3\,l_{2}^{2}+3\,l_{3}^{2}+\pi^{2}){\cal
L}_{2}^{a}-6\,l_{2}l_{3}{\cal L}_{2}^{b}$
$\displaystyle-2\,l_{1}(l_{1}^{2}+\pi^{2}){\cal
L}_{1}^{a}-l_{+}(-3\,l_{1}^{2}+l_{+}^{2}+\pi^{2}){\cal
L}_{1}^{c}-l_{-}(-3\,l_{1}^{2}+l_{-}^{2}+\pi^{2}){\cal L}_{1}^{d}\bigg{]},$
$\displaystyle g_{BB}$ $\displaystyle=-\frac{4}{9}\,l_{3}\bigg{[}-6\,({\cal
L}_{4}^{c}+{\cal L}_{4}^{d})-6\,l_{3}{\cal L}_{3}^{b}+6\,l_{1}({\cal
L}_{3}^{c}+{\cal L}_{3}^{d})+8\,h_{4}^{a}(2)$ (73)
$\displaystyle+6\,l_{1}l_{3}{\cal
L}_{2}^{b}+(-3\,l_{1}^{2}+3\,l_{+}l_{-}-\pi^{2})({\cal L}_{2}^{c}+{\cal
L}_{2}^{d})$ $\displaystyle+2\,l_{2}(l_{2}^{2}+\pi^{2}){\cal
L}_{1}^{a}-l_{3}(3\,l_{1}^{2}-3\,l_{2}^{2}+l_{3}^{2}+\pi^{2}){\cal
L}_{1}^{b}+l_{1}(l_{1}^{2}-3\,l_{+}l_{-}+\pi^{2})({\cal L}_{1}^{c}+{\cal
L}_{1}^{d})\bigg{]}.$

Here we define logarithms

$\displaystyle l_{1}$ $\displaystyle=\ln\alpha_{1},\quad
l_{2}=\ln\alpha_{2},\quad l_{3}=\ln\alpha_{3},\quad
l_{+}=\ln\alpha_{2}\alpha_{3},\quad l_{-}=\ln\frac{\alpha_{2}}{\alpha_{3}},$
(74)

and basis of polylogarithms

$\displaystyle{\cal L}_{n}^{a}$
$\displaystyle=\frac{1}{2}\left[\text{Li}_{n}\left(-\frac{\alpha_{1}}{\alpha_{2}\alpha_{3}}\right)+\text{Li}_{n}\left(-\frac{\alpha_{1}\alpha_{3}}{\alpha_{2}}\right)+\text{Li}_{n}\left(-\frac{\alpha_{1}\alpha_{2}}{\alpha_{3}}\right)+\text{Li}_{n}\left(-\alpha_{1}\alpha_{2}\alpha_{3}\right)\,\right]\;\pm\left(\alpha\leftrightarrow\frac{1}{\alpha}\right)$
$\displaystyle{\cal L}_{n}^{b}$
$\displaystyle=\frac{1}{2}\left[\text{Li}_{n}\left(-\frac{\alpha_{1}}{\alpha_{2}\alpha_{3}}\right)-\text{Li}_{n}\left(-\frac{\alpha_{1}\alpha_{3}}{\alpha_{2}}\right)-\text{Li}_{n}\left(-\frac{\alpha_{1}\alpha_{2}}{\alpha_{3}}\right)+\text{Li}_{n}\left(-\alpha_{1}\alpha_{2}\alpha_{3}\right)\,\right]\;\pm\left(\alpha\leftrightarrow\frac{1}{\alpha}\right)$
$\displaystyle{\cal L}_{n}^{c}$
$\displaystyle=\frac{1}{2}\left[\text{Li}_{n}\left(-\frac{\alpha_{1}}{\alpha_{2}\alpha_{3}}\right)-\text{Li}_{n}(-\alpha_{1}\alpha_{2}\alpha_{3})\,\right]\;\pm\left(\alpha\leftrightarrow\frac{1}{\alpha}\right)$
$\displaystyle{\cal L}_{n}^{d}$
$\displaystyle=\frac{1}{2}\left[\text{Li}_{n}\left(-\frac{\alpha_{1}\alpha_{3}}{\alpha_{2}}\right)-\text{Li}_{n}\left(-\frac{\alpha_{1}\alpha_{2}}{\alpha_{3}}\right)\,\right]\;\pm\left(\alpha\leftrightarrow\frac{1}{\alpha}\right)\,,$
(75)

where $\pm$ stands for $+$ if $n$ is odd, and $-$ if $n$ is even. By
construction, these basis are (anti-)symmetric under the transformation which
inverts all $\alpha$’s.

In addition, there are terms which depend on a single cusp angle. We group
them together and introduce the following set of functions

$\displaystyle h_{3}(i)$ $\displaystyle=-24\,li_{3}(i),$ $\displaystyle
h^{a}_{4}(i)$
$\displaystyle=6\,li_{4}(i)-4\,l_{i}\,li_{3}(i)+(l_{i}^{2}+\pi^{2})\,li_{2}(i)\,,\quad
h^{b}_{4}(i)=72\,li_{4}(i)+h_{4}^{a}(i),$ (76)

where

$\displaystyle
li_{n}(i)\equiv\frac{1}{2^{n}}\left[\text{Li}_{n}(\alpha_{i}^{2})-\zeta_{n}-\frac{2}{n-1}\ln\alpha_{i}\,(\text{Li}_{n-1}(\alpha_{i}^{2})-\zeta_{n-1})\right]\;\pm\left(\alpha_{i}\leftrightarrow\frac{1}{\alpha_{i}}\right)\,.$
(77)

Again, $\pm$ stands for $+$ if $n$ is odd, and $-$ if $n$ is even, $\forall
n\in\\{2,3,4\\}$. In the special case where $n=2$, the $\zeta_{n-1}$ term in
(77) should be set to 0.

Let us emphasize that the basis of polylogarithms given in (74),(4) and (77)
satisfy the following properties. Firstly, they are single-valued in the so-
called Euclidean region, where $\alpha_{1},\alpha_{2},\alpha_{3}>0$. In
particular the branch cuts on the positive real axis of $\alpha_{i}$ cancel in
(77). Secondly, each member in the set of function we introduced, i.e.

$\\{l_{1},l_{2},l_{3},{\cal L}_{n}^{a},{\cal L}_{n}^{b},{\cal
L}_{n}^{c}\pm{\cal L}_{n}^{d},li_{n}\\}$

is either symmetric or anti-symmetric under the transformation which inverts
any set of cusp angle variables
$(\alpha_{i}\rightarrow\alpha_{i}^{-1},\,\forall\,i)$. Moreover, the functions

$\\{l_{+},l_{-},{\cal L}_{n}^{a},{\cal L}_{n}^{b},{\cal L}_{n}^{c},{\cal
L}_{n}^{d}\\}$

are (anti-)symmetric under the exchange $(2\leftrightarrow 3)$. Given these
nice properties, the form of the final result is highly constrained.

##  5 Discussion and outlook

In this paper we proposed and tested a new regularization scheme for computing
the soft anomalous dimension matrix. We considered the ratio between the
correlation function for $n$ and $n-1$ Wilson-line operators, which is free
from IR divergences, provided that the $n-$th Wilson-line is dressed with an
offshellness regulator. Based on this observation, we defined a renormalized
quantity ${\cal S}_{0,\cdots,1}^{\text{fin.}}$ that satisfies the
renormalization group equation

$\displaystyle\left[\mu\frac{\partial}{\partial\mu}+\beta(\alpha)\frac{\partial}{\partial\alpha}\right]\mathcal{S}_{0,\cdots,1}^{\text{fin.}}={\bf\Gamma_{n-1}}\,\mathcal{S}_{0,\cdots,1}^{\text{fin.}}-\mathcal{S}_{0,\cdots,1}^{\text{fin.}}\,{\bf{\Gamma_{n}}}\,.$
(78)

This allows us to determine the soft anomalous dimension matrix for an
arbitrary number of massive particles, by sequentially computing
$\mathcal{S}_{0,1},\mathcal{S}_{0,0,1},\cdots,\mathcal{S}_{0,\cdots,1}$ and
extracting their UV poles.

We tested the novel method by demonstrating how the known one- and two-loop
results are obtained within this setup. Employing the differential equations
method, we note that calculations are simpler compared to the traditional
setup where all Wilson lines are regulated in the IR. As a proof of principle,
we evaluated one of the three-loop web diagrams that contributes to the three-
Wilson-line soft anomalous dimension matrix ${\bf\Gamma}_{3}$. A natural next
step, beyond the scope of this paper, is to evaluate the other relevant
Feynman diagrams contributing to this observable.

Our two-loop and sample three-loop calculations shed light on the new features
that might appear in ${\bf\Gamma}_{3}$. The final formula for the two-loop
soft anomalous dimension matrix depends on the following alphabet only,

$\displaystyle\\{1-x_{12},\,x_{12},\,1+x_{12},\,1-x_{13},\,x_{13},\,1+x_{13},\,1-x_{23},\,x_{23},\,1+x_{23}\\}\,.$
(79)

Although this depends on three variables, the dependence is factorized, which
implies that any function within this alphabet is a sum of products of single-
variable functions. In-contrast, the web diagrams we computed involve the
thirteen-letter alphabet (52). Compared to (79), this involves the following
new alphabet letters,

$\displaystyle\\{x_{12}+x_{13}x_{23},x_{13}+x_{12}x_{23},x_{23}+x_{12}x_{13},1+x_{12}x_{13}x_{23}\\}\,.$
(80)

This implies more complicated and richer structures for the transcendental
functions.111Also, the result of ref. Liu:2022elt  for the soft anomalous
dimension matrix for one massive and two massless particles may suggest the
appearance of further alphabet letters. (At the same time, it is interesting
that ${\cal F}_{122}$ is given by products of logarithms with classical
polylogarithms whose arguments are drawn from a small set of variables.) These
observations motivate a more systematic study on the function space of multi-
line soft anomalous dimension, which will be valuable input for bootstrap
approaches.

A further interesting direction is to investigate the web function in
$\mathcal{N}=4$ super Yang-Mills (sYM). At the three-loop order, the
difference between the corresponding soft anomalous dimension matrices is
captures by simpler, matter-dependent terms. Interestingly, those terms have
been observed to exhibit an universal iterative structure Grozin:2015kna  in
the two-line case, namely

$\displaystyle\Gamma_{\text{cusp}}(\phi)=C_{R}\frac{a}{\pi}\left[\Omega(\phi)+\frac{a}{\pi}C_{A}\Omega_{A}(\phi)+\left(\frac{a}{\pi}\right)^{2}C_{A}^{2}\Omega_{AA}(\phi)\right]+{\cal
O}(\alpha_{s}^{4})\,.$ (81)

The angle-dependent function $\Omega$’s are independent of the particle
content of the gauge theory, whereas $n_{f}-$ and $n_{s}-$dependence can be
associated with an effective coupling constant $a$. (This holds up to three
loops; at four loops terms proportional to the quartic Casimir color component
break this pattern Bruser:2020bsh .) It would be interesting to see whether
such a pattern exists for the tripole function $F_{1}^{(2)}$ and $F_{1}^{(3)}$
as well. Finally, in addition to considering the standard Wilson line
operator, it is possible to consider the (locally supersymmetric) Maldacena
operator, which includes a term coupling to scalars. This leads to additional
sets of angle parameters in the flavour space of scalars that the correlation
functions depend on, and that can be used to organize calculations and define
novel limits Correa:2012nk .

##  6 Acknowledgements

We would like to thank Einan Gardi for enlightening discussions. This research
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (grant agreement No
725110), Novel structures in scattering amplitudes. KY is supported by the
National Natural Science Foundation of China under Grant No. 12357077 and
would like to thank the sponsorship from Yangyang Development Fund. CM has
been partially supported by the Italian Ministry of University and Research
(MUR) through grant PRIN20172LNEEZ and by Compagnia di San Paolo through grant
TORPS1921EX-POST2101.

## References

  * (1) I. A. Korchemskaya and G. P. Korchemsky, _High-energy scattering in QCD and cross singularities of Wilson loops_ , _Nucl. Phys. B_ 437 (1995) 127 [hep-ph/9409446].
  * (2) L. J. Dixon, L. Magnea and G. F. Sterman, _Universal structure of subleading infrared poles in gauge theory amplitudes_ , _JHEP_ 08 (2008) 022 [0805.3515].
  * (3) T. Becher and M. Neubert, _Infrared singularities of scattering amplitudes in perturbative QCD_ , _Phys. Rev. Lett._ 102 (2009) 162001 [0901.0722].
  * (4) y. Almelid, C. Duhr and E. Gardi, _Three-loop corrections to the soft anomalous dimension in multileg scattering_ , _Phys. Rev. Lett._ 117 (2016) 172002 [1507.00047].
  * (5) J. M. Henn and B. Mistlberger, _Four-Gluon Scattering at Three Loops, Infrared Structure, and the Regge Limit_ , _Phys. Rev. Lett._ 117 (2016) 171601 [1608.00850].
  * (6) F. Caola, A. Chakraborty, G. Gambuti, A. von Manteuffel and L. Tancredi, _Three-loop helicity amplitudes for four-quark scattering in massless QCD_ , _JHEP_ 10 (2021) 206 [2108.00055].
  * (7) A. Grozin, J. M. Henn, G. P. Korchemsky and P. Marquard, _Three Loop Cusp Anomalous Dimension in QCD_ , _Phys. Rev. Lett._ 114 (2015) 062006 [1409.0023].
  * (8) A. Grozin, J. M. Henn, G. P. Korchemsky and P. Marquard, _The three-loop cusp anomalous dimension in QCD and its supersymmetric extensions_ , _JHEP_ 01 (2016) 140 [1510.07803].
  * (9) R. Brüser, C. Dlapa, J. M. Henn and K. Yan, _Full Angle Dependence of the Four-Loop Cusp Anomalous Dimension in QED_ , _Phys. Rev. Lett._ 126 (2021) 021601 [2007.04851].
  * (10) A. Mitov, G. F. Sterman and I. Sung, _The Massive Soft Anomalous Dimension Matrix at Two Loops_ , _Phys. Rev. D_ 79 (2009) 094015 [0903.3241].
  * (11) A. Ferroglia, M. Neubert, B. D. Pecjak and L. L. Yang, _Two-loop divergences of scattering amplitudes with massive partons_ , _Phys. Rev. Lett._ 103 (2009) 201601 [0907.4791].
  * (12) A. Mitov, G. F. Sterman and I. Sung, _Computation of the Soft Anomalous Dimension Matrix in Coordinate Space_ , _Phys. Rev. D_ 82 (2010) 034020 [1005.4646].
  * (13) Y.-T. Chien, M. D. Schwartz, D. Simmons-Duffin and I. W. Stewart, _Jet Physics from Static Charges in AdS_ , _Phys. Rev. D_ 85 (2012) 045010 [1109.6010].
  * (14) S. Catani and M. H. Seymour, _A General algorithm for calculating jet cross-sections in NLO QCD_ , _Nucl. Phys. B_ 485 (1997) 291 [hep-ph/9605323].
  * (15) A. Bassetto, M. Ciafaloni and G. Marchesini, _Jet Structure and Infrared Sensitive Quantities in Perturbative QCD_ , _Phys. Rept._ 100 (1983) 201.
  * (16) A. V. Kotikov, _Differential equations method: New technique for massive Feynman diagrams calculation_ , _Phys. Lett. B_ 254 (1991) 158.
  * (17) Z. Bern, L. J. Dixon and D. A. Kosower, _Dimensionally regulated pentagon integrals_ , _Nucl. Phys. B_ 412 (1994) 751 [hep-ph/9306240].
  * (18) T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys. B_ 580 (2000) 485 [hep-ph/9912329].
  * (19) J. M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601 [1304.1806].
  * (20) J. Henn and D. Simmons-Duffin, _Unpublished_ , 2011.
  * (21) C. W. Milloy, _Infrared divergences in scattering amplitudes from correlators of Wilson Lines_ , Ph.D. thesis, Edinburgh U., 2020.  10.7488/era/771.
  * (22) J. G. M. Gatheral, _Exponentiation of Eikonal Cross-sections in Nonabelian Gauge Theories_ , _Phys. Lett._ 133B (1983) 90.
  * (23) J. Frenkel and J. C. Taylor, _Nonabelian Eikonal Exponentiation_ , _Nucl. Phys._ B246 (1984) 231.
  * (24) G. F. Sterman, _Infrared divergences in perturbative QCD_ , _AIP Conf. Proc._ 74 (1981) 22.
  * (25) A. V. Manohar and M. B. Wise, _Heavy Quark Physics_ , Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology. Cambridge University Press, 2000, 10.1017/CBO9780511529351.
  * (26) G. Falcioni, E. Gardi, M. Harley, L. Magnea and C. D. White, _Multiple Gluon Exchange Webs_ , _JHEP_ 10 (2014) 010 [1407.3477].
  * (27) G. P. Korchemsky and A. V. Radyushkin, _Renormalization of the Wilson Loops Beyond the Leading Order_ , _Nucl. Phys._ B283 (1987) 342.
  * (28) A. J. N. Waelkens, _Calculation of webs in non-Abelian gauge theories using unitarity cuts_ , Ph.D. thesis, Edinburgh U., Edinburgh U., 2017.
  * (29) J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, _Integral reduction with Kira 2.0 and finite field methods_ , _Comput. Phys. Commun._ 266 (2021) 108024 [2008.06494].
  * (30) C. Meyer, _Algorithmic transformation of multi-loop master integrals to a canonical basis with CANONICA_ , _Comput. Phys. Commun._ 222 (2018) 295 [1705.06252].
  * (31) J. Henn, B. Mistlberger, V. A. Smirnov and P. Wasser, _Constructing d-log integrands and computing master integrals for three-loop four-particle scattering_ , _JHEP_ 04 (2020) 167 [2002.09492].
  * (32) C. Dlapa, J. Henn and K. Yan, _Deriving canonical differential equations for Feynman integrals from a single uniform weight integral_ , _JHEP_ 05 (2020) 025 [2002.02340].
  * (33) L. J. Dixon, E. Gardi and L. Magnea, _On soft singularities at three loops and beyond_ , _JHEP_ 02 (2010) 081 [0910.3653].
  * (34) R. N. Lee, _Libra: A package for transformation of differential systems for multiloop integrals_ , _Comput. Phys. Commun._ 267 (2021) 108058 [2012.00279].
  * (35) Z. L. Liu and N. Schalch, _Infrared Singularities of Multileg QCD Amplitudes with a Massive Parton at Three Loops_ , _Phys. Rev. Lett._ 129 (2022) 232001 [2207.02864].
  * (36) D. Correa, J. Henn, J. Maldacena and A. Sever, _The cusp anomalous dimension at three loops and beyond_ , _JHEP_ 05 (2012) 098 [1203.1019].