aainstitutetext: Bethe Center for Theoretical Physics, Universität Bonn,
D-53115, Germanybbinstitutetext: SLAC National Accelerator Laboratory,
Stanford University, Stanford, CA 94039, USA

# Lepton-pair production at hadron colliders at N3LO in QCD

Claude Duhr<EMAIL_ADDRESS>b Bernhard Mistlberger
<EMAIL_ADDRESS>
###### Abstract

We compute for the first time the complete corrections at N3LO in the strong
coupling constant to the inclusive neutral-current Drell-Yan process including
contributions from both photon and $Z$-boson exchange. Our main result is the
computation of the QCD corrections to the inclusive production cross section
of an axial-vector boson to third order in the strong coupling in a variant of
QCD with five massless quark flavours. Since the axial anomaly does not cancel
for an odd number of flavours, we also consistently include non-decoupling
effects in the top-quark mass through three loops. We perform a
phenomenological study of our results, and we present for the first time
predictions for the inclusive Drell-Yan process at the LHC at this order in
QCD perturbation theory.

###### Keywords:

N3LO, Drell-Yan production.

††preprint:  BONN-TH-2021-12, SLAC-PUB-17632

##  1 Introduction

The production of a pair of massless leptons at hadron colliders like the
Large Hadron Collider (LHC) at CERN – the so-called Neutral-Current Drell-Yan
(NCDY) process – is one of the most important and most studied hadron collider
processes. From an experimental perspective, its clean final-state signature
makes it an ideal candidate for luminosity measurements and detector
calibration. The DY process also plays a key role in the measurement of parton
distribution functions (PDFs) at the LHC and in many searches for physics
beyond the Standard Model (SM). From a theoretical perspective, the invariant-
mass distribution of the produced lepton pair is arguably the simplest hadron
collider observable, and it often serves as a template to understand the
structure of higher-order QCD corrections at hadron colliders more generally.
It is thus important to have a solid theoretical understanding for the NCDY
process, including higher-order corrections in both the strong and electroweak
coupling constants.

At leading order (LO) in the electroweak coupling constant, the massless
lepton pair is produced via the propagation of an intermediate, off-shell
photon or $Z$-boson Drell:1970wh . The next-to-leading order (NLO) QCD
corrections have been computed several decades ago Altarelli:1978id ;
Altarelli:1979ub . The next-to-next-to-leading order (NNLO) corrections were
computed in refs. Matsuura:1987wt ; Matsuura:1988nd ; Matsuura:1988sm ;
Hamberg:1990np ; Matsuura:1990ba ; vanNeerven:1991gh  in a version of the SM
where all quarks, including the top quark, are considered massless, and the
effects of finite quark masses were computed in refs. Dicus:1985wx ;
Gonsalves:1991qn ; Rijken:1995gi . Fully differential predictions for the
Drell-Yan process at NNLO are readily available, see for example refs.
Anastasiou:2003yy ; Cieri:2015rqa ; Gavin:2010az ; Boughezal:2016wmq ;
Grazzini:2017mhc . Electroweak corrections, including mixed QCD-electroweak
corrections are also available Delto:2019ewv ; Bonciani:2019nuy ;
Bonciani:2020tvf ; Bonciani:2021zzf ; Buccioni:2020cfi .

The uncertainty on the NNLO results due to missing higher orders in the QCD
perturbative expansion was estimated to be at the percent level by varying the
renormalisation and factorisation scales by a factor of two around a central
scale related to the invariant mass of the lepton pair. Very recently, also
next-to-next-to-next-to-leading order (N3LO) corrections in the strong
coupling have been computed, but only for the contributions from an
intermediate off-shell-photon Duhr:2020seh . Unlike for N3LO corrections to
inclusive Higgs production processes Anastasiou:2015ema ; Anastasiou:2016cez ;
Duhr:2019kwi ; Chen:2019lzz ; Duhr:2020kzd ; Dreyer:2018qbw ; Dreyer:2016oyx ,
it was found that for large ranges of invariant masses, the N3LO corrections
are sizeable and lower the value of the cross section by a few percent, which
is more than the missing higher order uncertainty estimated from varying the
perturbative scales at NNLO. This result was recently confirmed by the
independent computation from the double-differential computation for lepton-
pair production via an off-shell photon at N3LO in ref. Chen:2021vtu  and the
computation of the fiducial cross setion at N3LO in ref. Camarda:2021ict . A
similar behaviour was observed for the production of a lepton-neutrino pair at
N3LO Duhr:2020sdp . This shows that N3LO corrections are highly needed if we
want to reach a precision at the percent level for vector-boson production at
hadron colliders. In particular, a complete calculation of the NCDY process in
the SM at N3LO, including also the contributions from the exchange of the
$Z$-boson, is highly desired.

The computation of higher-order corrections including $Z$-bosons, however, is
much more challenging. Unlike the photon, the $Z$-boson has both vector and
axial couplings to fermions. Higher-order computations typically diverge, and
they are conventionally regularised using dimensional regularisation, where
the space-time dimension is analytically continued from $D=4$ to
$D=4-2\epsilon$ dimensions. Axial couplings involve a $\gamma^{5}$ matrix,
which is an intrinsically four-dimensional object, and so its treatment in
dimensional regularisation is ambiguous (see ref. Gnendiger:2017pys  and
references therein for a review). Moreover, it is well known that the axial
current is anomalous in QCD, and its divergence is described by the famous
Adler-Bell-Jackiw (ABJ) anomaly equation Adler:1969gk ; Bell:1969ts ;
Adler:1969er  (this anomaly is absent for the vector current, as a consequence
of Yang’s theorem Yang:1950rg ). In the complete SM, the axial anomaly
cancels. Calculations involving massive top quarks, however, are technically
extremely challenging, and therefore computations are usually done in an
effective theory where the top quark is infinitely heavy and is integrated
out. This procedure does no longer naively work in the presence of an axial
current, because the resulting effective theory is anomalous Collins:1978wz .

The goal of this paper is to present for the first time complete results for
the invariant mass distribution for the NCDY process in the SM at N3LO in the
strong coupling constant, including contributions from both an intermediate
photon and $Z$-boson. Our main contribution is the computation the inclusive
N3LO cross section for an axial-vector state in QCD. We treat the $\gamma^{5}$
matrix by working in the Larin scheme tHooft:1972tcz ; Breitenlohner:1977hr ;
Larin:1991tj ; Larin:1993tq ; Jegerlehner:2000dz  and we work in a theory with
five massless active quark flavours where the top quark has been integrated
out. Since this effective theory is anomalous, the top quark does not
completely decouple. We include non-decoupling effects via a Wilson
coefficient multiplying the axial current in the effective theory. As a
result, we obtain the complete N3LO cross section in the SM, up to terms that
are power-suppressed in the top quark mass. The finite top-mass effects are
known to be small at NNLO, and we expect the size of the missing power-
suppressed terms to be negligible at N3LO as well. As a result, we obtain for
the first time complete phenomenological predictions for the NCDY process to
third order in the strong coupling.

This paper is organised as follows: In section 2 we review the structure of
the QCD corrections to the NCDY process. In section 3 we present the main
ingredient of our computation, namely the cross section for axial-vector
production at N3LO in QCD, and we discuss our treatment of the $\gamma^{5}$
matrix in dimensional regularisation and the non-decoupling top-mass
contributions. In section 4 we discuss the phenomenological implications of
our results. In section 5 we draw our conclusions.

##  2 The Neutral-Current Drell-Yan process

The Neutral-Current Drell-Yan (NCDY) process describes the production of a
pair of (massless) leptons as the result of the collision of (anti-) protons:

$\text{P}(P_{1})\,+\,\text{P}(P_{2})\to l\,\bar{l}(Q)+X\,.$ (1)

Here, $P_{1}$ and $P_{2}$ are the momenta of the scattering protons, $Q$ is
the invariant mass of the lepton pair $l\,\bar{l}$, and $X$ collectively
denotes additional QCD radiation. The goal of this paper is to compute this
process to third order in the strong coupling constant and to discuss the
associated phenomenology. More precisely, we want to compute the cross section
for the process in eq. (1) differentially in the invariant mass $Q$. The
production cross section is mediated by a virtual photon or $Z$-boson which
subsequently decay to the final-state leptons. Our formalism is accurate up to
corrections that are suppressed in the electroweak coupling constant
$\alpha_{EW}$:

$Q^{2}\frac{\mathrm{d}\sigma_{P\,P\rightarrow
l\bar{l}+X}}{\mathrm{d}Q^{2}}=\sum\limits_{B,B^{\prime}\in\\{\gamma^{*},Z\\}}Q^{2}\frac{\mathrm{d}\sigma_{P\,P\to
B/B^{\prime}+X\to
l\bar{l}+X}}{\mathrm{d}Q^{2}}+\mathcal{O}\left(\alpha_{EW}^{3}\right).$ (2)

Figure 1:  Schematic depiction of the Drell-Yan process and its factorisation
into the production probability of a virtual gauge boson and the subsequent
decay to final-state leptons.

Figure 1 schematically shows the NCDY process for the production of a lepton
pair via an intermediate vector boson. The cross section can be written in a
factorised form,

$Q^{2}\frac{\mathrm{d}\sigma_{P\,P\to B/B^{\prime}+X\to
l\,\bar{l}+X}}{\mathrm{d}Q^{2}}=Q^{2}\frac{\mathrm{d}\sigma_{P\,P\rightarrow
B/B^{\prime}+X}}{\mathrm{d}Q^{2}}\times\Gamma_{B/B^{\prime}\to
l\,\bar{l}}(Q^{2})\times\text{BW}(Q^{2},m_{B},m_{B^{\prime}})\,.$ (3)

We now describe the factors in the right-hand side of eq. (3) in detail.

The propagation of the virtual gauge bosons is described by the Breit-Wigner
distribution,

$\text{BW}(Q^{2},m_{B},m_{B}^{\prime})=\frac{Q^{3}}{\pi}\mathtt{Re}\left[\frac{1}{(Q^{2}-m_{B}^{2})+im_{B}\Gamma_{B}}\frac{1}{(Q^{2}-m_{{B^{\prime}}}^{2})-im_{B^{\prime}}\Gamma_{B^{\prime}}}\right]\,.$
(4)

The decay of the gauge bosons into leptons is captured by the interference
contributions to the width:

$\displaystyle\Gamma_{B/B^{\prime}\to l\bar{l}}(Q,m_{l},m_{\bar{l}})=$ (5)
$\displaystyle\,\,\frac{1}{6Q}\int\frac{\mathrm{d}^{d}p_{l}}{(2\pi)^{d-1}}\delta_{+}(p_{l}^{2}-m_{l}^{2})\frac{\mathrm{d}^{d}p_{\bar{l}}}{(2\pi)^{d-1}}\delta_{+}(p_{\bar{l}}^{2}-m_{\bar{l}}^{2})(2\pi)^{d}\delta^{d}(q-p_{l}-p_{\bar{l}})\mathcal{M}_{B\to
l\bar{l}}\cdot\mathcal{M}_{B^{\prime}\to l\bar{l}}^{*}\,,$

We use the notation $\mathcal{M}_{1}\cdot\mathcal{M}_{2}^{*}$ to indicate that
the interferences are summed over the colours and spins of all external
particles. Equation (5) describes the interferences of the decay amplitudes of
gauge bosons into two final-state particles with masses $m_{l}$ and
$m_{\bar{l}}$ at Born level. It is chosen such that if the formula is
evaluated on-shell, $Q=m_{B}$ and $B=B^{\prime}$, then it corresponds to the
partial width $\Gamma_{B\to l\bar{l}}$ of the gauge boson $B$ in its rest-
frame. Since we are working to leading order in the electroweak coupling, we
only need to consider the tree-level contributions to the decay. For massless
leptons, we find:

$\Gamma_{B/B^{\prime}\to l\bar{l}}(Q^{2}):=\Gamma_{B/B^{\prime}\to
l\bar{l}}(Q^{2},0,0)=\frac{Qm_{W}^{2}}{3\pi
v^{2}}\left[Q_{V,l}^{B}Q_{V,l}^{B^{\prime}}+Q_{A,l}^{B}Q_{A,l}^{B^{\prime}}\right]\,.$
(6)

The charges in the above equation for a $Z$ boson or photon are defined as :

$\displaystyle Q_{V,f}^{\gamma^{*}}$ $\displaystyle=$ $\displaystyle
Q_{f}\sqrt{1-\frac{m_{W}^{2}}{m_{Z}^{2}}}\,,\hskip
28.45274ptQ_{A,f}^{\gamma^{*}}=0\,,$ (7) $\displaystyle Q_{V,f}^{Z}$
$\displaystyle=$
$\displaystyle\frac{m_{Z}}{2m_{W}}\left[T_{f}^{3}-2Q_{f}\left(1-\frac{m_{W}^{2}}{m_{Z}^{2}}\right)\right]\,,\hskip
28.45274ptQ_{A,f}^{Z}=\frac{m_{Z}}{2m_{W}}T_{f}^{3}\,.$

Here, $m_{W}$ and $m_{Z}$ are the masses of the $W$ and $Z$ bosons. In the SM
we have the following charge assignments for up-type and down-type quarks and
electrically-charged leptons:

$\begin{array}[]{c|c|c|c}\hline\cr\hline\cr&\hskip 14.22636ptu&\hskip
14.22636ptd&\hskip 14.22636ptl\\\ \hline\cr
Q_{f}&\frac{2}{3}&-\frac{1}{3}&-1\\\
T_{f}^{3}&\frac{1}{2}&-\frac{1}{2}&-\frac{1}{2}\\\
\hline\cr\hline\cr\end{array}$ (8)

Finally, the first factor in eq. (3) describes the production cross section
for the interference of two off-shell gauge bosons $B$ and $B^{\prime}$:

$\begin{split}Q^{2}\frac{\mathrm{d}\sigma_{P\,P\rightarrow
B/B^{\prime}+X}}{\mathrm{d}Q^{2}}=\tau\int_{0}^{1}\mathrm{d}x_{1}\,\mathrm{d}x_{2}\,\mathrm{d}z\,\delta(\tau-
x_{1}x_{2}z)\,\sum_{i,j}f_{i}(x_{1})\,f_{j}(x_{2})\,\eta_{\,ij\to
B/B^{\prime}+X}(z)\,.\end{split}$ (9)

In the previous equation we introduced the parton distribution functions
$f_{i}(x)$ that are convoluted with the partonic coefficient functions
$\eta_{\,ij\to B/B^{\prime}+X}$. We suppress the dependence of all quantities
on the factorisation scale to improve the readability. The variables $x_{1,2}$
represent the fraction of the momenta of the protons carried by the initial-
state partons. We introduce the variables

$S=(P_{1}+P_{2})^{2},\hskip 28.45274pt\tau=\frac{Q^{2}}{S}.$ (10)

The partonic coefficient functions are given by

$\eta_{\,ij\to
B/B^{\prime}+X}=\frac{\mathcal{N}_{ij}}{2Q^{2}}\sum_{X_{i}}\int\mathrm{d}\Phi_{B/B^{\prime}+X_{i}}\mathcal{M}_{ij\to
B+X_{i}}\cdot\mathcal{M}_{ij\to B^{\prime}+X_{i}}^{*}\,,$ (11)

where the factor $\mathcal{N}_{ij}$ represents a process-dependent averaging
over initial-state spins and colours, $\mathrm{d}\Phi_{B/B^{\prime}+X_{i}}$ is
the phase space measure for a particular final state, and $\mathcal{M}_{ij\to
B+X_{i}}$ is the scattering matrix element for the production of this final
state. In order to compute the partonic coefficient functions to a given order
in perturbation theory, we consider a perturbative expansion of the product of
matrix elements truncated at a fixed order in the coupling constant. The
partonic coefficient functions are the main ingredients to our computation, as
they incorporate the entirety of the higher-order QCD corrections. We will
study their structure in more detail in the remainder of this section.

The coupling of the $Z$-boson to fermions involves a vector and an axial-
vector part (see figure 2).

Figure 2:  Feynman rule coupling a Z boson to fermions.

Accordingly, we may split the partonic coefficient function into two terms
accounting for the coupling to the vector and axial-vector parts separately:

$\eta_{\,ij\to B/B^{\prime}+X}=\hat{\sigma}\left[\eta^{V}_{\,ij\to
B/B^{\prime}+X}+\eta^{A}_{\,ij\to B/B^{\prime}+X}\right],\hskip
28.45274pt\hat{\sigma}=\frac{4\pi m_{W}^{2}}{n_{c}Q^{2}v^{2}}.$ (12)

The inclusive cross section does not contain any non-zero partonic coefficient
functions that depend on both a vector and an axial-vector coupling. We are
interested in the perturbative coefficient functions computed up to a fixed
order in the expansion in the strong coupling constant:

$\eta^{Y}_{\,ij\to
B/B^{\prime}+X}=\sum_{i=0}^{\infty}a_{s}(\mu)^{i}\eta^{Y,\,(i)}_{\,ij\to
B/B^{\prime}+X},\hskip 28.45274ptY\in\\{V,A\\}\,,\hskip
28.45274pta_{s}(\mu)=\frac{\alpha_{S}(\mu)}{\pi}\,,$ (13)

where $\alpha_{S}(\mu)$ denotes the strong coupling constant in the
$\overline{\textrm{MS}}$ scheme. At leading order in perturbation theory only
the partonic coefficient functions with a quark – anti-quark pair in the
initial state are non-vanishing:

$\eta^{Y,\,(0)}_{q\bar{q}\to
B/B^{\prime}}=Q_{Y,q}^{B}Q_{Y,q}^{B^{\prime}}\delta\left(1-z\right),\hskip
28.45274ptY\in\\{V,A\\}.$ (14)

(a) $Q_{Y,i}^{B}Q_{Y,j}^{B^{\prime}}\delta_{ij}$

(b) $Q_{Y,i}^{B}Q_{Y,i}^{B^{\prime}}$

(c) $Q_{Y,i}^{B}Q_{Y,j}^{B^{\prime}}$

(d) $Q_{Y,j}^{B}\sum_{f}Q_{Y,f}^{B^{\prime}}$

(e) $\delta_{ij}Q_{Y,i}^{B}\sum_{f}Q_{Y,f}^{B^{\prime}}$

(f) $\sum_{f}Q_{Y,f}^{B}Q_{Y,f}^{B^{\prime}}$

Figure 3:  Representative diagrams for different charge structures appearing
in the NCDY cross section. The diagrams show interference diagrams where the
vertical dashed line represents the final state phase space cut. Solid, wavy
and curly lines represent propagating quarks, gauge bosons and gluons
respectively.

The partonic coefficient functions for the production of an electroweak gauge
boson can be split into 9 different contributions, corresponding to particular
combinations of the charges associated with different ways of coupling the
electroweak gauge bosons to the quarks. For the purpose of simplicity, we
write $\eta_{ij}=\eta^{Y}_{\,ij\to B/B^{\prime}+X}$. We then have the
following decomposition:

$\displaystyle\eta_{ij}$ $\displaystyle=$
$\displaystyle\eta_{ij}^{1}\left(\sum_{f}Q_{Y,f}^{B}\right)\left(\sum_{f}Q_{Y,f}^{B^{\prime}}\right)+\eta_{ij}^{2}\sum_{f}Q_{Y,f}^{B}Q_{Y,f}^{B^{\prime}}+\eta_{ij}^{3}Q_{Y,i}^{B}\sum_{f}Q_{Y,f}^{B^{\prime}}+\eta_{ij}^{4}Q_{Y,j}^{B}\sum_{f}Q_{Y,f}^{B^{\prime}}$
(15) $\displaystyle+$
$\displaystyle\eta_{ij}^{5}Q_{Y,i}^{B}Q_{Y,i}^{B^{\prime}}+\eta_{ij}^{6}Q_{Y,j}^{B}Q_{Y,j}^{B^{\prime}}+\eta_{ij}^{7}Q_{Y,i}^{B}Q_{Y,j}^{B^{\prime}}\delta_{ij}+\eta_{ij}^{8}Q_{Y,i}^{B}Q_{Y,j}^{B^{\prime}}+\eta_{ij}^{9}\delta_{ij}Q_{Y,i}^{B}\sum_{f}Q_{Y,f}^{B^{\prime}}\,.$

Figure 3 shows representative diagrams that illustrate how the different
charge structures arise. Figure 3(a) shows the tree-level diagram where the
initial-state $q\bar{q}$-pair annihilates and produces the gauge boson. The
flavors of the initial-state quarks are connected by the gauge boson. This is
made manifest by the Kronecker delta symbol $\delta_{ij}$ multiplying the
corresponding partonic coefficient function $\eta_{ij}^{7}$. In contrast,
there are classes of diagrams with a quark and anti-quark of the same flavor,
but the quark lines are not identical. In such diagrams the gauge boson may or
may not connect the different quark lines, as is illustrated in fig. 3(b) and
fig. 3(c). Their contributions are included in $\eta_{ij}^{5}$,
$\eta_{ij}^{6}$ and $\eta_{ij}^{8}$ respectively. Figures 3(d) and 3(e) show
contributions, denoted by by $\eta_{ij}^{3}$, $\eta_{ij}^{4}$ and
$\eta_{ij}^{9}$, where the gauge boson is connected only to one quark line of
an initial-state quark. Furthermore, it is possible that the gauge boson
couples to quark lines that are disconnected from the initial state, as shown
for example in fig. 3(f). These contributions are taken into account by
$\eta_{ij}^{1}$ and $\eta_{ij}^{2}$.

The partonic coefficient functions $\eta^{V}_{\,ij\to B/B^{\prime}+X}$ and
$\eta^{A}_{\,ij\to B/B^{\prime}+X}$ begin to differ starting from NNLO in QCD
perturbation theory. At this order, it is for the first time possible to
connect two different quark lines with the produced electroweak gauge boson.
As long as the gauge boson is twice connected to the same quark line and all
quarks are treated as massless, helicity conservation dictates that vector and
axial-vector partonic coefficient functions will be identical. In particular,
$\eta_{ij}^{2}$, $\eta_{ij}^{5}$, $\eta_{ij}^{6}$, $\eta_{ij}^{7}$ are
identical for vector and axial-vector production at all orders in perturbation
theory.

The partonic coefficient functions for the production of a vector boson
$\eta^{V}_{\,ij\to B/B^{\prime}+X}$ were computed at N3LO in QCD perturbation
theory in refs. Altarelli:1978id ; Altarelli:1979ub ; Matsuura:1987wt ;
Matsuura:1988nd ; Matsuura:1988sm ; Matsuura:1990ba ; Hamberg:1990np ;
vanNeerven:1991gh ; Harlander:2002wh ; Duhr:2020seh . The computation of the
partonic coefficient functions for an axial-vector boson $\eta^{A}_{\,ij\to
B/B^{\prime}+X}$ through N3LO in perturbative QCD is one of the main results
of this article. We present analytic formulæ for these partonic coefficient
functions in terms of ancillary files attached to the arXiv submission of this
article. In the next section we discuss in detail the computation and the
structure of the partonic coefficient functions for axial-vector production in
QCD.

##  3 Axial-vector production at N3LO in QCD

Let us consider N3LO corrections to the production cross section for a color-
singlet axial-vector state $Z_{A}$. In broad terms, our computational strategy
is similar to the computation of the N3LO corrections to the inclusive Higgs
boson production cross sections in gluon fusion Anastasiou:2015ema ;
Mistlberger:2018etf  and bottom-quark fusion Duhr:2019kwi , and the charged-
current and photon-only DY processes Duhr:2020sdp ; Duhr:2020seh . More
precisely, we generate Feynman diagrams using QGRAF Nogueira1993  and perform
the spinor and color algebra based on custom C++ algorithms using GiNaC
Bauer2000  and FORM Ruijl:2017dtg . We generate all required integrands for
interferences of matrix elements for this cross section up to N3LO. We then
reduce all real and virtual interference diagrams to a set of so-called master
integrals using IBP identities Tkachov1981 ; Chetyrkin1981 ; Laporta:2001dd
and the framework of reverse unitarity Anastasiou2002 ; Anastasiou2003 ;
Anastasiou:2002qz ; Anastasiou:2003yy ; Anastasiou2004a . The master integrals
were computed using the method of differential equations Kotikov:1990kg ;
Kotikov:1991hm ; Kotikov:1991pm ; Henn:2013pwa ; Gehrmann:1999as  in refs.
Anastasiou:2013mca ; Anastasiou:2013srw ; Anastasiou:2015yha ; Dulat:2014mda ;
Mistlberger:2018etf . We work in dimensional regularisation, and compute all
matrix elements in $D=4-2\epsilon$ dimensions. After inserting the master
integrals into our partonic coefficient functions, we renormalise UV
divergences in the $\overline{\text{MS}}$ scheme. Furthermore, we remove
collinear initial-state singularities via mass factorization counterterms
comprised of DGLAP Gribov:1972ri ; Altarelli:1977zs ; Dokshitzer:1977sg
splitting functions Moch:2004pa ; Vogt:2004mw ; Ablinger:2014nga ;
Ablinger:2017tan . Finally, we arrive at fully analytic formulæ for our finite
partonic coefficient functions for the production of an axial-vector boson
through N3LO in perturbative QCD. The functions are expressed in terms of the
class of iterated integrals introduced in ref. Mistlberger:2018etf . While the
general strategy is the same as for the N3LO cross section considered in ref.
Anastasiou:2015ema ; Mistlberger:2018etf ; Duhr:2019kwi ; Duhr:2020sdp , there
are several technical steps that distinguish axial-vector production from
those other processes, due to the ambiguities in how to treat $\gamma^{5}$ in
dimensional regularisation. In the remainder of this section we discuss our
treatment of $\gamma^{5}$ in detail.

Consider the operator describing the coupling of an axial current in QCD to
the axial-vector $Z_{A}$:

$\cal O_{A}=g_{A}\,Z_{A,\mu}\,\sum_{f=1}^{N_{f}}A_{f}\,J_{A,f}^{\mu}\,,\qquad
J_{A,f}^{\mu}=i\bar{q}_{f}\gamma^{\mu}\gamma^{5}q_{f}\,,$ (16)

where $N_{f}$ is the number of quark species. In the SM we have
$g_{A}=\frac{2m_{W}}{v\,\cos\theta_{W}}$, $A_{f}=\frac{1}{2}T^{3}_{f}$ and
$N_{f}=6$. Since the coupling of $Z_{A}$ to the axial current involves the
$\gamma^{5}$ matrix, which is only well-defined in four dimensions, care is
needed how this operator is analytically continued to arbitrary dimensions
when working in dimensional regularisation. Moreover, it is well-known that
the singlet axial current $J_{A,S}^{\mu}=\sum_{f=1}^{N_{f}}J_{A,f}^{\mu}$ is
anomalous in QCD and satisfies the Adler-Bell-Jackiw (ABJ) anomaly equation
Adler:1969gk ; Bell:1969ts ; Adler:1969er  (see also ref. Bos:1992nd ). In
four dimensions $\gamma^{\mu}$ and $\gamma^{5}$ anticommute, as one can easily
verify by writing $\gamma_{5}$ as in refs. tHooft:1972tcz ;
Breitenlohner:1977hr :

$\gamma^{5}=-\frac{i}{4!}\,\epsilon^{\mu\nu\rho\sigma}\,\gamma_{\mu}\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}\,.$
(17)

However, if the spacetime dimension is extended to $D$-dimension this no
longer holds true. Here we work in the Larin-scheme Larin:1991tj ;
Larin:1993tq ; Jegerlehner:2000dz  and define the axial-vector current
explicitly as

$J_{A,f}^{\mu}=\frac{1}{3!}\epsilon^{\mu\nu\rho\sigma}\bar{q}_{f}\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}q_{f}.$
(18)

For $D=4$, this definition is identical to the one of eq. (16) as on can
easily see by computing the anti-commutation relation

$\gamma^{\mu}\gamma_{5}=\frac{1}{2}\\{\gamma^{\mu},\gamma_{5}\\}=-\frac{i}{3!}\epsilon^{\mu\nu\rho\sigma}\gamma_{\nu}\gamma_{\rho}\gamma_{\sigma}.$
(19)

In our computation the $\gamma^{\mu}$ matrices are interpreted as
$D$-dimensional objects, and the Dirac algebra is performed in $D$-dimensions.
Any Feynman diagram in the Drell-Yan process mediated by an axial-vector
current will involve two insertions of the axial-vector current of eq. (18).
We use the identity

$\epsilon^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}\epsilon^{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}=\det\left(\begin{array}[]{cccc}g^{\mu_{1}\nu_{1}}&g^{\mu_{1}\nu_{2}}&g^{\mu_{1}\nu_{3}}&g^{\mu_{1}\nu_{4}}\\\
g^{\mu_{2}\nu_{1}}&g^{\mu_{2}\nu_{2}}&g^{\mu_{2}\nu_{3}}&g^{\mu_{2}\nu_{4}}\\\
g^{\mu_{3}\nu_{1}}&g^{\mu_{3}\nu_{2}}&g^{\mu_{3}\nu_{3}}&g^{\mu_{3}\nu_{4}}\\\
g^{\mu_{4}\nu_{1}}&g^{\mu_{4}\nu_{2}}&g^{\mu_{4}\nu_{3}}&g^{\mu_{4}\nu_{4}}\\\
\end{array}\right),$ (20)

which is valid in strictly $D=4$ dimensions, to contract Lorentz indices of
the Levi-Civita tensors, and we treat the metric tensors in the above equation
as $D$-dimensional. The above extension to $D$ space-time dimensions modifies
the computed cross sections at sub-leading order in the dimensional regulator.
In the presence of divergences, these modifications are propagated into finite
and singular terms of the bare partonic cross sections and their
renormalisation below.

It is well-known that the axial current $J_{A,f}^{\mu}$ is anomalous in QCD
and develops a UV-divergence starting from one-loop order. This UV-divergence
cancels in the complete SM, due to the relation

$\sum_{f=1}^{N_{f}=6}A_{f}=0\,.$ (21)

In variants of the SM with an odd number $N_{f}$ of fermion species, however,
eq. (21) is violated, and the UV-divergence does not cancel. The UV-divergence
can be removed by renormalising the axial current by the equation (valid
through at least three loops)

$\big{[}J_{A,f}^{\mu}\big{]}_{R}=Z_{ns}\,J_{A,f}^{\mu}+Z_{s}\,\sum_{f^{\prime}=1}^{N_{f}}J_{A,f^{\prime}}^{\mu}=Z_{ns}\,J_{A,f}^{\mu}+Z_{s}\,J_{A,S}^{\mu}\,,$
(22)

where $\big{[}J_{A,f}^{\mu}\big{]}_{R}$ is the renormalised axial current. We
see that the renormalisation mixes the axial currents for different flavors.
The non-singlet and singlet counterterms $Z_{ns}$ and $Z_{s}$ are not pure
$\overline{\textrm{MS}}$ counterterms, but they also include finite terms
whose purpose is to ensure that the renormalised singlet axial current
computed in the Larin-scheme satisfies the all-order ABJ anomaly equation
Adler:1969gk ; Bell:1969ts ; Adler:1969er :

$\partial_{\mu}\big{[}J_{A,S}^{\mu}\big{]}_{R}=\frac{\alpha_{S}(\mu^{2})}{8\pi}\,N_{f}\,\big{[}F\widetilde{F}\big{]}_{R}\,,$
(23)

with
$F\widetilde{F}=\epsilon_{\mu\nu\rho\sigma}\,\textrm{Tr}\big{(}F^{\mu\nu}F^{\rho\sigma}\big{)}$,
and the renormalised singlet axial current is

$\big{[}J_{A,S}^{\mu}\big{]}_{R}=Z_{S}\,J_{A,S}^{\mu}\,,\textrm{~{}~{}with~{}~{}}Z_{S}=Z_{ns}+N_{f}\,Z_{s}\,.$
(24)

Both the non-singlet and singlet counterterms are known to three loops in QCD
Larin:1993tq ; Ahmed:2015qpa ; Ahmed:2021spj ; Chen:2021rft :

$\displaystyle Z_{ns}$
$\displaystyle\,=1-a_{s}(\mu)\,C_{F}+a_{s}(\mu)^{2}\,C_{F}\,\left[\frac{11C_{A}-2N_{f}}{24\epsilon}+\frac{1}{144}\left(-107\
C_{A}+198C_{F}+2N_{f}\right)\right]$
$\displaystyle\,+a_{s}(\mu)^{3}\,C_{F}\,\Big{\\{}-\frac{\left(11C_{A}-2N_{f}\right)^{2}}{432\epsilon^{2}}+\frac{1}{2592\epsilon}\big{(}-416C_{A}N_{f}-2574C_{A}C_{F}+1789C_{A}^{2}$
$\displaystyle\,+360C_{F}N_{f}+4N_{f}^{2}\big{)}+\frac{1}{5184}\big{[}(4536\zeta_{3}-6441)\,C_{A}^{2}+(17502-12960\zeta_{3})\,C_{A}C_{F}$
$\displaystyle\,+(864\zeta_{3}+356)\,C_{A}N_{f}+(7776\zeta_{3}-9990)\,C_{F}^{2}-(864\zeta_{3}+186)\,C_{F}N_{f}+52N_{f}^{2}\big{]}\Big{\\}}$
$\displaystyle\,+\cal O(a_{s}(\mu)^{4})\,,$ (25) $\displaystyle Z_{s}$
$\displaystyle\,=a_{s}(\mu)^{2}\,C_{F}\left(\frac{3}{16\epsilon}+\frac{3}{32}\right)+a_{s}(\mu)^{3}\,C_{F}\Bigg{\\{}\frac{2N_{f}-11C_{A}}{96\epsilon^{2}}+\frac{109C_{A}-162\
C_{F}+2N_{f}}{576\epsilon}$
$\displaystyle\,+\frac{1}{3456}\big{[}(1404\zeta_{3}-326)\,C_{A}+(-1296\zeta_{3}+621)\,C_{F}+176N_{f}\big{]}\Bigg{\\}}+\cal
O(a_{s}(\mu)^{4})\,.$

We have computed the partonic coefficient functions in the Larin-scheme in a
variant of the SM with only $N_{f}=5$ massless active quark flavors, and the
top quark is considered infinitely-heavy and thus absent from the computation.
If only the strong coupling constant is UV-renormalised, then the partonic
coefficient functions still exhibit poles in the dimensional regulator
$\epsilon$ even after mass factorisation. However, all these poles cancel once
the axial current is renormalised as in eq. (22), and we obtain finite results
for all partonic coefficient functions. Besides the explicit analytic
cancellation of all ultraviolet and infrared poles, we performed several
additional checks to validate our results. In particular, we find agreement
with our previous computation of ref. Duhr:2020seh  of partonic coefficient
functions for vector boson production for the part of the partonic coefficient
functions that are identical between vector and axial-vector production.
Moreover, purely virtual corrections were computed in ref. Ahmed:2021spj ;
Gehrmann:2021ahy ; Gehrmann:2010ue ; Gehrmann:2010tu  and we find agreement.
Finally, we also checked that our results have the expected dependence on the
perturbative scale $\mu=\mu_{F}=\mu_{R}$. This check, however, involves a
subtle point, which we discuss in detail in the remainder of this section.

Our goal is to compute the coefficient functions $\eta_{ij}^{A}$, which
correspond to the axial-vector contribution from $Z$-boson exchange in the
limit where we drop all power-suppressed terms in the top-quark mass. If we
denote the partonic coefficient functions obtained from our procedure by
$\eta_{ij}^{A,5}$, we find that generically
$\eta_{ij}^{A,5}\neq\eta_{ij}^{A}$. To see this, consider the dependence of
the coefficient functions on the perturbative scale $\mu:=\mu_{F}=\mu_{R}$. In
the SM, this dependence is governed by the the well-known Dokshitzer-Gribov-
Lipatov-Altarelli-Parisi (DGLAP) equation Gribov:1972ri ; Altarelli:1977zs ;
Dokshitzer:1977sg :

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\log\mu^{2}}\eta^{A}_{ij}$
$\displaystyle\,=-P_{ik}\otimes\eta^{A}_{kj}-P_{jk}\otimes\eta^{A}_{ik}-\beta\,a_{s}\,\frac{\partial}{\partial
a_{s}}\eta^{A}_{ij}\,,$ (26)

where $\beta=-\sum_{n=0}^{\infty}\beta_{n}\,a_{s}^{n+1}$ denotes the QCD
$\beta$-function and $P_{kl}$ is the DGLAP splitting kernel, and we introduced
the convolution:

$(f\otimes
g)(z)=\int_{0}^{1}\mathrm{d}x_{1}\,\mathrm{d}x_{2}\,\delta(z-x_{1}x_{2})\,f(x_{1})\,g(x_{2})\,.$
(27)

The coefficient functions that we computed, however, satisfy the evolution
equation

$\begin{split}\frac{\mathrm{d}}{\mathrm{d}\log\mu^{2}}\eta^{A,5}_{ij}&\,=-P_{ik}\otimes\eta^{A,5}_{kj}-P_{jk}\otimes\eta^{A,5}_{ik}-\beta\,a_{s}\,\frac{\partial}{\partial
a_{s}}\eta^{A,5}_{ij}\\\ &\,-N_{A}\,\gamma_{J}\sum_{1\leq f_{1},f_{2}\leq
N_{f}}(A_{f_{1}}+A_{f_{2}})\,\eta^{A,5}_{ij,f_{1}f_{2}}\,,\end{split}$ (28)

where we defined implicitly

$\eta^{A,5}_{ij}=\sum_{1\leq f_{1},f_{2}\leq
N_{f}}A_{f_{1}}\,A_{f_{2}}\,\eta^{A,5}_{ij,f_{1}f_{2}}\,,$ (29)

and we find it useful to introduce the quantity
$N_{A}=\sum_{f=1}^{N_{f}}A_{f}$. Note that $N_{A}=0$ in the complete SM with
$N_{f}=6$ quark species. $\gamma_{J}$ is the anomalous dimension of the
(renormalised) axial current Larin:1993tq ; Ahmed:2021spj :

$(-\epsilon+\beta)\,a_{s}(\mu)\,{\bf
Z}_{5}^{-1}\frac{\mathrm{d}}{\mathrm{d}\log\mu^{2}}{\bf
Z}_{5}=\gamma_{J}\,{\bf E}_{N_{f}}+\cal O(\epsilon)\,.$ (30)

Here we defined the matrix

${\bf Z}_{5}=Z_{ns}\,{\bf I}_{N_{f}}+Z_{s}\,{\bf E}_{N_{f}}\,,$ (31)

where ${\bf I}_{N_{f}}$ is the $N_{f}\times N_{f}$ unit matrix and ${\bf
E}_{N_{f}}$ is the $N_{f}\times N_{f}$ matrix whose entries are all 1. We have

$\begin{split}\gamma_{J}&\,=-a_{s}(\mu)^{2}\,N_{A}\,\frac{3C_{F}}{8}-a_{s}(\mu)^{3}\,N_{A}\,\frac{C_{F}}{8}\,\Bigg{(}\frac{71}{12}\,C_{A}-\frac{9}{4}\,C_{F}-\frac{N_{f}}{6}\Bigg{)}+\cal
O(a_{s}(\mu)^{4})\,.\end{split}$ (32)

Clearly, the term proportional to $\gamma_{J}$ in eq. (28) is absent in the
complete SM, and so $\eta^{A,5}_{ij}$ cannot be identified with the
large-$m_{t}$ limit of the SM. Note that in the complete SM we have $N_{A}=0$,
and so $\gamma_{J}=0$.

The interpretation and resolution of this puzzle is as follows: If we start
from the complete SM, with $N_{f}=6$ and $m_{t}<\infty$, the UV-divergences
from the axial anomaly cancel. Naively, one would expect that by making the
top quark infinitely heavy, it decouples from the theory and we land on an
effective field theory with $N_{f}=5$ massless quarks, defined by simply
ignoring operators involving the top quark. This naive expectation, however,
is wrong: Since the UV divergences from the axial current do not cancel for
$N_{f}=5$, our EFT is anomalous and the top quark does not naively decouple
(see, e.g., ref. Collins:1978wz ). We include non-decoupling effects in the
form of a Wilson coefficient $C_{A,f}(\mu)=A_{f}+\cal O(a_{s}^{2})$:

$A_{f}\,J_{A,f}^{\mu}=iA_{f}\,\bar{q}_{f}\gamma^{\mu}\gamma^{5}q_{f}\to
J_{A,f,\textrm{eff}}^{\mu}=iC_{A,f}(\mu)\,\bar{q}_{f}\gamma^{\mu}\gamma^{5}q_{f}\,.$
(33)

The Wilson coefficient satisfies a renormalisation group equation which
precisely cancels the contribution from the anomalous dimension $\gamma_{J}$
in eq. (28):

$\frac{\mathrm{d}}{\mathrm{d}\log\mu^{2}}C_{A,f}(\mu)=-\gamma_{J}\,\sum_{f^{\prime}=1}^{N_{f}}C_{A,f^{\prime}}(\mu)\,.$
(34)

We see that the partonic coefficient functions $\eta^{A}_{ij}$ in the
large-$m_{t}$ limit of the SM are obtained from $\eta^{A,5}_{ij}$ by the
replacement $A_{f}\to C_{A,f}(\mu)$ (cf. eq. (29)):

$\eta_{ij}^{A}=\sum_{1\leq f_{1},f_{2}\leq
N_{f}}C_{A,f_{1}}(\mu)\,C_{A,f_{2}}(\mu)\,\eta^{A,5}_{ij,f_{1}f_{2}}\,.$ (35)

The Wilson coefficient admits the perturbative expansion:

$C_{A,f}(\mu)=A_{f}+N_{A}\,\sum_{n=2}^{\infty}\sum_{k=0}^{n-1}a_{s}(\mu)^{n}c_{A}^{(n,k)}\,\log^{k}\frac{m_{t}^{2}}{\mu^{2}}\,.$
(36)

It is easy to check that by combining the evolution equations (28) and (34),
the function $\eta_{ij}^{A}$ defined in eq. (35) satisfies the SM evolution
equation (26).

The Wilson coefficient depends explicitly on the (on-shell) top-quark mass,
making the non-decoupling of the top quark manifest. It was recently computed
through three loops in ref. Ju:2021lah ; Chen:2021rft . We have independently
obtained all the coefficients except for the non-logarithmic three-loop
coefficient $c^{(3,0)}_{A}$, and we find full agreement. Indeed, the
coefficients of the logarithms, i.e., the $c_{A}^{(n,k)}$ with $k\neq 0$, are
determined from the evolution equation (34). The non-logarithmic terms
$c_{A}^{(n,0)}$ are not constrained by the evolution equation and need to be
computed explicitly. We have recalculated the two-loop coefficient
$c^{(2,0)}_{A}$ by taking the large-$m_{t}$ limit on the known NNLO cross
section in the full SM, including all finite quark-mass effects Dicus:1985wx ;
Gonsalves:1991qn ; Rijken:1995gi . The two-loop coefficients are:

$\begin{split}c_{A}^{(2,1)}&\,=-\frac{3}{8}\,C_{F}\qquad\textrm{~{}~{}and~{}~{}}\qquad
c_{A}^{(2,0)}=-\frac{3}{16}\,C_{F}\,\,.\end{split}$ (37)

At three-loops, we have Ju:2021lah ; Chen:2021rft :

$\begin{split}c_{A}^{(3,2)}&\,=\frac{3}{8}\,C_{F}\,\beta_{0}=\frac{C_{F}}{32}(11\,C_{A}-2\,N_{f})\,,\\\
c_{A}^{(3,1)}&\,=-\frac{C_{F}}{96}\,(38C_{A}-27C_{F}+4N_{f})\,,\\\
c_{A}^{(3,0)}&\,=\frac{C_{F}}{1152}\big{[}C_{A}\,(1649-1512\,\zeta_{3})+27C_{F}\,(48\zeta_{3}+17)-374N_{f}-328\big{]}\,.\end{split}$
(38)

##  4 Phenomenological Results

In this section we present for the first time phenomenological results for the
invariant-mass distribution of a massless lepton pair at N3LO in QCD,

$\Sigma^{\textrm{N${}^{k}$LO}}(Q^{2})=Q^{2}\,\frac{\mathrm{d}\sigma^{\textrm{N${}^{k}$LO}}}{\mathrm{d}Q^{2}}\,.$
(39)

We work in the five-flavour scheme with $N_{f}=5$ active massless quark
flavors. The top-quark is considered infinitely heavy, though the non-
decoupling effects are included through the Wilson coefficient $C_{A,f}$. We
mention that finite quark-mass effects have been calculated through NNLO and
are known to be very small Dicus:1985wx ; Gonsalves:1991qn ; Rijken:1995gi .
We therefore expect our computation to give a very good estimate of the QCD
corrections in the full SM where the full finite top-mass effects are
retained. Our choice for the numerical values of the SM input parameters is
summarized in table 1. The strong coupling constant is evolved from
$\alpha_{S}(m_{Z})$ to the renormalisation scale $\mu_{R}$ using the four-loop
QCD beta function VanRitbergen:1997va ; Czakon:2004bu ; Baikov:2016tgj ;
Herzog:2017ohr  in the $\overline{\textrm{MS}}$-scheme using $N_{f}=5$ active
massless quark flavors. Unless stated otherwise, all results are obtained for
a proton-proton collider with $\sqrt{S}=13$ TeV using the zeroth member of the
combined `PDF4LHC15_nnlo_mc` set Butterworth:2015oua , and bands correspond to
varying the perturbative scales $\mu_{F}$ and $\mu_{R}$ by a factor of two
around the central scale $\mu_{F}=\mu_{R}=Q$ while respecting the constraint

$\frac{1}{2}\leq\frac{\mu_{R}}{\mu_{F}}\leq 2.$ (40)

Commonly, this is referred to as 7-point variation.

$\alpha_{EW}^{-1}$ | $\alpha_{S}(m_{Z})$ |  $v$ [GeV] |  $m_{Z}$ [GeV] |  $\Gamma_{Z}$ [GeV] |  $m_{W}$ [GeV] |  $m_{t}$ [GeV]
---|---|---|---|---|---|---
132.186 | 0.118 | 246.221 | 91.1876 | 2.4952 | 80.379 | 172.9
Table 1:  Numerical values of the SM input parameters used for our
phenomenological predictions.

###  4.1 Scale dependence and perturbative convergence

In table 2 we show values for a selection of representative invariant masses
$Q^{2}$ for the choice of central scale $\mu_{F}=\mu_{R}=Q$, together with the
corresponding QCD K-factors:

$\textrm{K}^{\text{N${}^{3}$LO}}(Q^{2})=\frac{\Sigma^{\textrm{N${}^{3}$LO}}(Q^{2})}{\Sigma^{\textrm{NNLO}}(Q^{2})}\,.$
(41)

We include an estimate of the residual perturbative uncertainty based on seven
point variation of the factorization and renormalisation scale, as well as the
uncertainties from PDFs and the value of the strong coupling constant (see
below).

$Q$ [GeV] |  $\Sigma^{\text{N3LO}}$ [pb] | $\textrm{K}^{\text{N${}^{3}$LO}}$ | $\delta(\text{PDF-$\alpha_{S}$})$ | $\delta(\text{PDF-TH})$
---|---|---|---|---
30 | $531.7^{+1.53\%}_{-2.54\%}$ | 0.952 | ${}^{+3.7\%}_{-3.8\%}$ | 2.8%
60 | $112.636^{+0.97\%}_{-1.29\%}$ | 0.97 | ${}^{+2.8\%}_{-2.5\%}$ | 2.5%
91.1876 | $21756.4^{+0.7\%}_{-0.86\%}$ | 0.977 | ${}^{+2.2\%}_{-2.1\%}$ | 2.5%
100 | $458.473^{+0.66\%}_{-0.79\%}$ | 0.979 | ${}^{+2.0\%}_{-1.8\%}$ | 2.5%
300 | $1.24661^{+0.26\%}_{-0.29\%}$ | 0.992 | ${}^{+1.9\%}_{-1.6\%}$ | 1.7%
Table 2: Cross section values for a selection of representative invariant mass
$Q$ for the choice of central scale $\mu_{F}=\mu_{R}=Q$, together with the
corresponding QCD K-factors and uncertainties (see main text for a detailed
description).

In figure 4 we show the inclusive cross section for the production of a
massless lepton pair as a function of $Q^{2}$.

Figure 4: The invariant-mass distribution $\Sigma(Q^{2})$ at the LHC with
$\sqrt{S}=13$ TeV at different orders in perturbation theory.

In figures 5 and 6 we show the NCDY cross section normalized to its value
computed at N3LO as a function of the invariant mass $Q$. We observe
qualitatively the same features as for the photon-only Duhr:2020seh ;
Chen:2021vtu  and the charged-current DY processes Duhr:2020sdp .
Specifically, we observe that in the range of invariant masses between $\sim
40$ GeV and $\sim 400$ GeV, the scale variation bands from NNLO and N3LO do
not overlap, indicating that conventional scale variation at NNLO
underestimated the true size of the N3LO corrections. We note, however, that
the size of the bands at NNLO was particularly small for the NCDY process,
often at the sub-percent level depending on the invariant masses considered.

In figure 7 we show the dependence of the cross section for $Q=100$ GeV on one
of the two perturbative scales with the other held fixed at some value in the
interval $[Q/2,2Q]$. We observe a very good reduction of the scale dependence
as we increase the perturbative order, with only a very mild scale dependence
at N3LO. Just like for the photon-only and $W$ cases, the bands from NNLO and
N3LO do not overlap. 111The leading order cross section does not depend on the
strong coupling constant and consequently does also not change with variation
of the renormalisation scale. As a result the right panel of fig. 7 does not
show any band for the leading order cross section.

Figure 5:  The K-factors
$\Sigma^{\textrm{N${}^{k}$LO}}/\Sigma^{\textrm{N${}^{3}$LO}}$ as a function of
invariant masses 10 GeV$\leq Q\leq$150 GeV for $k\leq 3$. The bands are
obtained by varying the perturbative scales by a factor of two around the
central $\mu_{\textrm{cent.}}=Q$.  Figure 6:  The K-factors
$\Sigma^{\textrm{N${}^{k}$LO}}/\Sigma^{\textrm{N${}^{3}$LO}}$ as a function of
invariant masses $Q\leq$1.800 GeV for $k\leq 3$. The bands are obtained by
varying the perturbative scales by a factor of two around the central
$\mu_{\textrm{cent.}}=Q$.

Figure 7:  Dependence of the invariant-mass distribution through N3LO on one
of the two perturbative scales with the other held fixed. The bands are
obtained by varying the other scale by a factor of two around the central
scale $Q=100$ GeV.

###  4.2 Ratios of K-factors

In the previous section we have seen that the NCDY process shows relatively
good perturbative stability as we increase the order. This behavior is very
reminiscent of the photon-only and charged-current DY processes considered in
refs. Duhr:2020seh ; Duhr:2020sdp . To investigate if and to what degree QCD
corrections differ among different DY-type processes, we study ratios of
K-factors between the NCDY process and the photon-only and charged-current DY
processes. In particular, in this section we refer to the ratio of the NnLO
cross section to its LO counter part as the NnLO K-factor, K(n). In contrast,
ratios of cross sections rather than K-factors involving the $W$ and
$\gamma^{*}$ cross sections were considered in ref. Duhr:2020sdp .

Before we discuss these ratios, we need to make a comment. Whenever one
studies ratios of cross sections, there is an ambiguity in how to choose the
perturbative scales. For example, if one believes that QCD corrections should
be similar between $W$, $\gamma^{*}$ and $Z$ production (as motivated for
example by the universality of certain limits, like the threshold limit), it
is natural to vary the scales in the numerator and the denominator in a
correlated way. Alternatively, one may choose the scales in an uncorrelated
way (e.g., because different partonic channels are weighted differently by the
PDFs for these processes, breaking the universality of the QCD corrections),
typically leading to larger scale variation bands. In ref. Duhr:2020sdp  it
was shown that for the $W$ and $\gamma^{*}$ cross sections the correlated
prescription leads to a vanishingly small scale dependence at the sub-permille
level at N3LO, while the uncorrelated prescription produces excessively large
scale variation bands at N3LO (much larger than the absolute shift from NNLO
to N3LO). As a consequence, neither the correlated nor the uncorrelated
prescriptions are expected to give reliable estimates for missing higher-order
terms for these ratios at N3LO. Therefore, in ref. Duhr:2020sdp  a new
prescription was considered, which uses the relative size of the last
considered order compared to the previous one as an estimator of the
perturbative uncertainty:

$\delta(\text{pert.})=\pm\left|1-\frac{\textrm{K}^{(n)}_{X}/\textrm{K}^{(n)}_{Y}}{\textrm{K}^{(n-1)}_{X}/\textrm{K}^{(n-1)}_{Y}}\right|\times
100\%\,.$ (42)

In the following we use this prescription to obtain uncertainty bands for the
ratios.

In figures 8 – 11 we show the ratios of the K-factors as a function of the
invariant mass $Q$. In all cases we observe a remarkable similarity of the
K-factors, which agree among themselves within at worst 2% for all ratios
considered. In particular, we see from figure 11 that the K-factors agree
within $1\%$ for the NCDY process computed with or without including the
contributions from the $Z$ boson. At the same time, we observe that there is a
dependence of the shape of the QCD corrections on the invariant mass, and
different invariant-mass regions may receive slightly different QCD
corrections. This shows that, if we want to reach a level of precision of 1%,
care is needed when using K-factors obtained from one process in one region of
phase space to reweight other processes or other regions of phase space.
However, our results also demonstrate that shape differences mainly introduced
at lower orders in perturbation theory and higher order corrections are
remarkably similar.

Figure 8:  Ratio of the K-factors for $W^{+}$ and $W^{-}$ production as a
function of $Q^{2}$ at different orders in perturbation theory.  Figure 9:
Ratio of the K-factors for NCDY and $W^{-}$ production as a function of
$Q^{2}$ at different orders in perturbation theory.  Figure 10:  Ratio of the
K-factors for NCDY and $W^{+}$ production as a function of $Q^{2}$ at
different orders in perturbation theory.  Figure 11:  Ratio of the K-factors
for NCDY and $\gamma^{*}$ production as a function of $Q^{2}$ at different
orders in perturbation theory.

###  4.3 Uncertainties related to PDFs

In order to assess the dependence of our predictions on the methodology of how
the PDFs are extracted, we follow the prescription of ref. Butterworth:2015oua
for the computation of PDF uncertainties $\delta(\text{PDF})$ using the Monte
Carlo method. The PDF set $\mathtt{PDF4LHC15\\_nnlo\\_mc}$ uses
$\alpha_{S}=0.118$ as a central value and two additional PDF sets are
available that allow for the correlated variation of the strong coupling
constant in the partonic cross section and the PDF sets to
$\alpha_{S}^{\text{up}}=0.1195$ and $\alpha_{S}^{\text{down}}=0.1165$. These
sets allow us to deduce an uncertainty $\delta(\alpha_{S})$ on our cross
section following the prescription of ref. Butterworth:2015oua . We combine
the PDF and strong coupling constant uncertainties in quadrature to give

$\delta(\text{PDF}+\alpha_{S})=\sqrt{\delta(\text{PDF})^{2}+\delta(\alpha_{S})^{2}}\,.$
(43)

Currently there is no available PDF set extracted from data with N3LO
accuracy, and so we are bound to use NNLO PDFs in our predictions. We estimate
the potential impact of this mismatch on our results using the prescription
introduced in ref. Anastasiou:2016cez . The PDF theory (PDF-TH) uncertainty is
then obtained by studying the variation of the NNLO cross section as NNLO- or
NLO-PDFs are used:

$\delta(\text{PDF-TH})=\frac{1}{2}\left|\frac{\Sigma^{\text{NNLO, NNLO-
PDFs}}(Q^{2})-\Sigma^{\text{NNLO, NLO-PDFs}}(Q^{2})}{\Sigma^{\text{NNLO, NNLO-
PDFs}}(Q^{2})}\right|.$ (44)

Here, the factor $\frac{1}{2}$ is introduced as it is expected that this
effect becomes smaller at N3LO compared to NNLO.

Figure 12:  Relative uncertainty of the NCDY process at N3LO due to incomplete
knowledge of parton distribution functions and the strong coupling constant as
a function of $Q$. $\delta(\textrm{PDF})$, $\delta(\textrm{PDF}+\alpha_{s})$
and the sum of $\delta(\textrm{PDF}+\alpha_{s})$ and $\delta(\textrm{PDF-TH})$
are shown in red, brown and green respectively.

In figure 12 we show the combined uncertainty from PDFs, the value of the
strong coupling constant $\alpha_{S}$ and the missing N3LO PDFs. The size of
these uncertainties is comparable to the uncertainties obtained in refs.
Duhr:2020seh ; Duhr:2020sdp  for the photon-only and charged-current DY
processes.

Figure 13:  Dependence of the NCDY process at N3LO on the choice of the PDF
set relative to the combined PDF4LHC15_nnlo_mc set. The red band corresponds
to the $\delta(\textrm{PDF})$ uncertainty.  Figure 14:  Dependence of the NCDY
process at N3LO on the choice of the PDF set relative to the combined
PDF4LHC15_nnlo_mc set. The red band corresponds to the $\delta(\textrm{PDF})$
uncertainty.

Figures 13 and 14 show the impact of evaluating the NCDY cross section with
different PDF sets. PDF4LHC15 is a combination of the CT14 Dulat:2015mca ,
MMHT14 Harland-Lang:2014zoa  and NNPDF3.0 Ball:2014uwa  PDF sets and we show
in fig. 13 predictions based on these individual sets relative to the
prediction based on the PDF4LHC15_nnlo_mc set. The red band in fig. 13
reflects the $\delta(\text{PDF})$ uncertainty of PDF4LHC15 and we observe that
the predictions based on the individual PDF sets are contained within this
band and that their spread is comparable in size to this band. Since the
publication of the PDF4LHC15 combination a plethora of developments and the
inclusion of LHC data into global PDF fits has led to updated PDF. In fig. 14
we study in particular the PDF sets ABMP16 Alekhin:2016uxn , CT18 Hou:2019efy
, MSHT20 Bailey:2020ooq , NNPDF3.1 Ball:2017nwa  and NNPDF4.0 Ball:2021leu .
We observe, that the spread among the newer sets is less than the the spread
of their predecessors. In particular, in the range where $Q$ is comparable to
the $Z$ and $W$-boson masses the different PDF sets seem to agree nicely. We
observe furthermore that NNPDF4.0 leads to a significant enhancement of the
NCDY cross section at large values of $Q$. It would be interesting to study
the impact of new PDF sets on high precision processes in further detail in
the future and we are looking forward to an updated version of PDF4LHC
Cridge:2021qjj .

###  4.4 Contributions from interferences

We conclude this phenomenological study by investigating some properties of
the NCDY cross section at N3LO and how it receives contributions from the
photon, the $Z$-boson and their interference.

Figure 15:  Decomposition of the invariant-mass distribution into the absolute
value of contributions from photon and $Z$ exchange and their interference as
a function of $Q$ at different order in perturbations theory.  Figure 16:
Relative contributions to the invariant-mass distribution from photon and $Z$
exchange and their interference as a function of $Q$ at different order in
perturbations theory.

In figures 15 and 16 we show the relative size of the three contributions. We
see that, as expected, the photon-only contribution dominates below the
$Z$-threshold. Above threshold all contributions are of similar size. Close to
the threshold the cross section is dominated entirely by the $Z$ boson, and
the interference changes sign when the threshold is crossed, as expected. The
observed pattern is hardly modified by perturbative corrections.

Computations with $\gamma^{5}$ matrices are technically more involved than for
purely vectorial couplings. For this reason, the axial-vector contribution is
frequently approximated by simply using the same partonic coefficient
functions as for the vector contribution and simply changing the coupling of
the quarks in an appropriate fashion. While this is correct for contributions
from Dirac traces with an even number of $\gamma^{5}$ matrices, a mismatch is
introduced starting from NNLO for contributions involving traces with an odd
number of $\gamma^{5}$ matrices. Moreover, effects from the axial-anomaly and
the resulting non-decoupling of the top-quark are neglected in this way.

Figure 17:  Relative difference of approximating the axial-vector production
cross section by a pure vector contribution as a function of $Q$ at different
order in perturbations theory.

In figure 17 we study the impact of replacing the axial-vector coefficient
functions $\eta^{A,i}_{\,ij\to A+X}$ in eq. (15) by their vector counterparts
$\eta^{V,i}_{\,ij\to Z+X}$ (while keeping the correct charges associated with
the axial-vector). The figure displays the ratio of the contribution to the DY
cross section due to the axial-vector current based on these ‘fake’ partonic
coefficient functions $\Sigma_{\text{Fake}}$ relativ to their correct
contribution $\Sigma_{\text{A}}$. We observe that at NNLO the mismatch
introduced is negligible, and it is well below the permille level for the
whole range of invariant masses considered. At N3LO instead, the mismatch
increases for small invariant masses, reaching the percent-level at invariant
masses below 50 GeV. However, the numerical impact grows with decreasing $Q$
as the the relative contribution of the $Z$-boson diminishes and the photon
exchange contribution dominates the DY production cross section, as can be
seen from fig. 16. We would like to point out that sufficiently differential
DY-type cross sections are afflicted by differences of the vector and axial-
vector contributions (see for example ref. Dixon:1997th ) that integrate to
zero in the inclusive cross section.

Figure 18:  Relative contribution from the non-decoupling top-mass effects as
a function of $Q$ at different order in perturbations theory.

In figure 18 we investigate the relative size of the non-decoupling of the
top-quark, i.e., the effect of including or not the Wilson coefficient
$C_{A,f}(\mu)$. We see that for the whole range of invariant masses
considered, the effect of including the Wilson coefficient or not is well
below the permille level.

##  5 Conclusions

In this paper we have computed for the first time the complete inclusive N3LO
corrections in the strong coupling constant for the production of a massless
lepton pair. This result has been made possible by combining our computation
for the N3LO corrections to vector production from ref. Duhr:2020seh  with the
computation of the axial-vector production cross section. A main ingredient in
our computation is the use of the Larin scheme to treat the $\gamma^{5}$
matrix in dimensional regularisation. Our computation has passed several non-
trivial checks: Apart from the explicit cancellation of all ultraviolet and
infrared poles, we have checked that our implementation of the Larin scheme
reproduces the contributions to the axial-vector cross section where both
$\gamma^{5}$ matrices appear in the same Dirac trace, and that our final
result satisfies the expected DGLAP renormalisation group equation once the
non-decoupling top-mass effects are included. We include all partonic
coefficient functions for the Drell-Yan production cross section as described
by eq. (15) as ancillary files together with the arXiv submission of this
article.

We have studied the phenomenological impact of our computation. We find that
the complete N3LO corrections to the NCDY process, including the axial-vector
contributions, has the same features as for the photon-only and $W$ cases of
refs. Duhr:2020seh ; Duhr:2020sdp . We find that the dependence of the
complete NCDY process on the perturbative scales is similar to the production
of a $\gamma^{*}$ or $W$. In particular, the bands obtained by varying the
perturbative scales by a factor of two do not overlap when going from NNLO to
N3LO, showing once more the importance of considering N3LO corrections for
precision LHC observables. We have also considered ratios of K-factors between
$\gamma^{*}$, $W$ and the complete NCDY process, and we find that, while the
K-factors are very similar, the ratio depends on the invariant mass of the
lepton pair, reaching a few percent depending on the invariant masses
considered. This shows that care is needed when taking K-factors from one
process and to apply it to another process, especially when aiming for
precision results. Finally, we have also studied the impact of the choice of
the PDF on our results. In particular, we have analysed how the missing N3LO
impact our predictions, which is one of the main bottlenecks to further
improve theoretical predictions for LHC observables. Since the NCDY process is
one of the main measurements used to constrain PDFs from LHC data, we expect
that our computation will play an important role in the future to determine
precisely the proton structure.

## Acknowledgements

We are grateful to W.-L. Ju, M. Schönherr, L. Chen, M. Czakon and M.
Niggetiedt for correspondence regarding the three-loop Wilson coefficient in
refs. Ju:2021lah ; Chen:2021rft . We are grateful to Lance Dixon for useful
discussions. BM was supported by the Department of Energy, Contract DE-
AC02-76SF00515.

## References

  * (1) S. D. Drell and T.-M. Yan, _Massive Lepton Pair Production in Hadron-Hadron Collisions at High-Energies_ , _Phys. Rev. Lett._ 25 (1970) 316.
  * (2) G. Altarelli, R. K. Ellis and G. Martinelli, _Leptoproduction and Drell-Yan Processes Beyond the Leading Approximation in Chromodynamics_ , _Nucl. Phys._ B143 (1978) 521.
  * (3) G. Altarelli, R. K. Ellis and G. Martinelli, _Large Perturbative Corrections to the Drell-Yan Process in QCD_ , _Nucl. Phys._ B157 (1979) 461.
  * (4) T. Matsuura and W. L. van Neerven, _Second Order Logarithmic Corrections to the Drell-Yan Cross-section_ , _Z. Phys._ C38 (1988) 623.
  * (5) T. Matsuura, S. C. van der Marck and W. L. van Neerven, _The Order $\alpha_{s}^{2}$ Contribution to the $K$ Factor of the Drell-Yan Process_, _Phys. Lett._ B211 (1988) 171.
  * (6) T. Matsuura, S. C. van der Marck and W. L. van Neerven, _The Calculation of the Second Order Soft and Virtual Contributions to the Drell-Yan Cross-Section_ , _Nucl. Phys._ B319 (1989) 570.
  * (7) R. Hamberg, W. L. van Neerven and T. Matsuura, _A complete calculation of the order $\alpha_{s}^{2}$ correction to the Drell-Yan $K$ factor_, _Nucl. Phys._ B359 (1991) 343.
  * (8) T. Matsuura, R. Hamberg and W. L. van Neerven, _The Contribution of the Gluon-gluon Subprocess to the Drell-Yan $K$ Factor_, _Nucl. Phys._ B345 (1990) 331.
  * (9) W. L. van Neerven and E. B. Zijlstra, _The $O(\alpha_{s}^{2})$ corrected Drell-Yan $K$ factor in the DIS and MS scheme_, _Nucl. Phys._ B382 (1992) 11.
  * (10) D. A. Dicus and S. S. D. Willenbrock, _Radiative Corrections to the Ratio of $Z$ and $W$ Boson Production_, _Phys. Rev._ D34 (1986) 148.
  * (11) R. J. Gonsalves, C.-M. Hung and J. Pawlowski, _Heavy quark triangle diagram contributions to Z boson production in hadron collisions_ , _Phys. Rev. D_ 46 (1992) 4930.
  * (12) P. J. Rijken and W. L. van Neerven, _Heavy flavor contributions to the Drell-Yan cross-section_ , _Phys. Rev. D_ 52 (1995) 149 [hep-ph/9501373].
  * (13) C. Anastasiou, L. J. Dixon, K. Melnikov and F. Petriello, _Dilepton rapidity distribution in the Drell-Yan process at NNLO in QCD_ , _Phys. Rev. Lett._ 91 (2003) 182002 [hep-ph/0306192].
  * (14) L. Cieri, F. Coradeschi and D. de Florian, _Diphoton production at hadron colliders: transverse-momentum resummation at next-to-next-to-leading logarithmic accuracy_ , _JHEP_ 06 (2015) 185 [1505.03162].
  * (15) R. Gavin, Y. Li, F. Petriello and S. Quackenbush, _FEWZ 2.0: A code for hadronic Z production at next-to-next-to-leading order_ , _Comput. Phys. Commun._ 182 (2011) 2388 [1011.3540].
  * (16) R. Boughezal, J. M. Campbell, R. K. Ellis, C. Focke, W. Giele, X. Liu et al., _Color singlet production at NNLO in MCFM_ , _Eur. Phys. J. C_ 77 (2017) 7 [1605.08011].
  * (17) M. Grazzini, S. Kallweit and M. Wiesemann, _Fully differential NNLO computations with MATRIX_ , _Eur. Phys. J. C_ 78 (2018) 537 [1711.06631].
  * (18) M. Delto, M. Jaquier, K. Melnikov and R. Roentsch, _Mixed QCD $\otimes$QED corrections to on-shell $Z$ boson production at the LHC_, _JHEP_ 01 (2020) 043 [1909.08428].
  * (19) R. Bonciani, F. Buccioni, N. Rana, I. Triscari and A. Vicini, _NNLO QCD $\times$EW corrections to Z production in the $q\bar{q}$ channel_, _Phys. Rev._ D101 (2020) 031301 [1911.06200].
  * (20) R. Bonciani, F. Buccioni, N. Rana and A. Vicini, _Next-to-Next-to-Leading Order Mixed QCD-Electroweak Corrections to on-Shell Z Production_ , _Phys. Rev. Lett._ 125 (2020) 232004 [2007.06518].
  * (21) R. Bonciani, L. Buonocore, M. Grazzini, S. Kallweit, N. Rana, F. Tramontano et al., _Mixed strong $-$electroweak corrections to the Drell$-$Yan process_, 2106.11953.
  * (22) F. Buccioni, F. Caola, M. Delto, M. Jaquier, K. Melnikov and R. Röntsch, _Mixed QCD-electroweak corrections to on-shell Z production at the LHC_ , _Phys. Lett. B_ 811 (2020) 135969 [2005.10221].
  * (23) C. Duhr, F. Dulat and B. Mistlberger, _The Drell-Yan cross section to third order in the strong coupling constant_ , _Phys. Rev. Lett._ 125 (2020) 172001 [2001.07717].
  * (24) C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger, _Higgs Boson Gluon-Fusion Production in QCD at Three Loops_ , _Phys. Rev. Lett._ 114 (2015) 212001 [1503.06056].
  * (25) C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog et al., _High precision determination of the gluon fusion Higgs boson cross-section at the LHC_ , _JHEP_ 05 (2016) 058 [1602.00695].
  * (26) C. Duhr, F. Dulat and B. Mistlberger, _Higgs Boson Production in Bottom-Quark Fusion to Third Order in the Strong Coupling_ , _Phys. Rev. Lett._ 125 (2020) 051804 [1904.09990].
  * (27) L.-B. Chen, H. T. Li, H.-S. Shao and J. Wang, _Higgs boson pair production via gluon fusion at N 3LO in QCD_, _Phys. Lett._ B803 (2020) 135292 [1909.06808].
  * (28) C. Duhr, F. Dulat, V. Hirschi and B. Mistlberger, _Higgs production in bottom quark fusion: matching the 4- and 5-flavour schemes to third order in the strong coupling_ , _JHEP_ 08 (2020) 017 [2004.04752].
  * (29) F. A. Dreyer and A. Karlberg, _Vector-Boson Fusion Higgs Pair Production at N 3LO_, _Phys. Rev._ D98 (2018) 114016 [1811.07906].
  * (30) F. A. Dreyer and A. Karlberg, _Vector-Boson Fusion Higgs Production at Three Loops in QCD_ , _Phys. Rev. Lett._ 117 (2016) 072001 [1606.00840].
  * (31) X. Chen, T. Gehrmann, N. Glover, A. Huss, T.-Z. Yang and H. X. Zhu, _Di-lepton Rapidity Distribution in Drell-Yan Production to Third Order in QCD_ , 2107.09085.
  * (32) S. Camarda, L. Cieri and G. Ferrera, _Drell-Yan lepton-pair production: $q_{T}$ resummation at N3LL accuracy and fiducial cross sections at N3LO_, 2103.04974.
  * (33) C. Duhr, F. Dulat and B. Mistlberger, _Charged current Drell-Yan production at N 3LO_, _JHEP_ 11 (2020) 143 [2007.13313].
  * (34) C. Gnendiger et al., _To ${d}$, or not to ${d}$: recent developments and comparisons of regularization schemes_, _Eur. Phys. J. C_ 77 (2017) 471 [1705.01827].
  * (35) S. L. Adler, _Axial vector vertex in spinor electrodynamics_ , _Phys. Rev._ 177 (1969) 2426.
  * (36) J. S. Bell and R. Jackiw, _A PCAC puzzle: $\pi^{0}\to\gamma\gamma$ in the $\sigma$ model_, _Nuovo Cim. A_ 60 (1969) 47.
  * (37) S. L. Adler and W. A. Bardeen, _Absence of higher order corrections in the anomalous axial vector divergence equation_ , _Phys. Rev._ 182 (1969) 1517.
  * (38) C.-N. Yang, _Selection Rules for the Dematerialization of a Particle Into Two Photons_ , _Phys. Rev._ 77 (1950) 242.
  * (39) J. C. Collins, F. Wilczek and A. Zee, _Low-Energy Manifestations of Heavy Particles: Application to the Neutral Current_ , _Phys. Rev. D_ 18 (1978) 242.
  * (40) G. ’t Hooft and M. J. G. Veltman, _Regularization and Renormalization of Gauge Fields_ , _Nucl. Phys. B_ 44 (1972) 189.
  * (41) P. Breitenlohner and D. Maison, _Dimensional Renormalization and the Action Principle_ , _Commun. Math. Phys._ 52 (1977) 11.
  * (42) S. A. Larin and J. A. M. Vermaseren, _The alpha-s**3 corrections to the Bjorken sum rule for polarized electroproduction and to the Gross-Llewellyn Smith sum rule_ , _Phys. Lett. B_ 259 (1991) 345.
  * (43) S. A. Larin, _The Renormalization of the axial anomaly in dimensional regularization_ , _Phys. Lett. B_ 303 (1993) 113 [hep-ph/9302240].
  * (44) F. Jegerlehner, _Facts of life with gamma(5)_ , _Eur. Phys. J. C_ 18 (2001) 673 [hep-th/0005255].
  * (45) R. V. Harlander and W. B. Kilgore, _Next-to-next-to-leading order Higgs production at hadron colliders_ , _Phys. Rev. Lett._ 88 (2002) 201801 [hep-ph/0201206].
  * (46) B. Mistlberger, _Higgs boson production at hadron colliders at N 3LO in QCD_, _JHEP_ 05 (2018) 028 [1802.00833].
  * (47) P. Nogueira, _Automatic Feynman Graph Generation_ , _Journal of Computational Physics_ 105 (1993) 279.
  * (48) C. W. Bauer, R. Kreckel and A. Frink, _Introduction to the GiNaC framework for symbolic computation within the C++ programming language_ , _J.Symb.Comput._ 33 (2000) 1.
  * (49) B. Ruijl, T. Ueda and J. Vermaseren, _FORM version 4.2_ , 1707.06453.
  * (50) F. Tkachov, _A theorem on analytical calculability of 4-loop renormalization group functions_ , _Physics Letters B_ 100 (1981) 65.
  * (51) K. Chetyrkin and F. Tkachov, _Integration by parts: The algorithm to calculate $\beta$-functions in 4 loops_, _Nuclear Physics B_ 192 (1981) 159.
  * (52) S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys._ A15 (2000) 5087 [hep-ph/0102033].
  * (53) C. Anastasiou and K. Melnikov, _Higgs boson production at hadron colliders in NNLO QCD_ , _Nuclear Physics B_ 646 (2002) 220 [0207004].
  * (54) C. Anastasiou and K. Melnikov, _Pseudoscalar Higgs boson production at hadron colliders in next-to-next-to-leading order QCD_ , _Physical Review D_ 67 (2003) 037501.
  * (55) C. Anastasiou, L. Dixon and K. Melnikov, _NLO Higgs boson rapidity distributions at hadron colliders_ , _Nuclear Physics B \- Proceedings Supplements_ 116 (2003) 193.
  * (56) C. Anastasiou, L. Dixon, K. Melnikov and F. Petriello, _High-precision QCD at hadron colliders: Electroweak gauge boson rapidity distributions at next-to-next-to leading order_ , _Physical Review D_ 69 (2004) 094008.
  * (57) A. Kotikov, _Differential equations method: New technique for massive Feynman diagrams calculation_ , _Phys. Lett. B_ 254 (1991) 158.
  * (58) A. Kotikov, _Differential equations method: The Calculation of vertex type Feynman diagrams_ , _Phys. Lett. B_ 259 (1991) 314.
  * (59) A. Kotikov, _Differential equation method: The Calculation of N point Feynman diagrams_ , _Phys. Lett. B_ 267 (1991) 123.
  * (60) J. M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601 [1304.1806].
  * (61) T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys._ B580 (2000) 485 [hep-ph/9912329].
  * (62) C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger, _Real-virtual contributions to the inclusive Higgs cross-section at $N^{3}LO$_, _JHEP_ 12 (2013) 088 [1311.1425].
  * (63) C. Anastasiou, C. Duhr, F. Dulat and B. Mistlberger, _Soft triple-real radiation for Higgs production at N3LO_ , _JHEP_ 07 (2013) 003 [1302.4379].
  * (64) C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, F. Herzog and B. Mistlberger, _Soft expansion of double-real-virtual corrections to Higgs production at N 3LO_, _JHEP_ 08 (2015) 051 [1505.04110].
  * (65) F. Dulat and B. Mistlberger, _Real-Virtual-Virtual contributions to the inclusive Higgs cross section at N3LO_ , 1411.3586.
  * (66) V. N. Gribov and L. N. Lipatov, _Deep inelastic e p scattering in perturbation theory_ , _Sov. J. Nucl. Phys._ 15 (1972) 438\.
  * (67) G. Altarelli and G. Parisi, _Asymptotic Freedom in Parton Language_ , _Nucl. Phys._ B126 (1977) 298.
  * (68) 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.
  * (69) S. Moch, J. A. M. Vermaseren and A. Vogt, _The Three loop splitting functions in QCD: The Nonsinglet case_ , _Nucl. Phys._ B688 (2004) 101 [hep-ph/0403192].
  * (70) A. Vogt, S. Moch and J. A. M. Vermaseren, _The Three-loop splitting functions in QCD: The Singlet case_ , _Nucl. Phys._ B691 (2004) 129 [hep-ph/0404111].
  * (71) J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, _The 3-loop pure singlet heavy flavor contributions to the structure function $F_{2}(x,Q^{2})$ and the anomalous dimension_, _Nucl. Phys. B_ 890 (2014) 48 [1409.1135].
  * (72) J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, _The three-loop splitting functions $P_{qg}^{(2)}$ and $P_{gg}^{(2,N_{F})}$_, _Nucl. Phys._ B922 (2017) 1 [1705.01508].
  * (73) M. Bos, _Explicit calculation of the renormalized singlet axial anomaly_ , _Nucl. Phys. B_ 404 (1993) 215 [hep-ph/9211319].
  * (74) T. Ahmed, T. Gehrmann, P. Mathews, N. Rana and V. Ravindran, _Pseudo-scalar Form Factors at Three Loops in QCD_ , _JHEP_ 11 (2015) 169 [1510.01715].
  * (75) T. Ahmed, L. Chen and M. Czakon, _Renormalization of the flavor-singlet axial-vector current and its anomaly in dimensional regularization_ , _JHEP_ 05 (2021) 087 [2101.09479].
  * (76) L. Chen, M. Czakon and M. Niggetiedt, _The complete singlet contribution to the massless quark form factor at three loops in QCD_ , 2109.01917.
  * (77) T. Gehrmann and A. Primo, _The three-loop singlet contribution to the massless axial-vector quark form factor_ , _Phys. Lett. B_ 816 (2021) 136223 [2102.12880].
  * (78) T. Gehrmann, E. W. N. Glover, T. Huber, N. Ikizlerli and C. Studerus, _Calculation of the quark and gluon form factors to three loops in QCD_ , _JHEP_ 06 (2010) 094 [1004.3653].
  * (79) T. Gehrmann, E. W. N. Glover, T. Huber, N. Ikizlerli and C. Studerus, _The quark and gluon form factors to three loops in QCD through to O(eps 2)_, _JHEP_ 11 (2010) 102 [1010.4478].
  * (80) W.-L. Ju and M. Schönherr, _The q T and $\Delta$$\phi$ spectra in W and Z production at the LHC at N3LL’+N2LO_, _JHEP_ 10 (2021) 088 [2106.11260].
  * (81) T. van Ritbergen, J. A. M. Vermaseren and S. A. Larin, _The Four loop beta function in quantum chromodynamics_ , _Phys. Lett._ B400 (1997) 379 [hep-ph/9701390].
  * (82) M. Czakon, _The Four-loop QCD beta-function and anomalous dimensions_ , _Nucl. Phys._ B710 (2005) 485 [hep-ph/0411261].
  * (83) P. A. Baikov, K. G. Chetyrkin and J. H. Kühn, _Five-Loop Running of the QCD coupling constant_ , _Phys. Rev. Lett._ 118 (2017) 082002 [1606.08659].
  * (84) F. Herzog, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt, _The five-loop beta function of Yang-Mills theory with fermions_ , _JHEP_ 02 (2017) 090 [1701.01404].
  * (85) J. Butterworth et al., _PDF4LHC recommendations for LHC Run II_ , _J. Phys._ G43 (2016) 023001 [1510.03865].
  * (86) S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky et al., _New parton distribution functions from a global analysis of quantum chromodynamics_ , _Phys. Rev._ D93 (2016) 033006 [1506.07443].
  * (87) L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne, _Parton distributions in the LHC era: MMHT 2014 PDFs_ , _Eur. Phys. J._ C75 (2015) 204 [1412.3989].
  * (88) NNPDF collaboration, R. D. Ball et al., _Parton distributions for the LHC Run II_ , _JHEP_ 04 (2015) 040 [1410.8849].
  * (89) S. Alekhin, J. Bluemlein, S.-O. Moch and R. Placakyte, _The new ABMP16 PDF_ , _PoS_ DIS2016 (2016) 016 [1609.03327].
  * (90) T.-J. Hou et al., _New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC_ , _Phys. Rev. D_ 103 (2021) 014013 [1912.10053].
  * (91) S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin and R. S. Thorne, _Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs_ , _Eur. Phys. J. C_ 81 (2021) 341 [2012.04684].
  * (92) NNPDF collaboration, R. D. Ball et al., _Parton distributions from high-precision collider data_ , _Eur. Phys. J._ C77 (2017) 663 [1706.00428].
  * (93) R. D. Ball et al., _The Path to Proton Structure at One-Percent Accuracy_ , 2109.02653.
  * (94) PDF4LHC21 combination group collaboration, T. Cridge, _PDF4LHC21: Update on the benchmarking of the CT, MSHT and NNPDF global PDF fits_ , in _28th International Workshop on Deep Inelastic Scattering and Related Subjects_ , 8, 2021, 2108.09099.
  * (95) L. J. Dixon and A. Signer, _Complete O (alpha-s**3) results for e+ e- — $>$ (gamma, Z) —$>$ four jets_, _Phys. Rev. D_ 56 (1997) 4031 [hep-ph/9706285].


# Dynamic Restrained Uncertainty Weighting Loss
for Multitask Learning of Vocal Expression

Meishu Song  Zijiang Yang  Andreas Triantafyllopoulos  Xin Jing  Vincent Karas
Xie Jiangjian  Zixing Zhang  Yamamoto Yoshiharu  Björn W. Schuller

###### Abstract

We propose a novel Dynamic Restrained Uncertainty Weighting Loss to
experimentally handle the problem of balancing the contributions of multiple
tasks on the ICML ExVo 2022 Challenge. The multitask aims to recognise
expressed emotions and demographic traits from vocal bursts jointly. Our
strategy combines the advantages of Uncertainty Weight and Dynamic Weight
Average, by extending weights with a restraint term to make the learning
process more explainable. We use a lightweight multi-exit CNN architecture to
implement our proposed loss approach. The experimental H-Mean score (0.394)
shows a substantial improvement over the baseline H-Mean score (0.335).

multitask learning, CNN

##  1 Introduction

“Transfer should always be useful”; any pair of distributions underlying a
pair of tasks must have something in common (Mahmud, 2009). Many deep learning
works attempt to improve their performance using Multitask Learning (MTL) – a
form of learning which constitutes the learning of several related tasks
simultaneously. MTL has been applied successfully in a variety of fields,
including emotion recognition (Eyben et al., 2012), (Chen et al., 2017), (Deng
et al., 2020), (Shen et al., 2021), visual scene understanding (Liu et al.,
2019), (Kendall et al., 2018) and automatic speech recognition (Krishna et
al., 2018). The multitasking learning strategy has several advantages compared
with single task learning. It allows the model to discover general
representations that are useful for multiple tasks, giving it greater
contextual knowledge and improving performance in the individual tasks
compared to single-task learning (Kendall et al., 2018). In addition, the
shared parts of the network architecture make this approach more
computationally efficient than training a separate model for each task, which
is important in settings where the model has to run in real time (Kendall et
al., 2018).

Even though multitasking has demonstrable theoretical advantages, there are
two main issues in applying multitasking to real works (Liu et al., 2019). 1\.
Define what to share. Traditional MTL algorithms usually assume that all the
tasks are related, which means tasks are sharing the same embeddings. However,
the same low-level features or high-level features cannot always perform well
for all tasks. 2\. Define how to balance tasks. Prior MTL works treated all
tasks with equal weights, which leads to easier tasks overfitting while tough
tasks are still training. In other words, model performance is extremely
sensitive to weight selection for balancing the individual losses.

Various prior works have attempted to address these issues. One existing
approach is using multi-exit architectures for different objectives (Phuong &
Lampert, 2019), in which a stack of processing output layers is interleaved
with early output layers. The work (Liu et al., 2019) introduced a CNN-based
framework that involves multiple levels of a shared network, in order to learn
a global feature set and the contributions of features with varying complexity
for each specific task. Both of these methods can help determine the most
appropriate output layer for different tasks. In this regard, we propose a
multi-exit CNN to solve the first issue.

Loss weighting strategies in MTL can be classified based on task gradients and
loss: (1) pure loss strategies including Uncertainty Weighting (UW) (Kendall
et al., 2018) and Dynamic Weight Average (DWA) (Liu et al., 2019), and (2)
loss based on gradients including Gradient Normalisation (Chen et al., 2018),
Multiple-Gradient Descent Algorithm (Sener & Koltun, 2018), Projecting
Conflicting Gradient (Yu et al., 2020), Gradient Sign Dropout (Chen et al.,
2020), Impartial Multi-Task Learning (Liu et al., 2021), and Gradient Vaccine
(Wang et al., 2020). In this work, we limit our analysis to pure loss
strategies since conflicting gradient signals coming from different tasks can
degrade model performance (Gong et al., 2019). In comparison, loss strategies,
such as UW and DWA, are easy to be implemented, and continue to be popular in
the literature (Gong et al., 2019). We aim to combine the advantages of UW and
DWA by extending UW with a term restraining the loss weights, which improves
performance while also making the learning process more explainable.

Our contributions are two-fold: We use a multi-exit CNN to learn data
representations at multiple layers, resulting in a simple, effective, and
lightweight model that achieves competitive performance. We propose a novel
loss weighting strategy, namely the Dynamic Restrained Uncertainty Weighting
loss for balancing multiple tasks, which incorporates both UW’s and DWA’s
benefits.

##  2 Dataset: ExVo Multitask Challenge

The dataset111We use the original release of processed .wav files not the
follow-up release of unprocessed ones. used for the ExVo Multitask challenge
is the HUME-VB competition dataset (Baird et al., 2022), a large-scale dataset
of emotional human vocal bursts, which was collected in four countries with
broadly different cultures: China, South Africa, USA, and Venezuela. There are
59,299 recordings totalling 36:50:40 (HH:MM:SS) of audio data from 1702
speakers, aged from 20 to 39 years old. Furthermore, the data was recorded in
speakers’ homes with their own microphones. Each vocal burst had been self-
annotated by each speaker in terms of the intensity of ten different expressed
emotions, each on a 1-100 scale, Amusement, Awe, Awkwardness, Distress,
Excitement, Fear, Horror, Sadness, Surprise, and Triumph. In this paper, we
utilise this dataset to train a multi-task model to jointly predict the
expression of ten emotions along with the age and native-country of the
speaker.

##  3 Methodology

We now introduce our multi-exit CNN for MTL, which addresses the “what to
share” issue and Dynamic Restraint Uncertainty Weighting (DRUW) loss to
address the “how to balance task weights” issue.

###  3.1 Multi-Exit CNN Architecture

Figure 1: Overview of the multi-exit CNN architecture.

CNN architectures are naturally hierarchical (Zhu & Bain, 2017). In (Zeiler &
Fergus, 2014), it is shown that lower layers in CNNs usually capture the low-
level features, while higher layers are likely to extract high-level features.
As a consequence, a possible way to embed different tasks into a single CNN
model is to obtain predictions from different CNN layers as data flow through
it (Zhu & Bain, 2017), thus effectively transforming a single, monolithic
architecture to a scalable, multi-exit one. Inspired by the flexibility of
controlling output layers of a CNN, we propose a multi-exit CNN model by
adding different output layers to learn different data representations. The
overview of the architecture is shown in Figure 1. Our network is composed of
five hierarchical CNN blocks which contains a total of ten CNN layers, and
each block has two layers. The Age prediction is branched off after the first
Conv block, the Country prediction after the third Conv block, while the
emotion prediction after the fifth Conv block. One important aspect of this
work is determining the depth of the exit branch for each task. This is done
experimentally by selecting the depth which gives the best performance. Due to
the layer number being relatively low, we apply a grid search technique based
on the validation set performance to find the proper exits for each task.

###  3.2 Loss Weighting Strategies

####  3.2.1 Uncertainty Weighting Loss

The performance of MTL models is largely determined by shared weights,
according to (Kendall et al., 2018), but these weights are difficult to train.
This leads to the concept of uncertainty to measure the loss of different
tasks, making it possible to learn different types of tasks simultaneously
(Kendall et al., 2018). UW is defined as in Equation 1:

$\displaystyle\mathcal{L}(w,\alpha_{1},\alpha_{2})$
$\displaystyle=\frac{1}{\alpha_{1}^{2}}\mathcal{L}_{1}(w)+\frac{1}{\alpha_{2}^{2}}\mathcal{L}_{2}(w)$
(1) $\displaystyle\quad+log\alpha_{1}+log\alpha_{2}.$

The authors interpret $\alpha_{1}$ and $\alpha_{2}$ as learning the relative
weight of the losses $\mathcal{L}_{1}(w)$ and $\mathcal{L}_{2}(w)$ adaptively
(Kendall et al., 2018), based on the data. For a task with a larger loss and
higher uncertainty, the strategy can effectively avoid the model to “take a
big step” towards a “blind” direction. Thus, it should update $w$ with a
smaller gradient. At the same time, the strategy can avoid the problem of a
task with larger loss dominating the overall loss. Recently, the Revised
Uncertainty Weighting (RUW) Loss (Liebel & Körner, 2018) was proposed to solve
a major disadvantage of the UW loss: experimentally, because the log part in
Equation 1 leads the whole loss to a negative value (Liebel & Körner, 2018),
the RUW adapts the regularisation term $log\alpha$ to $log(1+log\alpha^{2})$
in order to enforce positive values. The concept can be represented as in
Equation 2:

$\displaystyle\mathcal{L}(w,\alpha_{1},\alpha_{2})$
$\displaystyle=\frac{1}{\alpha_{1}^{2}}\mathcal{L}_{1}(w)+\frac{1}{\alpha_{2}^{2}}\mathcal{L}_{2}(w)$
(2) $\displaystyle\quad+log(1+log\alpha_{1}^{2})+log(1+log\alpha_{2}^{2}).$

####  3.2.2 Dynamic Weight Average

In (Liu et al., 2019), Dynamic Weight Average (DWA) is proposed as a simple,
yet effective adaptive weighting method. Inspired by GradNorm (Chen et al.,
2018), this learns to average task weighting over time by evaluating the rate
of change of individual loss. In detail, a weighting $\lambda_{k}$ is defined
for each task $k$ as:

$\lambda_{k}(t)=\frac{\mathcal{K}exp(\frac{\mathcal{L}_{k}(t-1)}{\mathcal{L}_{k}(t-2)}/\mathcal{T})}{\sum_{i}exp(\frac{\mathcal{L}_{k}(t-1)}{\mathcal{L}_{k}(t-2)}/\mathcal{T})},$
(3)

where $t$ is an iteration index and $w_{k}$ represents the relative descending
rate in the range (0, +$\infty$). First, $\mathcal{T}$ plays the role of
smoothing task weights; the larger $\mathcal{T}$ is, the more evenly
distributed the weights of different tasks are. Finally, the softmax operator
multiplies by $\mathcal{K}$ to ensure that $\sum_{i}\lambda_{i}(t)=K.$ In this
setting, the total loss value is the sum of DWA weights multiply each task
loss.

Table 1: ExVo-MultiTask task validation and test results. Mean concordance correlation coefficient (CCC) across the ten (Emo)tional classes, unweighted average recall (UAR) for the four class (Cou)try task (chance level 0.25), and mean absolute error (MAE) for the Age regression task are all reported. H-Mean score means harmonic mean between these metrics. Emphasised results indicate best scores. Besides official baselines, we report Single Task (ST) baselines and MTL baselines using our proposed multi-exit CNN model. For comparison, we report weighting strategies: Equal Weighting (EW), Uncertainty Weighting (UW), Revised Uncertainty Weighting (RUW), our Restrained Uncertainty Weighting (RRUW), Dynamic Weight Average (DWA), and our Dynamic Restrained Uncertainty Weighting (DRUW). |  | Validation | Test
---|---|---|---
| Weighting | Emo-CCC | Cou-UAR | Age-MAE | H-Mean |  | H-Mean
Official Baselines ComPARE | EW | 0.416 | 0.506 | 4.222 | 0.349 |  | 0.335
Official Baselines eGeMAPS | EW | 0.353 | 0.423 | 4.011 | 0.324 |  | 0.314
Our CNN ST Baselines | — | 0.645 | 0.588 | 3.926 | 0.418 |  | 0.393
Our MTL Baselines | EW | 0.633 | 0.525 | 3.928 | 0.405 |  | 0.382
Our MTL | UW | 0.615 | 0.575 | 4.024 | 0.406 |  | 0.391
RUW | 0.629 | 0.539 | 3.798 | 0.414 |  | 0.391
Our RRUW | 0.635 | 0.576 | 3.803 | 0.421 |  | 0.392
DWA | 0.637 | 0.545 | 3.754 | 0.419 |  | 0.389
Our DRUW | 0.635 | 0.570 | 3.763 | 0.423 |  | 0.394

####  3.2.3 Proposed Loss Strategy

We firstly use all previously described approaches for our task. We notice
that both, the UW and RUW methods, lead to a similar outcome: since the
weights are also trainable, if all weights turn to decrease simultaneously,
the whole loss value is also decreasing, thus leading to a trivial solution.
Thus, we propose to extend the UW and RUW with a constraint. Similarly with
DWA, it uses a “softmax” output that converts the logits of the loss change
ratio into a weight for each task. Also, in the original paper of UW (Kendall
et al., 2018), the UW formula can be “scaled up”. In more detail, using
softmax, all losses can be summed up to one. “Scaled” is a term used to
describe a situation where all losses can be added up to a customised positive
value. Differently, this is suggested in a classification task (Kendall et
al., 2018). In our case, the task is a heterogeneous multitask problem which
involves regression and classification tasks in parallel. To constrain the
weights in a controlled manner, we use a different strategy with softmax, and
we simply add the following $\varphi$ with RUW, namely Restrained Revised
Uncertainty Weighting (RRUW) Loss:

$\displaystyle\mathcal{L}(w,\alpha_{1},\alpha_{2})$
$\displaystyle=\frac{1}{\alpha_{1}^{2}}\mathcal{L}_{1}(w)+\frac{1}{\alpha_{2}^{2}}\mathcal{L}_{2}(w)$
(4) $\displaystyle\quad+log(1+log\alpha_{1}^{2})+log(1+log\alpha_{2}^{2})$
$\displaystyle\quad+\left|\varphi-(\left|log\alpha_{1}\right|+\left|log\alpha_{2}\right|)\right|.$

Our method regularises the weights by guiding their sum towards a fixed
positive value $\varphi$, which can constrain the sum of all weights so that
we avoid the problem of the weights degenerating to a trivial solution. To
maintain a positive final loss, we additionally use the absolute value of
$\varphi$ for the loss. Unlike softmax, which adds up all weights into one,
our method controls the weight sum into a ‘flexible’ positive value of
$\varphi$. The term ‘flexible’ refers to the fact that the sum of weights can
be indefinitely near to $\varphi$, the fixed value. Experimentally, DWA is
shown to effectively decrease over-fitting effects of overly complex networks
for easy tasks, while RRUW can uncover good data representations for complex
tasks. This leads to the question: “how can we best combine both strategies?”.
We extend our idea to get sum values of RRUW and DWA weights for each task,
which means for each training, there are two weight suggestion values. The
joint loss is given as:

$\displaystyle\mathcal{L}(w,\alpha_{1},\alpha_{2})$
$\displaystyle=\frac{\mathcal{K}exp(\frac{\mathcal{L}_{1}(t-1)}{\mathcal{L}_{1}(t-2)}/\mathcal{T})}{\sum_{i}exp(\frac{\mathcal{L}_{1}(t-1)}{\mathcal{L}_{1}(t-2)}/\mathcal{T})}\mathcal{L}_{1}+\frac{1}{\alpha_{1}^{2}}\mathcal{L}_{1}(w)$
(5)
$\displaystyle\quad+\frac{\mathcal{K}exp(\frac{\mathcal{L}_{2}(t-1)}{\mathcal{L}_{2}(t-2)}/\mathcal{T})}{\sum_{i}exp(\frac{\mathcal{L}_{2}(t-1)}{\mathcal{L}_{2}(t-2)}/\mathcal{T})}\mathcal{L}_{2}+\frac{1}{\alpha_{2}^{2}}\mathcal{L}_{2}(w)$
$\displaystyle\quad+log(1+log\alpha_{1}^{2})+log(1+log\alpha_{2}^{2})$
$\displaystyle\quad+\left|\varphi-(\left|log\alpha_{1}\right|+\left|log\alpha_{2}\right|)\right|.$

###  3.3 Experimental Set-Up

We extracted Mel-Spectrogram features by randomly cutting raw audios into 2.5
seconds length. The size of the Mel-Spectrogram is: 64*512. In the training
session, we use AdamW as our optimiser. The initial learning rate is 0.001;
the batch size is 32. All of our models are trained on Nvidia RTX3090 and
Nvidia A40 GPUs, and the number of training epochs are limited to 60 on all
tasks. In loss setting, the $\varphi$ value is 1 and $\mathcal{T}$ is 10. In
our implementation, all models are trained without data augmentation and no
ensemble methods are adopted.

##  4 Results

For the baseline results (Baird et al., 2022) in Table 1, we present four
types including two official baselines – Single Task (ST) and MTL baselines –,
along with our proposed multi-exit CNN networks. On single Emotions and the
Country task (CCC: 0.645 and UAR: 0.588), they outperform MTL (CCC: 0.635 and
UAR: 0.570), supported by previous work (Zhang et al., 2019), MTL does not
always guarantee a better performance. However, on the Age task, MTL (MAE:
0.3763) is superior to the one on ST (MAE: 3.928). Furthermore, MTL
(Validation: 0.423, Test: 0.394) presents better H-Mean scores compared with
ST (Validation: 0.418, Test: 0.393). Overall, our results confirm the
advantage of MTL over ST – it can build well-shared data representations with
limited training resources (Gong et al., 2019).

We also present MTL results with various loss weighting schemes. First,
comparing to EW, all the strategies show improvements on H-Mean scores. This
suggests that the two main stream losses, as well as our proposed losses, can
all be effective in balancing task weights during training. Second, it is
noted that our proposed multi-exit CNN with RRUW weighting strategy achieves
the better H-Mean scores (Validation: 0.421 and Test: 0.392) than the previous
strategies (UW, RUW, and DWA). A rational behind this is that our “Restraint”
component can assist the model to balance the weight in a regulated manner and
provide easily explainable features of the model. In addition, we observe
that, when integrating with the DRUW (Validation H-Mean: 0.423 and Test
H-Mean: 0.394) strategy, the system performance results in better H-Mean
scores than the RRUW strategy (Validation H-Mean: 0.421 and Test H-Mean:
0.392). This indicates that as we expected, the DRUW method can keep task
weights well balanced and combine the advantages of UW and DWA, as well as
keep UW in a controlled way due to the “Restraint” component and generate an
easily explainable data representation.

##  5 Conclusion

In this work, we contributed a novel loss weighting strategy: DRUW, which
applies on a light and effective multi-exit CNN network. Our proposed method
combines the previous loss strategies: UW and DWA, and also gained a
restrained feature for a clear explainable purpose. The experimental results
demonstrate a considerable improvement over the baseline results, which
indicates the effectiveness of our methods on the ExVo Multitask challenge.
Future work includes automatic adjustment of CNN multi exits on each task. In
addition, we also intend to explore the comparison between gradient oriented
loss strategies and our loss strategies. Furthermore, the automatic balancing
of the weights of RRUW and DWA is also promising in the broader MTL field.

## References

  * Baird et al. (2022) Baird, A., Tzirakis, P., Gidel, G., Jiralerspong, M., Muller, E. B., Mathewson, K., Schuller, B., Cambria, E., Keltner, D., and Cowen, A.  The icml 2022 expressive vocalizations workshop and competition: Recognizing, generating, and personalizing vocal bursts.  _arXiv preprint arXiv:2205.01780_ , 2022.

  * Chen et al. (2017) Chen, S., Jin, Q., Zhao, J., and Wang, S.  Multimodal multi-task learning for dimensional and continuous emotion recognition.  In _Proc. Annual Workshop on Audio/Visual Emotion Challenge_ , pp. 19–26, New York, NY, USA, 2017.
  * Chen et al. (2018) Chen, Z., Badrinarayanan, V., Lee, C.-Y., and Rabinovich, A.  Gradnorm: Gradient normalization for adaptive loss balancing in deep multitask networks.  In _Proc. International Conference on Machine Learning_ , pp. 794–803, Stockholm, Sweden, 2018.
  * Chen et al. (2020) Chen, Z., Ngiam, J., Huang, Y., Luong, T., Kretzschmar, H., Chai, Y., and Anguelov, D.  Just pick a sign: Optimizing deep multitask models with gradient sign dropout.  _Advances in Neural Information Processing Systems_ , 33:2039–2050, 2020.
  * Deng et al. (2020) Deng, D., Chen, Z., and Shi, B. E.  Multitask emotion recognition with incomplete labels.  In _Proc. Automatic Face and Gesture Recognition_ , pp. 592–599, Virtual, 2020.
  * Eyben et al. (2012) Eyben, F., Wöllmer, M., and Schuller, B.  A multitask approach to continuous five-dimensional affect sensing in natural speech.  _ACM Trans. Interact. Intell. Syst._ , 2(1), March 2012\.  ISSN 2160-6455.
  * Gong et al. (2019) Gong, T., Lee, T., Stephenson, C., Renduchintala, V., Padhy, S., Ndirango, A., Keskin, G., and Elibol, O. H.  A comparison of loss weighting strategies for multi task learning in deep neural networks.  _IEEE Access_ , 7:141627–141632, 2019.
  * Kendall et al. (2018) Kendall, A., Gal, Y., and Cipolla, R.  Multi-task learning using uncertainty to weigh losses for scene geometry and semantics.  In _Proc. Computer Vision and Pattern Recognition_ , pp. 7482–7491, Salt Lake City, UT, USA, 2018.
  * Krishna et al. (2018) Krishna, K., Toshniwal, S., and Livescu, K.  Hierarchical multitask learning for ctc-based speech recognition.  _CoRR_ , abs/1807.06234, 2018.  URL http://arxiv.org/abs/1807.06234.
  * Liebel & Körner (2018) Liebel, L. and Körner, M.  Auxiliary tasks in multi-task learning.  _arXiv preprint arXiv:1805.06334_ , 2018.
  * Liu et al. (2021) Liu, L., Li, Y., Kuang, Z., Xue, J., Chen, Y., Yang, W., Liao, Q., and Zhang, W.  Towards impartial multi-task learning.  In _Proc. International Conference on Learning Representations_ , Vienna, Austria, 2021.
  * Liu et al. (2019) Liu, S., Johns, E., and Davison, A. J.  End-to-end multi-task learning with attention.  In _Proc. Computer Vision and Pattern Recognition_ , Long Beach, CA, USA, 2019.
  * Mahmud (2009) Mahmud, M. H.  On universal transfer learning.  _Theoretical Computer Science_ , 410(19):1826–1846, 2009.
  * Phuong & Lampert (2019) Phuong, M. and Lampert, C. H.  Distillation-based training for multi-exit architectures.  In _Proc. International Conference on Computer Vision_ , pp. 1355–1364, Seoul, South Korea, 2019.
  * Sener & Koltun (2018) Sener, O. and Koltun, V.  Multi-task learning as multi-objective optimization.  2018\.
  * Shen et al. (2021) Shen, J., Zheng, J., and Wang, X.  Mmtrans-mt: A framework for multimodal emotion recognition using multitask learning.  In _Proc. International Conference on Advanced Computational Intelligence_ , pp. 52–59, Doha, Qatar, 2021.
  * Wang et al. (2020) Wang, Z., Tsvetkov, Y., Firat, O., and Cao, Y.  Gradient vaccine: Investigating and improving multi-task optimization in massively multilingual models.  _arXiv preprint arXiv:2010.05874_ , 2020.
  * Yu et al. (2020) Yu, T., Kumar, S., Gupta, A., Levine, S., Hausman, K., and Finn, C.  Gradient surgery for multi-task learning.  _Advances in Neural Information Processing Systems_ , 33:5824–5836, 2020.
  * Zeiler & Fergus (2014) Zeiler, M. D. and Fergus, R.  Visualizing and understanding convolutional networks.  In _Proc. European conference on computer vision_ , pp. 818–833. Springer, 2014.
  * Zhang et al. (2019) Zhang, Z., Wu, B., and Schuller, B.  Attention-augmented end-to-end multi-task learning for emotion prediction from speech.  In _Proc. International Conference on Acoustics, Speech and Signal Processing_ , pp. 6705–6709, Brighton, UK, 2019.
  * Zhu & Bain (2017) Zhu, X. and Bain, M.  B-cnn: branch convolutional neural network for hierarchical classification.  _arXiv preprint arXiv:1709.09890_ , 2017.