A massive variable flavour number scheme for the Drell-Yan process

R. Gauld1

1 Nikhef, Science Park 105, NL-1098 XG Amsterdam, The Netherlands

*<EMAIL_ADDRESS>
27th August 2024

## Abstract

The prediction of differential cross-sections in hadron-hadron scattering
processes is typically performed in a scheme where the heavy-flavour quarks
($c,b,t$) are treated either as massless or massive partons. In this work, a
method to describe the production of colour-singlet processes which combines
these two approaches is presented. The core idea is that the contribution from
power corrections involving the heavy-quark mass can be numerically isolated
from the rest of the massive computation. These power corrections can then be
combined with a massless computation (where they are absent), enabling the
construction of differential cross-section predictions in a massive variable
flavour number scheme. As an example, the procedure is applied to the low-mass
Drell-Yan process within the LHCb fiducial region, where predictions for the
rapidity and transverse-momentum distributions of the lepton pair are
provided. To validate the procedure, it is shown how the $n_{f}$-dependent
coefficient of a massless computation can be recovered from the massless limit
of the massive one. This feature is also used to differentially extract the
massless $\text{N}^{3}\text{LO}$ coefficient of the Drell-Yan process in the
gluon-fusion channel.

###### Contents

  1. 1 Introduction
  2. 2 (De)constructing the massive calculation
  3. 3 Heavy-quark mass slicing to $\mathcal{O}(\alpha_{s}^{3})$
  4. 4 Intrinsic charm contributions
  5. 5 Constructing the M-VFNS
  6. 6 Computational set-up.
  7. 7 Numerical validation and $\mathcal{O}(\alpha_{s}^{3})$ results
  8. 8 Differential distributions
  9. 9 Conclusions

##  1 Introduction

The prediction of high-energy scattering processes which involve initial-state
hadrons is crucial for understanding the physics of hadron collisions in both
controlled environments (such as the LHC) as well as a range of naturally
occurring scattering processes (such as those involving cosmic rays). In
general, the starting point for the theoretical description that describes
high-energy interactions in these collisions is a factorisation theorem [1] of
the form

$\displaystyle\frac{{\rm d}\sigma}{{\rm d}Q^{2}\,{\rm d}X}\sim$
$\displaystyle\sum_{a,b}\int{\rm d}\xi_{A}\,{\rm
d}\xi_{B}f_{a/A}(\xi_{A},\mu)f_{b/B}(\xi_{B},\mu)$ $\displaystyle\times{\rm
d}\hat{\sigma}_{ab}\left(\frac{x_{A}}{\xi_{A}},\frac{x_{B}}{\xi_{B}},Q;\alpha_{s}(\mu),\frac{\mu}{Q}\right)\,.$
(1)

This theorem separates the full scattering process into a partonic scattering
process involving the scattering of the hadron constituents $a,b$ (e.g. quarks
and gluons), and a set of parton distribution functions (PDFs) $f(x,Q)$ which
describe the probability distribution of the internal content of the hadron as
a function of hadron momentum-fraction and virtuality carried by the
constituent particle. The energy-scale $Q$ denotes a representative scale of
the scattering process, e.g. the dilepton invariant mass in the Drell-Yan (DY)
process [2], and $X$ is a hadronic level observable such as the rapidity of
the dilepton system.

A primary consideration when applying a factorisation theorem of this form is
the treatment of heavy-flavour quarks (e.g. charm and beauty). For example,
if/when it is a good approximation to consider these quarks as massless
partons, or whether to retain the exact mass dependence of the heavy quarks.
When treated as a massless parton, the heavy-flavour quark can be considered
as an active parton in the perturbative evolution of PDFs as well as the
strong-coupling $\alpha_{\mathrm{s}}$. This approach is often convenient as,
through this evolution, it allows to account for (to all orders) a class of
logarithmic corrections to the scattering process of the form
$\alpha_{\mathrm{s}}^{i}\ln[m/Q]^{j}$ for $i\geq j$, where $m$ is the heavy-
flavour quark mass. Instead, when considered as a massive parton, the impact
of the heavy-quark mass can be incorporated exactly up to the known
perturbative (fixed order) accuracy of the partonic cross-section ${\rm
d}\hat{\sigma}$. This allows for the computation of the same logarithmic
corrections as in the massless case outlined above (limited to fixed-order
accuracy only), and in addition power corrections of the form $m/Q$ which are
absent in the massless calculation.

Alternatively one can develop a scheme which combines these approaches,
providing a uniform description of the scattering process across arbitrary
energy scales—such a description is provided by a massive variable flavour
number scheme. This topic has been studied in various contexts in the past [3,
4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21], and with
particular focus on the process of lepton-nucleon scattering [22, 23, 24, 25,
26, 27, 28, 29, 30, 31]. In the latter case, it is well understood how to
apply such a formalism to the description of nucleon structure functions. Due
to the more rich structure of hadron-hadron scattering processes, the
development and application of such a formalism is relatively less mature. It
has been discussed for identified-hadron production [6, 32, 8, 20], processes
with flavoured-jets [33, 34, 35], inclusive quantities [36, 15, 37, 38, 39],
exclusive quantities in the framework of SCET [40], and also within the
context of Parton Showers [34, 41, 42].

The goal of this work is to revisit this topic for the production of colour-
singlet processes, focussing on the neutral-current DY process. LHC
measurements of this process have now reached per-mille level accuracy [43],
massless $\text{N}^{3}$LO predictions of such processes have now been obtained
[44, 45, 46], and precise computations of the transverse-momentum spectrum are
available at fixed-order [47, 48] and beyond [49, 50, 51, 52, 53, 54, 55].
Given this progress, I believe it is important to unambiguously assess the
importance of heavy-quark mass effects for fully differential collider physics
predictions. To do so, I develop a method that allows to numerically extract
(at the differential level) the contribution of massive power-corrections to
the hadronic cross-section. The method is applicable to arbitrary processes,
provided the considered observables are inclusive in QCD radiation and/or
infrared and collinear safe. In anticipation of a measurement, and as an
application and validation of the procedure, the low-mass Drell-Yan process is
considered within the fiducial volume of the LHCb experiment. As a by-product
of this work, I also show how the presented method can be used to obtain
differential information on the DY cross-section, which is used to extract the
$\text{N}^{3}$LO contribution to this process in the $gg$-channel.

##  2 (De)constructing the massive calculation

The general structure of the prediction of a differential hadronic-level
cross-section involving a single massive quark can be written as

${\rm d}\sigma^{\rm M}={\rm d}\sigma^{\rm m=0,n_{f}}+{\rm
d}\sigma^{\rm\ln[m]}+{\rm d}\sigma^{\rm pc}\,.$ (2)

The three contributions on RHS of Eq. (2) are: the $n_{f}$-dependent part of
the calculation which is present when $m=0$; those terms which depend
logarithmically on $m$ which diverge in the limit $m\to 0$; and all remaining
contributions that take the form of power corrections (labelled ‘pc’), and
which vanish in the limit $m\to 0$. The first two terms on the RHS of Eq. (2)
are also present in a massless calculation (as they define the $m\to 0$ limit
of the massive calculation), while the power corrections are uniquely
described by the massive calculation. The core idea of this work is that the
contribution ${\rm d}\sigma^{\rm pc}$ can be numerically isolated by directly
calculating all other terms appearing in Eq. (2). At fixed-order accuracy,
this isolation procedure should be applicable to arbitrarily differential
observables, provided they are inclusive with respect to QCD radiation and/or
are infrared and collinear safe such that the zero-mass limit is well defined.
This includes the differential description of a colour-singlet system (which
will be the focus of this work), but also applies to processes involving
hadronic jets (including those with identified flavour [56]). The application
to identified hadron production is slightly different (due to the presence of
final-state mass singularities), and has been discussed in the past [6, 32, 8,
20].

To illustrate how the procedure is performed, the neutral-current DY process
(i.e. ${\mathrm{p}}{\mathrm{p}}\to{\mathrm{\ell\bar{\ell}}}+X$) will be
considered, and a description of how to evaluate each of the terms appearing
in Eq. (2) is given. The current availability and the perturbative accuracy of
these terms is also described.

Massive computation, ${\rm d}\sigma^{\rm M}$. The cross-section ${\rm
d}\sigma^{\rm M}$ appearing on the LHS of Eq. (2) denotes the contribution
from a single massive quark with mass $m$ to the hadronic scattering process.
For the DY process, the presence of a massive quark alters the calculation
starting at $\mathcal{O}(\alpha_{s}^{2})$. The mass enters the calculation
explicitly in subprocesses of the form
$ab\to{\mathrm{\ell\bar{\ell}}}+Q\bar{Q}$ (where $a,b$ denote massless partons
and $Q$ the massive quark), but also implicitly enters the lower multiplicity
subprocesses $ab\to{\mathrm{\ell\bar{\ell}}}(+c)$ either in double-virtual
corrections or through the definition of UV renormalisation counter-terms—see
the Appendix of [57] for a detailed discussion. It should be clear that it is
necessary to consider all contributions of the massive quark (whether they
appear explicitly or not).

A massive calculation of the DY process is not available at
$\mathcal{O}(\alpha_{s}^{3})$. This requires perturbative ingredients, such as
various two and three-loop corrections involving a closed massive fermion
loop, which are currently unknown.

Zero mass computation, ${\rm d}\sigma^{\rm m=0,n_{f}}$. When considered
massless, the quark $Q$ still contributes to the same subprocesses as in the
massive computation described above (but with zero mass). This contribution
can be computed directly after extracting the $n_{f}$-dependent part of the
massless partonic cross-section at this order. Due to the presence of single-
and double-unresolved emissions in the differential calculation, this
extraction must also be applied to the (integrated) subtraction/slicing terms
which are required to regulate these emissions.

While first differential results for the massless DY cross-section at
$\mathcal{O}(\alpha_{s}^{3})$ have been presented [46] (relying on the NNLO
QCD calculation for $Zj$ [47] reported in [52]), a careful (and lengthy)
computation is required to extract the $\mathcal{O}(\alpha_{s}^{3}n_{f})$
component. A differential calculation of the massless
$\mathcal{O}(\alpha_{s}^{3}n_{f})$ contribution to DY is therefore currently
unavailable.

Logarithmic computation, ${\rm d}\sigma^{\rm\ln[m]}$. Provided that QCD
inclusive and/or infrared- and collinear-safe observables are considered, the
logarithmic dependence of the massive cross-section on the heavy-quark mass
$m$ is of collinear origin. This behaviour is universal, and it can be
described with knowledge of a set of decoupling relations which describe how
parameters (e.g. $\alpha_{s}$ and PDFs) in a theory with a massive quark are
mapped (at fixed-order accuracy) into an effective theory where that quark is
treated as massless. With this information, it becomes possible to construct
the logarithmic behaviour of the differential cross-section using only
massless inputs.

This construction requires the decoupling relation for $\alpha_{\mathrm{s}}$
(and $m$) which is known analytically to high perturbative-order [58], and
also available with public software (see for example [59]). The corresponding
relations for the PDFs are provided in the form of massive Operator Matrix
Elements (OMEs, and denoted $\hat{A}_{ab}$) which describe the transition
between the partonic states $a\to b$. The perturbative structure of these
objects has been studied at great length in the past, and calculations of the
massive OMEs are available at two-loop [60, 61, 62, 63, 24, 10, 64, 65], and
three-loop [66, 67, 68, 69, 70] order. Notably, these OMEs also define the
matching conditions/decoupling relations which allow to construct a VFNS for
PDFs—see for example Eq. (12-15) of [70] (and originally [24]).

To construct the logarithmic cross-section for the DY process, one has to
consider convolutions of the form

$\displaystyle\hat{A}^{(i)}_{ab}\otimes\hat{A}^{(j)}_{cd}\otimes{\rm
d}\hat{\sigma}^{(k),m=0}_{bd\to{\mathrm{\ell\bar{\ell}}}+X}\,,$ (3)

where the superscripts $(i-k)$ denote the perturbative order of the OMEs and
the massless partonic scattering cross-section (${\rm
d}\hat{\sigma}^{(k),m=0}_{bd\to{\mathrm{\ell\bar{\ell}}}+X})$. All of the
$\hat{A}^{(i)}_{ab}$ inputs required to construct ${\rm d}\sigma^{\ln[m]}$ up
to $\mathcal{O}(\alpha_{s}^{2})$ (i.e. $i+j\leq 2$) have been presented in
[24]. The results presented in [66, 67, 68, 69, 70] should also allow to
extend this calculation to $\mathcal{O}(\alpha_{s}^{3})$. For $(k)\geq 1$, the
decoupling relation for the strong coupling
$\Delta^{(i)}_{n_{f}}(\alpha_{\mathrm{s}})$ is also required, which can be
applied as a multiplicative factor to Eq. (3). The perturbative expansion for
$\Delta^{(i)}_{n_{f}}(\alpha_{\mathrm{s}})$ is reported in Eq. (20) of [59].
Working in the $\overline{\mathrm{MS}}$ scheme for the strong coupling, and
defining the heavy quark mass in the on-shell scheme (which is consistent with
that of the OME calculation in [24]), the expansion is

$\displaystyle\Delta_{n_{f}}(\alpha_{\mathrm{s}})=1+\frac{\alpha_{\mathrm{s}}}{2\pi}\left(-\frac{1}{3}\ln\left[\frac{\mu^{2}}{m^{2}}\right]\right)+\left(\frac{\alpha_{\mathrm{s}}}{2\pi}\right)^{2}\left(-\frac{7}{6}-\frac{19}{6}\ln\left[\frac{\mu^{2}}{m^{2}}\right]+\frac{1}{9}\ln^{2}\left[\frac{\mu^{2}}{m^{2}}\right]\right)+\mathcal{O}(\alpha_{s}^{3})\,.$
(4)

Here, $\alpha_{\mathrm{s}}$ denotes the strong coupling defined in the
massless scheme (where the heavy quark is included in the running) at the
scale $\mu$. Up to $\mathcal{O}(\alpha_{s}^{2})$, the relevant expansion which
is required to build the (partonic) logarithmic cross-section is

$\displaystyle{\rm d}\hat{\sigma}^{\rm\ln[m],(2)}_{q\bar{q}}$
$\displaystyle=\left(\hat{A}^{(2)}_{qq,Q}+\hat{A}^{(2)}_{\bar{q}\bar{q},Q}\right)\otimes{\rm
d}\hat{\sigma}^{(0)}_{q\bar{q}}+\Delta_{n_{f}}^{(1)}(\alpha_{\mathrm{s}})\cdot{\rm
d}\hat{\sigma}^{(1)}_{q\bar{q}}\,,$ $\displaystyle{\rm
d}\hat{\sigma}^{\rm\ln[m],(2)}_{gg}$
$\displaystyle=\hat{A}_{gQ}^{(1)}\otimes\hat{A}_{g\bar{Q}}^{(1)}\otimes{\rm
d}\hat{\sigma}^{(0)}_{Q\bar{Q}}+\hat{A}_{gQ}^{(1)}\otimes\left({\rm
d}\hat{\sigma}^{(1)}_{Qg}+{\rm
d}\hat{\sigma}^{(1)}_{gQ}\right)+\left(Q\leftrightarrow\bar{Q}\right)\,.$ (5)

Constructed in this way (i.e. using massless inputs) the logarithmic
calculation will also contain those terms which are independent of $m$. They
are generated by the constant terms contained in $\hat{A}_{ab}$ and
$\Delta_{n_{f}}(\alpha_{\mathrm{s}})$—i.e. those which define the de-coupling
across heavy-flavour thresholds in a variable flavour number scheme. It is
necessary to account for these terms as they are part of the massive
calculation (i.e. they appear on the LHS of Eq. (2)), but are not generated
when ${\rm d}\sigma^{\rm m=0,n_{f}}$ is computed with inputs (PDFs and
$\alpha_{\mathrm{s}}$) defined in the massive scheme (e.g. $n_{f}^{\rm max}=4$
for the $b$-quark).

##  3 Heavy-quark mass slicing to $\mathcal{O}(\alpha_{s}^{3})$

Following from the discussion in the previous Section, it is clear that all
ingredients required to extract the power corrections for the DY process are
known up to $\mathcal{O}(\alpha_{s}^{2})$. This extraction is achieved
numerically by evaluating the first three terms appearing in Eq. (2) and
solving for ${\rm d}\sigma^{\rm pc}$. To validate this procedure, is it
important to test that the extracted power corrections vanish in the limit
$m\to 0$, which can be done by performing the extraction for decreasing values
of $m$. Extracted in this way, the power corrections will only vanish provided
that ${\rm d}\sigma^{\ln[m]}$ reproduces the logarithmic behaviour of the
massive cross-section in the limit $m\to 0$, and that the calculation of the
constant term ${\rm d}\sigma^{\rm m=0,n_{f}}$ is correct.

Viewed in another way, one can also use the small-mass limit to numerically
extract ${\rm d}\sigma^{\rm m=0,n_{f}}$ when it is unknown. This can be
achieved if both massive and logarithmic calculations are known at the desired
perturbative order, by performing a fit to the constant difference $\left({\rm
d}\sigma^{M}-{\rm d}\sigma^{\ln[m]}\right)$ in the limit $m\to 0$. Once this
constant is known, the power corrections can also be extracted at the physical
value of the heavy-quark mass. In practice, this corresponds to a global
slicing method, where the heavy-quark mass parameter $m$ acts as a collinear
regulator.

This technique is noteworthy, as it can be used to extract differential
results for the DY cross-section at $\mathcal{O}(\alpha_{s}^{3})$ in the
gluon-fusion channel. This is possible because both the massive and
logarithmic calculations are available at this order. The massive calculation
is simply the NLO QCD correction to the subprocess
$gg\to\ell\bar{\ell}Q\bar{Q}$ which can be obtained with automated codes such
as aMC@NLO [71, 72], and the logarithmic cross-section can be constructed
using the two-loop OMEs given in [24]. The (partonic) logarithmic cross-
section at $\mathcal{O}(\alpha_{s}^{3})$ is constructed following the same
procedure which resulted in Eq. (2), but expanded to one order higher. At this
order new technical features are encountered, which are briefly discussed
below.

At $\mathcal{O}(\alpha_{s}^{3})$, it is necessary to consider convolutions of
the form

$\displaystyle\hat{A}^{(1)}_{gQ}\otimes{\rm
d}\hat{\sigma}^{(2),m=0}_{Qg\to{\mathrm{\ell\bar{\ell}}}+X}\,,$ (6)

where ${\rm d}\hat{\sigma}^{(2),m=0}_{Qg\to{\mathrm{\ell\bar{\ell}}}+X}$ is
the second-order massless partonic cross-section and $\hat{A}^{(1)}_{gQ}$ the
one-loop OME in the $Qg$ channel. The partonic cross-section appearing in this
convolution contains contributions from double-real, real-virtual, and double-
virtual phase-space configurations—with (integrated) subtraction terms
appearing at each level. It is therefore necessary to convolute
$\hat{A}^{(1)}_{gQ}$ with all terms at all levels, which at the real-virtual
and virtual-virtual level includes iterated convolutions with integrated
subtraction terms. To better clarify this issue, it is sufficient to consider
the construction of the logarithmic cross-section at finite
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ values where the partonic cross-
section $\hat{\sigma}^{(2),m=0}_{Qg\to{\mathrm{\ell\bar{\ell}}}+X}$
effectively becomes NLO accurate (there can be maximally a single unresolved
QCD emission at this order). Practically, in this case, one has to consider
terms of the form

$\displaystyle+\hat{A}^{(1)}_{gQ}$ $\displaystyle\otimes\left({\rm
d}\hat{\sigma}^{R,m=0}_{Qg\to{\mathrm{\ell\bar{\ell}}}+X}-{\rm
d}\hat{\sigma}^{Rs,m=0}_{Qg\to{\mathrm{\ell\bar{\ell}}}+X}\right)$
$\displaystyle+\hat{A}^{(1)}_{gQ}$ $\displaystyle\otimes\left({\rm
d}\hat{\sigma}^{V,m=0}_{Qg\to{\mathrm{\ell\bar{\ell}}}+X}-{\rm
d}\hat{\sigma}^{Vs,m=0}_{Qg\to{\mathrm{\ell\bar{\ell}}}+X}\right)\,,$ (7)

where the superscripts $R,Rs,V,Vs$ denote real, real-subtraction, virtual, and
virtual-subtraction terms which are present in a differential NLO calculation.
The convolution with $\hat{A}^{(1)}_{gQ}$ must be applied to all terms. As the
virtual-subtraction cross-section contains mass-factorisation terms and
integrated subtraction terms which are distributions in a collinear variable,
the convolution of these terms with $\hat{A}^{(1)}_{gQ}$ leads to the
aforementioned iterated convolutions. To obtain a result which is valid at all
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ values, the above procedure is
applied to the NNLO accurate partonic cross-section.

At $\mathcal{O}(\alpha_{s}^{3})$, it is also important to consider the scheme
dependence of the various inputs which enter the construction in Eq. (3). For
example, the OME calculation of [24] is valid when the input PDFs are
renormalised in the $\overline{\mathrm{MS}}$ scheme with $n_{f}-1$ light
flavours while the strong coupling is renormalised with $n_{f}$ light
flavours. Also the first-order partonic cross-section ${\rm
d}\hat{\sigma}^{(1)}_{gQ}$ appearing in line 2 of Eq. (2) is defined with
$n_{f}-1$ light flavours (i.e. for $\alpha_{\mathrm{s}}$ and the gluon PDF
which is not convoluted with the OME). Practically, a scheme conversion is
applied such that all inputs are defined with inputs (PDFs and
$\alpha_{\mathrm{s}}$) that are defined in the $\overline{\mathrm{MS}}$ scheme
with $n_{f}$ light flavours. The relevant scheme conversions are

$\displaystyle\alpha_{\mathrm{s}}^{[n_{f}-1]}$
$\displaystyle=\alpha_{\mathrm{s}}^{[n_{f}]}\left(1+\Delta^{(1)}(\alpha_{\mathrm{s}}^{[n_{f}]})\right)+\mathcal{O}(\alpha_{s}^{2})\,,$
$\displaystyle g^{[n_{f}-1]}(x,\mu_{F}^{2})$
$\displaystyle=g^{[n_{f}]}(x,\mu_{F}^{2})-\frac{\alpha_{\mathrm{s}}^{[n_{f}]}}{2\pi}\hat{A}_{gg,Q}^{S,(1)}\otimes
g^{[n_{f}]}(x,\mu_{F}^{2})+\mathcal{O}(\alpha_{s}^{2})\,,$ (8)

with the definition

$\displaystyle\hat{A}_{gg,Q}^{S,(1)}\left(z,\mu_{F}^{2}/m^{2}\right)=-\delta(1-z)\frac{1}{3}\ln\left[\frac{\mu_{F}^{2}}{m^{2}}\right]\,.$
(9)

Notice that when this conversion is applied to equal powers of
$\alpha_{\mathrm{s}}$ and the gluon PDF, the dependence on $m$ vanishes and a
logarithm of the ratio $\mu_{F}/\mu_{R}$ remains.

##  4 Intrinsic charm contributions

So far, the discussion has implicitly assumed that the contribution from
massive initial-state quarks is absent. This is consistent with the assumption
that there are no intrinsic heavy-flavour PDFs, which is the set-up of most
modern global PDF fitting groups. In contrast, the NNPDF collaboration have
relaxed this assumption [73] (see also [74]), and now fits for an intrinsic
charm quark PDF as part of the nominal fit. In this case, the formalism
outlined above can also be applied to extract the massive power-corrections
associated to initial-state charm quarks. This is done in this work, extending
the previous results for DIS [30, 75] and inclusive observables [21], to the
fully differential level. This requires the use of the OMEs for massive
initial states originally computed in [4] which have also been presented in
the Appendix of [21].

As an aside, I note that the general factorisation theorem for the computation
of hadronic-level cross-section predictions involving massive-initial state
quarks is known to be violated at $\mathcal{O}(\alpha_{s}^{2})$. This topic
has been studied in the past [76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87],
and has received recent attention in [88] in the context of the Drell-Yan
process. A deeper theoretical understanding of factorisation theorems for
massive-initial states remains desirable today.

##  5 Constructing the M-VFNS

The massive variable flavour number scheme (M-VFNS) is constructed by
combining a massless calculation with that of the massive power-corrections
outlined above, for differential predictions, according to

${\rm d}\sigma^{\rm M-VFNS}={\rm d}\sigma^{\rm m=0}+\sum^{n_{f}^{\rm
max}}_{i=c,b,...}{\rm d}\sigma_{i}^{\rm pc}\,.$ (10)

In this matching formula, the first term ${\rm d}\sigma^{\rm m=0}$ is the
massless computation (i.e. that in a zero mass variable flavour number scheme
with $n_{f}^{\rm max}$ flavours, ZM-VFNS), and the second term denotes the
power corrections which are obtained by re-arranging Eq. (2). The power
corrections can be evaluated separately for each of the heavy-flavour quarks
(at higher-orders, one could also extend the formalism to deal with the
presence of two-mass contributions). The master formula Eq. (10) is similar to
those which have been presented for DIS Structure Functions (e.g. [29]), where
${\rm d}\sigma^{\rm pc}$ is often written as the difference ${\rm
d}\sigma^{\rm pc}=\left({\rm d}\sigma^{M}-{\rm d}\sigma^{M\to 0}\right)$.

With respect to either a massive or a massless approach, the benefits of this
construction are that a resummation of a class of collinear logarithms
involving the heavy-quark mass $m$ (through PDF and $\alpha_{s}$
Renormalisation Group evolution) are included to all orders, and the exact
heavy-quark mass dependence is included to the fixed-order accuracy to which
${\rm d}\sigma_{i}^{\rm pc}$ is known.

Details of the computational set-up used for this work are provided in the
following Section, before providing a numerical validation of the procedure
and phenomenological results relevant for a measurement by the LHCb
collaboration.

##  6 Computational set-up.

Theoretical implementation. The predictions of the various differential cross-
sections which enter the construction of the M-VFNS for the DY process are
provided with a specialised Monte Carlo programme. It was originally purposed
to enable the construction of a M-VFNS for the
${\mathrm{p}}{\mathrm{p}}\to{\mathrm{Z}}+b\text{-jet}$ process [35]. The
programme has since been extended to contain all ingredients which are
required for the computation of the process
${\mathrm{p}}{\mathrm{p}}\to\ell\bar{\ell}+X$ up to
$\mathcal{O}(\alpha_{s}^{2}$), which may involve involve massless or massive
heavy-flavour QCD partons. This includes those axial contributions arising due
the presence of heavy-quark triangle diagrams, see for example [89, 90, 91].
Processes involving massive initial states are instead limited to
$\mathcal{O}(\alpha_{s}$).

These computations are performed using a combination of Dipole subtraction
[92] to treat the presence of single unresolved emissions (see [93, 18] for
massive initial-states), and N-jettiness slicing for double unresolved
emissions [94]. This implementation relies on many existing results, which
include: amplitudes [95, 90, 91, 96, 97, 57, 98]; N-jettiness inputs [99, 100,
101, 102, 103, 104, 105, 106]; several OpenLoops libraries for tree-level
amplitudes [107]; as well as a number of results manually computed with the
aid of FeynArts [108] and FormCalc [109]. Beyond fixed-order, the programme
also facilitates the computation of resummed predictions of the
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ spectrum at NNLL accuracy using a
combination of results from [110, 111, 112]. The numerical integration of all
contributions in the Monte Carlo programme are performed with the VEGAS
algorithm as implemented in CUBA [113].

Numerical inputs. All predictions are provided with the NNPDF3.1 NNLO PDF set
[73] with $\alpha_{s}(M_{\mathrm{Z}})=0.118$ (with $n_{f}^{\rm max}=5$), where
the PDF and $\alpha_{s}$ values are accessed via LHAPDF [114]. These PDFs are
used as an input to all calculations, which requires the application of a
renormalisation scheme change to some of the inputs which enter the
$\mathcal{O}(\alpha_{s}^{3})$ calculation in the $gg$-channel. The use of
these PDFs grids (and the corresponding $\alpha_{\mathrm{s}}$ values)
practically defines how the resummation is implemented within the first term
appearing in Eq. (10). The values of the on-shell heavy quark masses in this
PDF set are $m_{c,b,t}^{\rm pdf}=1.51,4.92,172.5\leavevmode\nobreak\
\mathrm{GeV}$. It should be noted that the boundary condition for the PDF set
is defined at $Q_{0}=1.65\leavevmode\nobreak\ \mathrm{GeV}$, which is larger
than $m_{c}^{\rm pdf}$. This information is relevant as it is used to derive
the the static charm-quark PDF $f_{c}(x)$, according to the de-coupling
relations calculated in [4].

All calculations are performed in the Complex Mass Scheme [115], with
Electroweak inputs defined in the $G_{\mu}$-scheme following [107]. The
following values for the numerical inputs are used
$M_{{\mathrm{Z}}}^{\mathrm{os}}=91.1876\leavevmode\nobreak\ \mathrm{GeV}$,
$\Gamma_{{\mathrm{Z}}}^{\mathrm{os}}=2.4952\leavevmode\nobreak\ \mathrm{GeV}$,
$M_{{\mathrm{W}}}^{\mathrm{os}}=80.379\leavevmode\nobreak\ \mathrm{GeV}$,
$\Gamma_{{\mathrm{W}}}^{\mathrm{os}}=2.085\leavevmode\nobreak\ \mathrm{GeV},$
and $G_{\mu}=1.16638\times 10^{-5}\leavevmode\nobreak\ \mathrm{GeV}^{-2}$. The
massless DY computations at $\mathcal{O}(\alpha_{s}^{2})$ use the N-jettiness
slicing method with a technical parameter of $\tau_{\rm
cut}=10^{-3}\leavevmode\nobreak\ \mathrm{GeV}$. Such a small value was chosen
(at substantial CPU cost) to suppress the impact of missing power corrections
beyond those included via [106] (see also [116] for a recent discussion).

For the results shown in the following Sections, an uncertainty due to the
impact of missing higher-order corrections is assessed by varying the values
of $\mu_{R}$ and $\mu_{F}$ by a factor of two around the dynamical scale
$\mu_{0}\equiv E_{\mathrm{T},{{\mathrm{\ell\bar{\ell}}}}}$ (the transverse
mass of the dilepton pair), with the constraint that
$\frac{1}{2}\leq\mu_{F}/\mu_{R}\leq 2$. When the M-VFNS is constructed
according to Eq. (10), the scale uncertainties are correlated between the
power corrections and the massless computations. Where shown, PDF
uncertainties have been obtained from individual replica predictions ($i$)
calculated in the following way:

$\displaystyle{\rm d}\sigma_{i}=K\,{\rm
d}\sigma_{i}[\mathcal{O}(\alpha_{s})]\,,\qquad K=\frac{{\rm
d}\sigma_{0}[\mathcal{O}(\alpha_{s}^{2})]}{{\rm
d}\sigma_{0}[\mathcal{O}(\alpha_{s})]}\,.$ (11)

That is to say that a differential K-factor is calculated for the central PDF
member at $\mathcal{O}(\alpha_{s}^{2})$, and then applied to each of the
individual replica cross-sections which are computed at
$\mathcal{O}(\alpha_{s})$.

LHCb fiducial definition. In anticipation of a measurement of the process
${\mathrm{p}}{\mathrm{p}}\to\ell\bar{\ell}+X$ by the LHCb collaboration at
$\sqrt{S}=13\leavevmode\nobreak\ \mathrm{TeV}$, the procedure outlined in the
previous Section is both validated and applied in the fiducial region of the
LHCb experiment. The predictions will be performed double differentially with
respect to the invariant mass of the dilepton pair within the range
$m_{\ell\bar{\ell}}\in[4,200]\leavevmode\nobreak\ \mathrm{GeV}$ and either the
transverse momentum ($p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$) or rapidity
($y_{{\mathrm{\ell\bar{\ell}}}}$) of the dilepton pair. The following set of
cuts are applied to the charged leptons:
$p_{\mathrm{T},\ell}>1.5\leavevmode\nobreak\
\mathrm{GeV},|p_{\ell}|>20\leavevmode\nobreak\
\mathrm{GeV},2.0<\eta_{\ell}<4.5$.

Predictions have been generated at fixed-order accuracy for both
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ and $y_{{\mathrm{\ell\bar{\ell}}}}$
distributions in 20 invariant mass bins (a total of 400 bins). In this work, I
have chosen to present the results within the invariant mass region of
$m_{\ell\bar{\ell}}\in[12.5,13.5]\leavevmode\nobreak\ \mathrm{GeV}$. This
region is of phenomenological interest as it provides sensitivity to the input
PDFs at small-$x$, without being overwhelmed by perturbative uncertainties
which grow in the very low $m_{{\mathrm{\ell\bar{\ell}}}}$ regime. At the same
time, this invariant mass region is sufficiently small that the fixed-order
predictions within the range of
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}\in[2.5,11.0]\leavevmode\nobreak\
\mathrm{GeV}$ (which will be the focus of the
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ measurement) are expected to be
reliable. At higher values of $m_{\ell\bar{\ell}}$, a resummation of Sudakov
logarithms of the form
$\ln[p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}/m_{\ell\bar{\ell}}]$ is
necessary.

The numerical validation of the theoretical procedure (to follow) will be
performed in the following kinematic regimes:

$\displaystyle{\rm Fiducial}$ $\displaystyle:\quad
p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}\leavevmode\nobreak\ {\rm
inclusive}\,,$ $\displaystyle p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm
low}$ $\displaystyle:\quad p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}\leq
2.5\leavevmode\nobreak\ \mathrm{GeV}\,,$ $\displaystyle
p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm high}$ $\displaystyle:\quad
p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}\geq 2.5\leavevmode\nobreak\
\mathrm{GeV}\,.$ (12)

The inclusion of this additional restriction in
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ is relevant as this defines the bin
edge of the on-going $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ measurement.
Clearly, the Fiducial volume is the sum of the latter contributions.

##  7 Numerical validation and $\mathcal{O}(\alpha_{s}^{3})$ results

Validation at $\mathcal{O}(\alpha_{s}^{2})$. To first validate the procedure,
the small-mass limit of the massive cross-section ${\rm d}\sigma^{\rm M}$
appearing in Eq. (2) is considered within the LHCb Fiducial region. This
contribution is computed at $\mathcal{O}(\alpha_{s}^{2})$ and is compared to
the logarithmic cross-section ${\rm d}\sigma^{\rm\ln[m]}$ at the same order.
In both cases, the scales are set to
$\mu_{\mathrm{F}}=\mu_{\mathrm{R}}=E_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$
and the contribution from $b$-quarks are considered. The results for
$c$-quarks are qualitatively similar, and differ in magnitude due to the
coupling of the quarks with the exchanged gauge-boson.

The results are shown in Fig. 1 with a breakdown into $q\bar{q}$\- and
$gg$-initiated channels ($q$ indicating a light-flavour quark). In the lower-
panel, the cross-section difference $\left({\rm d}\sigma^{\rm M}-{\rm
d}\sigma^{\rm\ln[m]}\right)$ is shown alongside the direct calculation of
${\rm d}\sigma^{\rm m=0,n_{f}}$. The direct and indirect (obtained via the
$m\to 0$ limit) methods of calculating ${\rm d}\sigma^{\rm m=0,n_{f}}$
coincide, confirming the structure of the massive calculation presented in Eq.
(2).

It is important to highlight that the massive calculation necessarily includes
a sum over all contributions of the massive quark. To understand why this is
critical, one can consider the double-real subprocess
$q\bar{q}\to{\mathrm{\ell\bar{\ell}}}+Q\bar{Q}$. When the heavy-quark pair is
emitted in a soft and double-collinear configuration, a triple-logarithmic
contribution of the form $\alpha_{s}^{2}\ln[m]^{3}$ is generated. This triple
logarithm is cancelled (at the level of the differential cross-section) by the
exchange of a virtual soft and double-collinear massive quark-pair which is
present in the two-loop form factor for a massless quark pair. Which is to
say, without including all contributions involving the massive quark, the
cross-section prediction will contain logarithmic sensitivity to the heavy-
quark mass which is not described by universal structure of the OMEs and
$\alpha_{\mathrm{s}}$ decoupling relations. This is effectively a statement of
the KLN theorem [117, 118].

Figure 1: Upper panel: absolute cross-section for the full massive and
logarithmic calculations at $\mathcal{O}(\alpha_{s}^{2})$. Lower panel: the
cross-section difference (massive-logarithmic) as compared to the direct
computation of the massless $n_{f}$-dependent cross-section. A breakdown into
partonic channels is provided.

Extension to $\mathcal{O}(\alpha_{s}^{3})$. As discussed, the perturbative
ingredients required to extend the procedure to $\mathcal{O}(\alpha_{s}^{3})$
are available only for the $gg$-channel. In this case, no direct calculation
of ${\rm d}\sigma^{\rm m=0,n_{f}}$ is available, but it can be extracted from
a numerical fit to the difference $\left({\rm d}\sigma^{\rm M}-{\rm
d}\sigma^{\rm\ln[m]}\right)$ in the $m\to 0$ limit.

This is done by generating data for the quantity $\left({\rm d}\sigma^{\rm
M}-{\rm d}\sigma^{\rm\ln[m]}\right)$ for several values of $m$ in the range of
$m\in[0.5,12]\leavevmode\nobreak\ \mathrm{GeV}$, and subsequently performing a
numerical fit. By inspecting Eq. (2), the resultant distribution should be
equal to the sum of the contributions from power corrections and the zero mass
computation. An ansatz of the following form is therefore used for the
numerical fit

$\displaystyle
f(m)=a_{0,0}+\sum_{i=1,j=0}a_{i,j}\,\left(m^{2}\right)^{i}\,\ln^{j}[m]\,.$
(13)

The form of this ansatz is motivated by the behaviour of the squared matrix-
element and phase space which both contain corrections of the form
$m^{2}/Q^{2}$. The integer $j$ is limited to $2(3)$ when the
$\alpha_{s}^{2(3)}$ coefficient is fitted, and a maximum value of $i=2$ is
considered in each case. The choice for $j$ is guided by the powers of
collinear logarithms which may be present at each order, whereas increasing
$i$ beyond 2 had little impact on the fit. The $m$-independent constant
$a_{0,0}$ is equivalent to ${\rm d}\sigma^{\rm m=0,n_{f}}$, while the
remaining terms describe the power corrections. Fitted in this way, all
$m$-independent information (such as dependence on $\mu$, which is chosen as
the dynamic scale $E_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$) is absorbed into
the $a_{i,j}$ coefficients.

Figure 2: The same as the lower panel of Fig. 1, focussing on the $gg$-channel
in two $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ regions. The result of a
numerical fit to the difference between the massive and logarithmic cross-
sections is shown, and compared to the direct calculation of the zero mass
constant at $\mathcal{O}(\alpha_{s}^{2})$.

The results are shown in Fig. 2, where a total of four fitted curves are
displayed corresponding to the two kinematic regimes of
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm low}$ and
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm high}$ at
$\mathcal{O}(\alpha_{s}^{2})$ and $\mathcal{O}(\alpha_{s}^{3})$. In addition
to the fitted central value, the fitted value of $a_{0,0}$ and its
corresponding uncertainty (indicated by a solid filled band) is displayed for
each of the curves. As noted, $a_{0,0}$ should correspond to ${\rm
d}\sigma^{\rm m=0,n_{f}}$, and it is therefore compared to the direct
computation of this quantity available at $\mathcal{O}(\alpha_{s}^{2})$. The
fit leads to a result which is consistent with the direct calculation,
providing confidence that the fitting procedure leads to reliable results.

From the fitted values of $a_{0,0}$ (which were produced for $b$ quarks), it
is possible to construct the full $n_{f}$-dependent massless cross-section in
the $gg$-channel. This is done by multiplying these results by a factor of
$F=n_{d}+n_{u}(Q_{u}/Q_{d})^{2}$, where $n_{u}$ and $n_{d}$ are the number of
down- and up-type quarks. This relation holds (at the per-mille level) for the
$gg$-channel as there is a direct coupling of the heavy-quark line to the
gauge-boson (which is dominated by photon exchange for
$m_{\ell\bar{\ell}}\in[12.5,13.5]\leavevmode\nobreak\ \mathrm{GeV}$). The
fitted (and, where available, direct computation) are summarised in Table 1.

| Order | Fiducial $[\mathrm{pb}]$ | $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm low}[\mathrm{pb}]$ | $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm high}[\mathrm{pb}]$
---|---|---|---|---
${\rm d}\sigma^{\rm m=0}_{gg}$ | $\alpha_{s}^{2}$ | +0.51(2) | +3.96(2) | -3.45(0)
${\rm d}\sigma^{\rm m=0}_{gg,{\rm fit}}$ | $\alpha_{s}^{2}$ | +0.49(4) | +3.97(4) | -3.49(4)
${\rm d}\sigma^{\rm m=0}_{gg,{\rm fit}}$ | $\alpha_{s}^{3}$ | +0.10(6) | +4.11(5) | -4.01(4)
Table 1: Predicted and fitted values of the coefficient of the zero-mass
computation up to $\mathcal{O}(\alpha_{s}^{3})$ in the $gg$-channel. The
results are for the central scale, and include an uncertainty due to the
fitting procedure and the statistical error.

Impact of massive power-corrections. To provide another validation of the
procedure, it is also useful to directly show the contribution from the
massive power-corrections ${\rm d}\sigma^{\rm pc}$ (including scale
variation).

A selection of such results are shown in Fig. 3, indicating that the power
corrections vanish in the limit $m\to 0$ for arbitrary scale choices. These
results include those from $q\bar{q}$ and $gg$ channels (where the heavy-
flavour quarks are produced in the final state) at
$\mathcal{O}(\alpha_{s}^{2})$, as well as $\mathcal{O}(\alpha_{s})$
contributions from massive initial-state charm quark contributions in the $cg$
channel. In the $gg$ and $cg$ channels, the results are displayed for the
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm low}$ and
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm high}$ regions to indicate
large cancellations which occur for the power corrections when integrated in
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$. It is worth noting that the power
corrections in the $gg$-channel have a different sign at the on-shell value of
the $c$\- and $b$-quark mass, which leads to an additional source of
cancellation. The results from massive $c\bar{c}$-initiated states were
negligibly small (due to the small value of the static charm PDF), and are not
shown here.

Figure 3: Massive power-corrections to the DY cross-section within the LHCb
fiducial region with $m_{\ell\bar{\ell}}\in[12.5,13.5]\leavevmode\nobreak\
\mathrm{GeV}$. The partonic-channels, perturbative orders, and considered
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ regions (unless inclusive) are
highlighted. The scale uncertainty of each of these predictions is shown.

So far, the results presented as a function of $m$ are used to validate the
general procedure introduced in this paper. The actual impact of the power-
corrections should be studied at the values of $m=m_{c,b}^{\rm pdf}$. As an
additional note, the construction of the logarithmic cross-section should be
performed in the same way that it is done for the massless calculation. For
example, the logarithmic contribution generated by the PDFs and $\alpha_{s}$
is only present when evaluated above heavy-flavour threshold (typically
$m^{\rm pdf}$). This approach ensures that in the limit of $m=m^{\rm
pdf}\to\infty$ (where the massive cross-section vanishes), the massive power
corrections obtained via Eq. (2) exactly cancel those contributions from the
zero mass computation.

To place the results of Fig. 3 in context, the magnitude of the power
corrections (evaluated at $m_{c,b}^{\rm pdf}$) should be compared to that of
the (total) massless computation. The results of this comparison are
summarised in Table 2. The first column indicates results within the Fiducial
region, while the second column shows results in the
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm high}$ region which includes
the constraint $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}\geq
2.5\leavevmode\nobreak\ \mathrm{GeV}$. The reference calculation ${\rm
d}\sigma^{\rm m=0}$ (which includes scale uncertainties) is NNLO QCD accurate
within the Fiducial region, and NLO QCD accurate at finite
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$. The $\alpha_{s}^{2}$ and
$\alpha_{s}^{3}$ coefficients of the massive power-corrections are shown in
the same kinematic regions, where a breakdown into those contributions from
charm and beauty quarks is given. Note that the charm-quark contributions in
the second row are $\mathcal{O}(\alpha_{s}^{2})$, as they include the
intrinsic contributions from both the born and $\alpha_{s}$ coefficients. The
$\alpha_{s}^{3}$ results are obtained from the functional fits shown in Fig.
2, evaluated at the on-shell mass values for each quark (with an appropriate
normalisation correction for the up-type quark).

Prediction | Order | Fiducial $[\mathrm{pb}]$ | $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{\rm high}[\mathrm{pb}]$
---|---|---|---
${\rm d}\sigma^{\rm m=0}$ | $\mathcal{O}(\alpha_{s}^{2})$ | $59.9^{+2.0}_{-5.6}$ | $46.2^{+5.1}_{-8.1}$
${\rm d}\sigma_{c}^{\rm pc}$ | $\mathcal{O}(\alpha_{s}^{2})$ | $+0.07$ | $-0.04$
${\rm d}\sigma_{b}^{\rm pc}$ | $\alpha_{s}^{2}$ | $+0.23$ | $+0.19$
${\rm d}\sigma_{c}^{\rm pc}(gg)$ | $\alpha_{s}^{3}$ | $+0.09(2)$ | $-0.04(2)$
${\rm d}\sigma_{b}^{\rm pc}(gg)$ | $\alpha_{s}^{3}$ | $+0.05(1)$ | $+0.10(0)$
Table 2: Predictions for the DY cross-section within the LHCb fiducial region
with $m_{\ell\bar{\ell}}\in[12.5,13.5]\leavevmode\nobreak\ \mathrm{GeV}$. The
contributions from the massless calculation and the massive power-corrections
are shown for the central scale. The scale uncertainties of the massless
$\mathcal{O}(\alpha_{s}^{2})$ prediction are indicated, and the uncertainties
in parenthesis correspond to the fit uncertainties of the $\alpha_{s}^{3}$
coefficients.

When compared to the massless prediction, the massive power-corrections
introduce a correction which is typically at the level of $\approx 0.5\%$.
These corrections are small in general, and there are a number of
cancellations which occur between different partonic channels, different quark
flavours (charm vs. beauty), and also across the
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ spectrum (see Fig. 2 and 3).
Overall, the sum of these corrections is negligible compared to the size the
perturbative uncertainty of the massless calculation. Similar behaviour was
found to persist for the entire $m_{\ell\bar{\ell}}$ range up to
$200\leavevmode\nobreak\ \mathrm{GeV}$.

##  8 Differential distributions

This work focusses on improving our understanding the role of massive quarks
in hadron-hadron scattering processes. However, in view of the on-going
measurement of low-mass DY at LHCb, I also take this opportunity to provide
some phenomenological results and recommendations.

As highlighted in Table 2, the overall normalisation of the massless cross-
section has a uncertainty due to scale variation as large as 9%. It is
therefore useful to consider normalised differential measurements where this
theoretical uncertainty is reduced, but sensitivity to the input PDFs is
retained. Such an example is

$\displaystyle\frac{1}{\sigma}\frac{{\rm d}\sigma}{{\rm
d}y_{{\mathrm{\ell\bar{\ell}}}}}\,,$ (14)

where $\sigma$ is the rapidity integrated cross-section for a given
$m_{\ell\bar{\ell}}$ region. This observable is also experimentally well
motivated as systematic uncertainties due to lepton reconstruction are
strongly correlated in rapidity (at fixed $m_{\ell\bar{\ell}}$).

Theoretical predictions for this quantity are shown in upper panel of Fig. 4
at $\mathcal{O}(\alpha_{s}^{2})$. In the lower panel, the various predictions
and uncertainties are shown normalised to that of the central M-VFNS
prediction constructed using Eq. (10). The power corrections are small, and
further cancel when constructing the normalised cross-section (as the
corrections are approximately flat in rapidity) resulting in a negligible
contribution. The PDF uncertainties are dominant in the region of forward-
rapidity, which is driven by PDF sampling in the region of small-$x$. An
improved description of PDFs in this region has important consequences for
neutrino astronomy [119, 120, 121, 122, 123, 124, 125], and may also provide a
cross-check of those results which have been obtained using forward $D$\- and
$B$-hadron production data [126, 119, 127, 124]. It is therefore recommended
that the experiment publishes a correlation matrix for the rapidity
distributions which also includes the rapidity-integrated distribution as an
entry (for a given $m_{\ell\bar{\ell}}$ region). As a final observation, the
NNLO correction for this quantity is $\approx 5\%$ at large
$y_{{\mathrm{\ell\bar{\ell}}}}$ which may indicate the contribution of large
$\ln[x]$ corrections. It could be interesting to investigate the impact of
resumming these corrections, such as in [128, 129] for the DIS process, using
the formalism presented in [130, 131, 132, 133].

Figure 4: Differential prediction for the normalised
$y_{{\mathrm{\ell\bar{\ell}}}}$ distribution within the LHCb fiducial region
for $m_{\ell\bar{\ell}}\in[12.5,13.5]\leavevmode\nobreak\ \mathrm{GeV}$.

The theoretical study of the $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$
distribution is a little more delicate (particularly at large
$m_{\ell\bar{\ell}}$) as a reliable description of the kinematic region of
small $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ relies on the resummation of
Sudakov logarithms of the form
$\frac{1}{p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}^{2}}\ln^{n}[p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}/m_{{\mathrm{\ell\bar{\ell}}}}]$.
The situation is tricky because the massive power-corrections obtained from
applying Eq. (2) contain contributions which have the same form as this, and
diverge in the limit $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}\to 0$. This
feature prohibits a straightforward matching of the fixed-order M-VFNS
prediction with a (Sudakov) resummed calculation (and potentially also joint
small-$x$ resummation [134]). This issue has been addressed and will be
detailed in future work, where a dedicated study of the impact of heavy-quark
mass effects on the transverse-momentum distributions of gauge-bosons will be
presented.

In the region of small $m_{\ell\bar{\ell}}$ the fixed-order results are likely
sufficient, and certainly useful to indicate the phenomenological
(ir)relevance of the massive power-corrections. These results are shown in
Fig. 5 for $m_{\ell\bar{\ell}}\in[12.5,13.5]\leavevmode\nobreak\
\mathrm{GeV}$, where the $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$
distribution is shown normalised to the integrated cross-section. The impact
of the massive power-corrections can be inferred by comparing the central
prediction of the M-VFNS (dash-dotted red) compared to that of the massless
calculation (dashed blue). The corrections amount to $\approx 1\%$ at small
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$, leading to a slight change in the
slope of the normalised $p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$
distribution. Overall, these effects are small compared to the either PDF or
scale uncertainties (which were, for visual clarity, not shown in the lower
panel).

Figure 5: Differential prediction for the normalised
$p_{\mathrm{T},{\mathrm{\ell\bar{\ell}}}}$ distribution within the LHCb
fiducial region for $m_{\ell\bar{\ell}}\in[12.5,13.5]\leavevmode\nobreak\
\mathrm{GeV}$.

##  9 Conclusions

The main goal of this work is to provide a deeper theoretical understanding of
the treatment and role of massive quarks in predicting hadron-hadron
scattering processes.

This has been achieved by studying the general structure of calculations which
involve a single massive quark, and presenting a formalism to construct
differential cross-section predictions in a massive variable flavour number
scheme. The formalism can be applied to colour-singlet production processes as
well as those involving (flavoured) hadronic jets, provided the differential
observables are inclusive with respect to QCD radiation and/or are infrared
and collinear safe. Hopefully, these developments will help to clarify several
issues regarding heavy-quark mass effects in hadron-hadron scattering
processes.

As a practical application, results have been provided for the low-mass DY
rapidity distribution within the LHCb fiducial region, and I have demonstrated
how a normalised distribution can give important information on the structure
of the proton at low-x.

Finally, the formalism presented here represents an important step towards
constructing a scheme which can be applied to predict the transverse-momentum
distribution of gauge bosons, that includes the impact of heavy-flavour
massive power-corrections and a resummation of Sudakov logarithoms in a
consistent way. This is a critical development towards reducing the theory
systematic related to the modelling of the $p_{\mathrm{T},\ell}$ distribution
in the charged-current DY process, which will in turn improve the sensitivity
of LHC measurements to extract the $W$ boson mass [135].

## Acknowledgements

I am grateful to Eric Laenen for comments at the initial stages of this
project. The work presented here has benefited from previous/on-going
collaboration with Adrian Rodrigues Garcia, Aude Gehrmann De–Ridder, Thomas
Gehrmann, Nigel Glover, and Alexander Huss as well as Marco Bonvini, Tommaso
Giani and Simone Marzani. I additionally thank Valerio Bertone, Davide
Napoletano, and Ben Pecjak for previous discussions about heavy-quark mass
effects, and Stephen Farry and Phil Ilten for comments/correspondence about
the LHCb measurement. A big thanks to Alex for providing constructive comments
on this manuscript, and to the careful eyes of Lucian Harland-Lang.

### Funding information

This research is supported by the Dutch Organisation for Scientific Research
(NWO) through the VENI grant 680-47-461.

## References

  * [1] J. C. Collins, D. E. Soper and G. F. Sterman,  _Factorization of Hard Processes in QCD_ ,  Adv. Ser. Direct. High Energy Phys. 5, 1 (1989),  10.1142/9789814503266_0001,  hep-ph/0409313.
  * [2] S. D. Drell and T.-M. Yan,  _Massive Lepton Pair Production in Hadron-Hadron Collisions at High-Energies_ ,  Phys. Rev. Lett. 25, 316 (1970),  10.1103/PhysRevLett.25.316,  [Erratum: Phys.Rev.Lett. 25, 902 (1970)].
  * [3] J. C. Collins, F. Wilczek and A. Zee,  _Low-Energy Manifestations of Heavy Particles: Application to the Neutral Current_ ,  Phys. Rev. D18, 242 (1978),  10.1103/PhysRevD.18.242.
  * [4] S. Kretzer and I. Schienbein,  _Heavy quark initiated contributions to deep inelastic structure functions_ ,  Phys. Rev. D 58, 094035 (1998),  10.1103/PhysRevD.58.094035,  hep-ph/9805233.
  * [5] J. C. Collins,  _Hard scattering factorization with heavy quarks: A General treatment_ ,  Phys. Rev. D58, 094002 (1998),  10.1103/PhysRevD.58.094002,  hep-ph/9806259.
  * [6] M. Cacciari, M. Greco and P. Nason,  _The P(T) spectrum in heavy flavor hadroproduction_ ,  JHEP 9805, 007 (1998),  10.1088/1126-6708/1998/05/007,  hep-ph/9803400.
  * [7] W.-K. Tung, S. Kretzer and C. Schmidt,  _Open heavy flavor production in QCD: Conceptual framework and implementation issues_ ,  J. Phys. G28, 983 (2002),  10.1088/0954-3899/28/5/321,  hep-ph/0110247.
  * [8] B. A. Kniehl, G. Kramer, I. Schienbein and H. Spiesberger,  _Collinear subtractions in hadroproduction of heavy quarks_ ,  Eur. Phys. J. C41, 199 (2005),  10.1140/epjc/s2005-02200-7,  hep-ph/0502194.
  * [9] A. Mitov and S. Moch,  _The Singular behavior of massive QCD amplitudes_ ,  JHEP 05, 001 (2007),  10.1088/1126-6708/2007/05/001,  hep-ph/0612149.
  * [10] I. Bierenbaum, J. Blumlein and S. Klein,  _The Gluonic Operator Matrix Elements at O(alpha(s)**2) for DIS Heavy Flavor Production_ ,  Phys. Lett. B 672, 401 (2009),  10.1016/j.physletb.2009.01.057,  0901.0669.
  * [11] M. Guzzi, P. M. Nadolsky, H.-L. Lai and C.-P. Yuan,  _General-Mass Treatment for Deep Inelastic Scattering at Two-Loop Accuracy_ ,  Phys.Rev. D86, 053005 (2012),  10.1103/PhysRevD.86.053005,  1108.5112.
  * [12] R. Harlander, M. Kramer and M. Schumacher,  _Bottom-quark associated Higgs-boson production: reconciling the four- and five-flavour scheme approach_ (2011),  1112.3478.
  * [13] F. Maltoni, G. Ridolfi and M. Ubiali,  _b-initiated processes at the LHC: a reappraisal_ ,  JHEP 07, 022 (2012),  10.1007/JHEP04(2013)095, 10.1007/JHEP07(2012)022,  [Erratum: JHEP04,095(2013)],  1203.6393.
  * [14] A. Behring, I. Bierenbaum, J. Blümlein, A. De Freitas, S. Klein and F. Wißbrock,  _The logarithmic contributions to the $O(\alpha^{3}_{s})$ asymptotic massive Wilson coefficients and operator matrix elements in deeply inelastic scattering_,  Eur. Phys. J. C 74(9), 3033 (2014),  10.1140/epjc/s10052-014-3033-x,  1403.6356.
  * [15] M. Bonvini, A. S. Papanastasiou and F. J. Tackmann,  _Resummation and matching of b-quark mass effects in $b\overline{b}H$ production_,  JHEP 11, 196 (2015),  10.1007/JHEP11(2015)196,  1508.03288.
  * [16] A. H. Hoang, P. Pietrulewicz and D. Samitz,  _Variable Flavor Number Scheme for Final State Jets in DIS_ ,  Phys. Rev. D93(3), 034034 (2016),  10.1103/PhysRevD.93.034034,  1508.04323.
  * [17] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, C. Schneider and F. Wißbrock,  _Three Loop Massive Operator Matrix Elements and Asymptotic Wilson Coefficients with Two Different Masses_ ,  Nucl. Phys. B 921, 585 (2017),  10.1016/j.nuclphysb.2017.05.017,  1705.07030.
  * [18] F. Krauss and D. Napoletano,  _Towards a fully massive five-flavor scheme_ ,  Phys. Rev. D 98(9), 096002 (2018),  10.1103/PhysRevD.98.096002,  1712.06832.
  * [19] V. Bertone, A. Glazov, A. Mitov, A. Papanastasiou and M. Ubiali,  _Heavy-flavor parton distributions without heavy-flavor matching prescriptions_ ,  JHEP 04, 046 (2018),  10.1007/JHEP04(2018)046,  1711.03355.
  * [20] I. Helenius and H. Paukkunen,  _Revisiting the D-meson hadroproduction in general-mass variable flavour number scheme_ ,  JHEP 05, 196 (2018),  10.1007/JHEP05(2018)196,  1804.03557.
  * [21] S. Forte, T. Giani and D. Napoletano,  _Fitting the b-quark PDF as a massive-b scheme: Higgs production in bottom fusion_ ,  Eur. Phys. J. C 79(7), 609 (2019),  10.1140/epjc/s10052-019-7119-3,  1905.02207.
  * [22] M. A. G. Aivazis, F. I. Olness and W.-K. Tung,  _Leptoproduction of heavy quarks. 1. General formalism and kinematics of charged current and neutral current production processes_ ,  Phys. Rev. D 50, 3085 (1994),  10.1103/PhysRevD.50.3085,  hep-ph/9312318.
  * [23] M. A. G. Aivazis, J. C. Collins, F. I. Olness and W.-K. Tung,  _Leptoproduction of heavy quarks. 2. A Unified QCD formulation of charged and neutral current processes from fixed target to collider energies_ ,  Phys. Rev. D50, 3102 (1994),  10.1103/PhysRevD.50.3102,  hep-ph/9312319.
  * [24] M. Buza, Y. Matiounine, J. Smith and W. L. van Neerven,  _Charm electroproduction viewed in the variable flavor number scheme versus fixed order perturbation theory_ ,  Eur. Phys. J. C1, 301 (1998),  10.1007/BF01245820,  hep-ph/9612398.
  * [25] J. C. Collins,  _Proof of factorization for diffractive hard scattering_ ,  Phys. Rev. D 57, 3051 (1998),  10.1103/PhysRevD.61.019902,  [Erratum: Phys.Rev.D 61, 019902 (2000)],  hep-ph/9709499.
  * [26] R. S. Thorne and R. G. Roberts,  _An Ordered analysis of heavy flavor production in deep inelastic scattering_ ,  Phys. Rev. D57, 6871 (1998),  10.1103/PhysRevD.57.6871,  hep-ph/9709442.
  * [27] M. Kramer, F. I. Olness and D. E. Soper,  _Treatment of heavy quarks in deeply inelastic scattering_ ,  Phys. Rev. D62, 096007 (2000),  10.1103/PhysRevD.62.096007,  hep-ph/0003035.
  * [28] R. S. Thorne,  _A Variable-flavor number scheme for NNLO_ ,  Phys. Rev. D73, 054019 (2006),  10.1103/PhysRevD.73.054019,  hep-ph/0601245.
  * [29] S. Forte, E. Laenen, P. Nason and J. Rojo,  _Heavy quarks in deep-inelastic scattering_ ,  Nucl. Phys. B834, 116 (2010),  10.1016/j.nuclphysb.2010.03.014,  1001.2312.
  * [30] R. D. Ball, M. Bonvini and L. Rottoli,  _Charm in Deep-Inelastic Scattering_ ,  JHEP 11, 122 (2015),  10.1007/JHEP11(2015)122,  1510.02491.
  * [31] S. Alekhin, J. Blümlein and S. Moch,  _Heavy-flavor PDF evolution and variable-flavor number scheme uncertainties in deep-inelastic scattering_ ,  Phys. Rev. D 102(5), 054014 (2020),  10.1103/PhysRevD.102.054014,  2006.07032.
  * [32] M. Cacciari, S. Frixione and P. Nason,  _The p(T) spectrum in heavy flavor photoproduction_ ,  JHEP 0103, 006 (2001),  10.1088/1126-6708/2001/03/006,  hep-ph/0102134.
  * [33] A. Banfi, G. P. Salam and G. Zanderighi,  _Accurate QCD predictions for heavy-quark jets at the Tevatron and LHC_ ,  JHEP 07, 026 (2007),  10.1088/1126-6708/2007/07/026,  0704.2999.
  * [34] D. Figueroa, S. Honeywell, S. Quackenbush, L. Reina, C. Reuschle and D. Wackeroth,  _Electroweak and QCD corrections to $Z$-boson production with one $b$ jet in a massive five-flavor scheme_,  Phys. Rev. D 98(9), 093002 (2018),  10.1103/PhysRevD.98.093002,  1805.01353.
  * [35] R. Gauld, A. Gehrmann-De Ridder, E. W. N. Glover, A. Huss and I. Majer,  _Predictions for $Z$ -Boson Production in Association with a $b$-Jet at $\mathcal{O}(\alpha_{s}^{3})$_,  Phys. Rev. Lett. 125(22), 222002 (2020),  10.1103/PhysRevLett.125.222002,  2005.03016.
  * [36] S. Forte, D. Napoletano and M. Ubiali,  _Higgs production in bottom-quark fusion in a matched scheme_ ,  Phys. Lett. B751, 331 (2015),  10.1016/j.physletb.2015.10.051,  1508.01529.
  * [37] S. Forte, D. Napoletano and M. Ubiali,  _Higgs production in bottom-quark fusion: matching beyond leading order_ ,  Phys. Lett. B763, 190 (2016),  10.1016/j.physletb.2016.10.040,  1607.00389.
  * [38] S. Forte, D. Napoletano and M. Ubiali,  _$Z$ boson production in bottom-quark fusion: a study of $b$-mass effects beyond leading order_,  Eur. Phys. J. C78(11), 932 (2018),  10.1140/epjc/s10052-018-6414-8,  1803.10248.
  * [39] 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_ (2020),  2004.04752.
  * [40] P. Pietrulewicz, D. Samitz, A. Spiering and F. J. Tackmann,  _Factorization and Resummation for Massive Quark Effects in Exclusive Drell-Yan_ ,  JHEP 08, 114 (2017),  10.1007/JHEP08(2017)114,  1703.09702.
  * [41] E. Bagnaschi, F. Maltoni, A. Vicini and M. Zaro,  _Lepton-pair production in association with a $b\overline{b}$ pair and the determination of the $W$ boson mass_,  JHEP 07, 101 (2018),  10.1007/JHEP07(2018)101,  1803.04336.
  * [42] S. Hoche, J. Krause and F. Siegert,  _Multijet Merging in a Variable Flavor Number Scheme_ ,  Phys. Rev. D100(1), 014011 (2019),  10.1103/PhysRevD.100.014011,  1904.09382.
  * [43] G. Aad _et al._ ,  _Measurement of the transverse momentum distribution of Drell–Yan lepton pairs in proton–proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_,  Eur. Phys. J. C 80(7), 616 (2020),  10.1140/epjc/s10052-020-8001-z,  1912.02844.
  * [44] C. Duhr, F. Dulat and B. Mistlberger,  _Drell-Yan Cross Section to Third Order in the Strong Coupling Constant_ ,  Phys. Rev. Lett. 125(17), 172001 (2020),  10.1103/PhysRevLett.125.172001,  2001.07717.
  * [45] C. Duhr, F. Dulat and B. Mistlberger,  _Charged current Drell-Yan production at N 3LO_,  JHEP 11, 143 (2020),  10.1007/JHEP11(2020)143,  2007.13313.
  * [46] S. Camarda, L. Cieri and G. Ferrera,  _Drell-Yan lepton-pair production: $q_{T}$ resummation at N3LL accuracy and fiducial cross sections at N3LO_ (2021),  2103.04974.
  * [47] A. Gehrmann-De Ridder, T. Gehrmann, E. W. N. Glover, A. Huss and T. A. Morgan,  _Precise QCD predictions for the production of a Z boson in association with a hadronic jet_ ,  Phys. Rev. Lett. 117(2), 022001 (2016),  10.1103/PhysRevLett.117.022001,  1507.02850.
  * [48] R. Boughezal, J. M. Campbell, R. K. Ellis, C. Focke, W. T. Giele, X. Liu and F. Petriello,  _Z-boson production in association with a jet at next-to-next-to-leading order in perturbative QCD_ ,  Phys. Rev. Lett. 116(15), 152001 (2016),  10.1103/PhysRevLett.116.152001,  1512.01291.
  * [49] L. Cieri, G. Ferrera and G. F. R. Sborlini,  _Combining QED and QCD transverse-momentum resummation for Z boson production at hadron colliders_ ,  JHEP 08, 165 (2018),  10.1007/JHEP08(2018)165,  1805.11948.
  * [50] W. Bizoń, X. Chen, A. Gehrmann-De Ridder, T. Gehrmann, N. Glover, A. Huss, P. F. Monni, E. Re, L. Rottoli and P. Torrielli,  _Fiducial distributions in Higgs and Drell-Yan production at N 3LL+NNLO_,  JHEP 12, 132 (2018),  10.1007/JHEP12(2018)132,  1805.05916.
  * [51] T. Becher and M. Hager,  _Event-Based Transverse Momentum Resummation_ ,  Eur. Phys. J. C 79(8), 665 (2019),  10.1140/epjc/s10052-019-7136-2,  1904.08325.
  * [52] W. Bizon, A. Gehrmann-De Ridder, T. Gehrmann, N. Glover, A. Huss, P. F. Monni, E. Re, L. Rottoli and D. M. Walker,  _The transverse momentum spectrum of weak gauge bosons at N 3 LL + NNLO_,  Eur. Phys. J. C 79(10), 868 (2019),  10.1140/epjc/s10052-019-7324-0,  1905.05171.
  * [53] A. Bacchetta, V. Bertone, C. Bissolotti, G. Bozzi, F. Delcarro, F. Piacenza and M. Radici,  _Transverse-momentum-dependent parton distributions up to N 3LL from Drell-Yan data_,  JHEP 07, 117 (2020),  10.1007/JHEP07(2020)117,  1912.07550.
  * [54] M. A. Ebert, J. K. L. Michel, I. W. Stewart and F. J. Tackmann,  _Drell-Yan $q_{T}$ resummation of fiducial power corrections at N3LL_,  JHEP 04, 102 (2021),  10.1007/JHEP04(2021)102,  2006.11382.
  * [55] E. Re, L. Rottoli and P. Torrielli,  _Fiducial Higgs and Drell-Yan distributions at N 3LL′+NNLO with RadISH_ (2021),  2104.07509.
  * [56] A. Banfi, G. P. Salam and G. Zanderighi,  _Infrared safe definition of jet flavor_ ,  Eur. Phys. J. C47, 113 (2006),  10.1140/epjc/s2006-02552-4,  hep-ph/0601139.
  * [57] A. Behring, W. Bizoń, F. Caola, K. Melnikov and R. Röntsch,  _Bottom quark mass effects in associated $WH$ production with the $H\to b\bar{b}$ decay through NNLO QCD_,  Phys. Rev. D 101(11), 114012 (2020),  10.1103/PhysRevD.101.114012,  2003.08321.
  * [58] K. G. Chetyrkin, B. A. Kniehl and M. Steinhauser,  _Strong coupling constant with flavor thresholds at four loops in the MS scheme_ ,  Phys. Rev. Lett. 79, 2184 (1997),  10.1103/PhysRevLett.79.2184,  hep-ph/9706430.
  * [59] K. G. Chetyrkin, J. H. Kuhn and M. Steinhauser,  _RunDec: A Mathematica package for running and decoupling of the strong coupling and quark masses_ ,  Comput. Phys. Commun. 133, 43 (2000),  10.1016/S0010-4655(00)00155-7,  hep-ph/0004189.
  * [60] M. Buza, Y. Matiounine, J. Smith, R. Migneron and W. L. van Neerven,  _Heavy quark coefficient functions at asymptotic values $Q\leavevmode\nobreak\ 2\gg m\leavevmode\nobreak\ 2$_,  Nucl. Phys. B472, 611 (1996),  10.1016/0550-3213(96)00228-3,  hep-ph/9601302.
  * [61] I. Bierenbaum, J. Blumlein and S. Klein,  _Two-Loop Massive Operator Matrix Elements and Unpolarized Heavy Flavor Production at Asymptotic Values Q**2 $>$$>$ m**2_,  Nucl. Phys. B 780, 40 (2007),  10.1016/j.nuclphysb.2007.04.030,  hep-ph/0703285.
  * [62] M. Buza, Y. Matiounine, J. Smith and W. L. van Neerven,  _O (alpha-s**2) corrections to polarized heavy flavor production at Q**2 is $>$$>$ m**2_,  Nucl. Phys. B 485, 420 (1997),  10.1016/S0550-3213(96)00606-2,  hep-ph/9608342.
  * [63] I. Bierenbaum, J. Blumlein and S. Klein,  _Two-loop massive operator matrix elements for polarized and unpolarized deep-inelastic scattering_ ,  In _15th International Workshop on Deep-Inelastic Scattering and Related Subjects_ ,  10.3204/proc07-01/148 (2007), 0706.2738.
  * [64] M. Buza and W. L. van Neerven,  _O(alpha(s)**2) contributions to charm production in charged-current deep-inelastic lepton hadron scattering_ ,  Nucl. Phys. B500, 301 (1997),  10.1016/S0550-3213(97)00327-1,  hep-ph/9702242.
  * [65] J. Blümlein, A. Hasselhuhn and T. Pfoh,  _The $O(\alpha_{s}^{2})$ heavy quark corrections to charged current deep-inelastic scattering at large virtualities_,  Nucl. Phys. B 881, 1 (2014),  10.1016/j.nuclphysb.2014.01.023,  1401.4352.
  * [66] J. Ablinger, J. Blumlein, S. Klein, C. Schneider and F. Wissbrock,  _The $O(\alpha_{s}^{3})$ Massive Operator Matrix Elements of $O(n_{f})$ for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_,  Nucl. Phys. B 844, 26 (2011),  10.1016/j.nuclphysb.2010.10.021,  1008.3347.
  * [67] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round, C. Schneider and F. Wißbrock,  _The 3-Loop Non-Singlet Heavy Flavor Contributions and Anomalous Dimensions for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_,  Nucl. Phys. B 886, 733 (2014),  10.1016/j.nuclphysb.2014.07.010,  1406.4654.
  * [68] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, C. Schneider and K. Schönwald,  _The three-loop single mass polarized pure singlet operator matrix element_ ,  Nucl. Phys. B 953, 114945 (2020),  10.1016/j.nuclphysb.2020.114945,  1912.02536.
  * [69] A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, K. Schönwald and C. Schneider,  _The polarized transition matrix element $A_{gq}(N)$ of the variable flavor number scheme at $O(\alpha^{3}_{s})$_,  Nucl. Phys. B 964, 115331 (2021),  10.1016/j.nuclphysb.2021.115331,  2101.05733.
  * [70] J. Blümlein, A. De Freitas, M. Saragnese, C. Schneider and K. Schönwald,  _The Logarithmic Contributions to the Polarized $O(\alpha^{3}_{s})$ Asymptotic Massive Wilson Coefficients and Operator Matrix Elements in Deeply Inelastic Scattering_ (2021),  2105.09572.
  * [71] V. Hirschi, R. Frederix, S. Frixione, M. V. Garzelli, F. Maltoni and R. Pittau,  _Automation of one-loop QCD corrections_ ,  JHEP 05, 044 (2011),  10.1007/JHEP05(2011)044,  1103.0621.
  * [72] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni _et al._ ,  _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ ,  JHEP 1407, 079 (2014),  10.1007/JHEP07(2014)079,  1405.0301.
  * [73] R. D. Ball _et al._ ,  _Parton distributions from high-precision collider data_ ,  Eur. Phys. J. C77(10), 663 (2017),  10.1140/epjc/s10052-017-5199-5,  1706.00428.
  * [74] T.-J. Hou, S. Dulat, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, C. Schmidt, J. Winter, K. Xie and C. P. Yuan,  _CT14 Intrinsic Charm Parton Distribution Functions from CTEQ-TEA Global Analysis_ ,  JHEP 02, 059 (2018),  10.1007/JHEP02(2018)059,  1707.00657.
  * [75] R. D. Ball, V. Bertone, M. Bonvini, S. Forte, P. Groth Merrild, J. Rojo and L. Rottoli,  _Intrinsic charm in a matched general-mass scheme_ ,  Phys. Lett. B 754, 49 (2016),  10.1016/j.physletb.2015.12.077,  1510.00009.
  * [76] R. Doria, J. Frenkel and J. C. Taylor,  _Counter Example to Nonabelian Bloch-Nordsieck Theorem_ ,  Nucl. Phys. B 168, 93 (1980),  10.1016/0550-3213(80)90278-3.
  * [77] C. Di’Lieto, S. Gendron, I. G. Halliday and C. T. Sachrajda,  _A Counter Example to the Bloch-Nordsieck Theorem in Nonabelian Gauge Theories_ ,  Nucl. Phys. B 183, 223 (1981),  10.1016/0550-3213(81)90554-X.
  * [78] J. Frenkel, J. G. M. Gatheral and J. C. Taylor,  _IS QUARK - ANTI-QUARK ANNIHILATION INFRARED SAFE AT HIGH-ENERGY?_ ,  Nucl. Phys. B 233, 307 (1984),  10.1016/0550-3213(84)90418-8.
  * [79] S. Catani, M. Ciafaloni and G. Marchesini,  _NONCANCELLING INFRARED DIVERGENCES IN QCD COHERENT STATE_ ,  Nucl. Phys. B 264, 588 (1986),  10.1016/0550-3213(86)90500-6.
  * [80] S. Catani,  _Violation of the Bloch-nordsieck Mechanism in General Nonabelian Theories and SUSY QCD_ ,  Z. Phys. C 37, 357 (1988),  10.1007/BF01578128.
  * [81] A. Andrasi, M. Day, R. Doria, J. Frenkel and J. C. Taylor,  _Soft Divergences in Perturbative QCD_ ,  Nucl. Phys. B 182, 104 (1981),  10.1016/0550-3213(81)90460-0.
  * [82] C. A. Nelson,  _Avoidance of Counter Example to Nonabelian Bloch-Nordsieck Conjecture by Using Coherent State Approach_ ,  Nucl. Phys. B 186, 187 (1981),  10.1016/0550-3213(81)90099-7.
  * [83] I. Ito,  _Cancellation of Infrared Divergence and Initial Degenerate State in QCD_ ,  Prog. Theor. Phys. 65, 1466 (1981),  10.1143/PTP.65.1466.
  * [84] N. Yoshida,  _Diagrammatical Display of the Counter Example to Nonabelian Bloch-nordsieck Conjecture_ ,  Prog. Theor. Phys. 66, 269 (1981),  10.1143/PTP.66.269.
  * [85] N. Yoshida,  _CANCELLATION OF THE INFRARED SINGULARITY THROUGH THE UNITARITY RELATION_ ,  Prog. Theor. Phys. 66, 1803 (1981),  10.1143/PTP.66.1803.
  * [86] T. Muta and C. A. Nelson,  _Role of Quark - Gluon Degenerate States in Perturbative QCD_ ,  Phys. Rev. D 25, 2222 (1982),  10.1103/PhysRevD.25.2222.
  * [87] B. F. L. Ward,  _Quark masses and resummation in precision QCD theory_ ,  Phys. Rev. D 78, 056001 (2008),  10.1103/PhysRevD.78.056001,  0707.2101.
  * [88] F. Caola, K. Melnikov, D. Napoletano and L. Tancredi,  _Noncancellation of infrared singularities in collisions of massive quarks_ ,  Phys. Rev. D 103(5), 054013 (2021),  10.1103/PhysRevD.103.054013,  2011.04701.
  * [89] D. A. Dicus and S. S. D. Willenbrock,  _Radiative Corrections to the Ratio of $Z$ and $W$ Boson Production_,  Phys. Rev. D 34, 148 (1986),  10.1103/PhysRevD.34.148.
  * [90] 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, 4930 (1992),  10.1103/PhysRevD.46.4930.
  * [91] P. J. Rijken and W. L. van Neerven,  _Heavy flavor contributions to the Drell-Yan cross-section_ ,  Phys. Rev. D 52, 149 (1995),  10.1103/PhysRevD.52.149,  hep-ph/9501373.
  * [92] S. Catani and M. H. Seymour,  _A General algorithm for calculating jet cross-sections in NLO QCD_ ,  Nucl. Phys. B485, 291 (1997),  10.1016/S0550-3213(96)00589-5, 10.1016/S0550-3213(98)81022-5,  [Erratum: Nucl. Phys.B510,503(1998)],  hep-ph/9605323.
  * [93] S. Dittmaier,  _A General approach to photon radiation off fermions_ ,  Nucl. Phys. B 565, 69 (2000),  10.1016/S0550-3213(99)00563-5,  hep-ph/9904440.
  * [94] J. Gaunt, M. Stahlhofen, F. J. Tackmann and J. R. Walsh,  _N-jettiness Subtractions for NNLO QCD Calculations_ ,  JHEP 09, 058 (2015),  10.1007/JHEP09(2015)058,  1505.04794.
  * [95] B. A. Kniehl,  _Two Loop QED Vertex Correction From Virtual Heavy Fermions_ ,  Phys. Lett. B 237, 127 (1990),  10.1016/0370-2693(90)90474-K.
  * [96] Z. Bern, L. J. Dixon and D. A. Kosower,  _One loop amplitudes for e+ e- to four partons_ ,  Nucl. Phys. B 513, 3 (1998),  10.1016/S0550-3213(97)00703-7,  hep-ph/9708239.
  * [97] 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, 094 (2010),  10.1007/JHEP06(2010)094,  1004.3653.
  * [98] J. Blümlein, G. Falcioni and A. De Freitas,  _The Complete $O(\alpha_{s}^{2})$ Non-Singlet Heavy Flavor Corrections to the Structure Functions $g_{1,2}^{ep}(x,Q^{2})$, $F_{1,2,L}^{ep}(x,Q^{2})$, $F_{1,2,3}^{\nu(\bar{\nu})}(x,Q^{2})$ and the Associated Sum Rules_,  Nucl. Phys. B 910, 568 (2016),  10.1016/j.nuclphysb.2016.06.018,  1605.05541.
  * [99] C. W. Bauer, S. Fleming, D. Pirjol and I. W. Stewart,  _An Effective field theory for collinear and soft gluons: Heavy to light decays_ ,  Phys. Rev. D 63, 114020 (2001),  10.1103/PhysRevD.63.114020,  hep-ph/0011336.
  * [100] I. W. Stewart, F. J. Tackmann and W. J. Waalewijn,  _Factorization at the LHC: From PDFs to Initial State Jets_ ,  Phys. Rev. D 81, 094035 (2010),  10.1103/PhysRevD.81.094035,  0910.0467.
  * [101] I. W. Stewart, F. J. Tackmann and W. J. Waalewijn,  _N-Jettiness: An Inclusive Event Shape to Veto Jets_ ,  Phys. Rev. Lett. 105, 092002 (2010),  10.1103/PhysRevLett.105.092002,  1004.2489.
  * [102] R. Kelley, M. D. Schwartz, R. M. Schabinger and H. X. Zhu,  _The two-loop hemisphere soft function_ ,  Phys. Rev. D 84, 045022 (2011),  10.1103/PhysRevD.84.045022,  1105.3676.
  * [103] P. F. Monni, T. Gehrmann and G. Luisoni,  _Two-Loop Soft Corrections and Resummation of the Thrust Distribution in the Dijet Region_ ,  JHEP 08, 010 (2011),  10.1007/JHEP08(2011)010,  1105.4560.
  * [104] A. Hornig, C. Lee, I. W. Stewart, J. R. Walsh and S. Zuberi,  _Non-global Structure of the $\mathcal{O}({\alpha}^{2}_{s})$ Dijet Soft Function_,  JHEP 08, 054 (2011),  10.1007/JHEP08(2011)054,  [Erratum: JHEP 10, 101 (2017)],  1105.4628.
  * [105] J. R. Gaunt, M. Stahlhofen and F. J. Tackmann,  _The Quark Beam Function at Two Loops_ ,  JHEP 04, 113 (2014),  10.1007/JHEP04(2014)113,  1401.5478.
  * [106] I. Moult, L. Rothen, I. W. Stewart, F. J. Tackmann and H. X. Zhu,  _Subleading Power Corrections for N-Jettiness Subtractions_ ,  Phys. Rev. D 95(7), 074023 (2017),  10.1103/PhysRevD.95.074023,  1612.00450.
  * [107] F. Buccioni, J.-N. Lang, J. M. Lindert, P. Maierhöfer, S. Pozzorini, H. Zhang and M. F. Zoller,  _OpenLoops 2_ ,  Eur. Phys. J. C 79(10), 866 (2019),  10.1140/epjc/s10052-019-7306-2,  1907.13071.
  * [108] T. Hahn,  _Generating Feynman diagrams and amplitudes with FeynArts 3_ ,  Comput. Phys. Commun. 140, 418 (2001),  10.1016/S0010-4655(01)00290-9,  hep-ph/0012260.
  * [109] T. Hahn and M. Perez-Victoria,  _Automatized one loop calculations in four-dimensions and D-dimensions_ ,  Comput. Phys. Commun. 118, 153 (1999),  10.1016/S0010-4655(98)00173-8,  hep-ph/9807565.
  * [110] S. Catani, D. de Florian and M. Grazzini,  _Universality of nonleading logarithmic contributions in transverse momentum distributions_ ,  Nucl. Phys. B 596, 299 (2001),  10.1016/S0550-3213(00)00617-9,  hep-ph/0008184.
  * [111] A. Banfi, M. Dasgupta, S. Marzani and L. Tomlinson,  _Probing the low transverse momentum domain of Z production with novel variables_ ,  JHEP 01, 044 (2012),  10.1007/JHEP01(2012)044,  1110.4009.
  * [112] S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini,  _Vector boson production at hadron colliders: hard-collinear coefficients at the NNLO_ ,  Eur. Phys. J. C 72, 2195 (2012),  10.1140/epjc/s10052-012-2195-7,  1209.0158.
  * [113] T. Hahn,  _CUBA: A Library for multidimensional numerical integration_ ,  Comput. Phys. Commun. 168, 78 (2005),  10.1016/j.cpc.2005.01.010,  hep-ph/0404043.
  * [114] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page _et al._ ,  _LHAPDF6: parton density access in the LHC precision era_ ,  Eur.Phys.J. C75(3), 132 (2015),  10.1140/epjc/s10052-015-3318-8,  1412.7420.
  * [115] A. Denner, S. Dittmaier, M. Roth and L. H. Wieders,  _Electroweak corrections to charged-current e+ e- — $>$ 4 fermion processes: Technical details and further results_,  Nucl. Phys. B 724, 247 (2005),  10.1016/j.nuclphysb.2011.09.001,  [Erratum: Nucl.Phys.B 854, 504–507 (2012)],  hep-ph/0505042.
  * [116] S. Alekhin, A. Kardos, S. Moch and Z. Trócsányi,  _Precision studies for Drell-Yan processes at NNLO_ (2021),  2104.02400.
  * [117] T. Kinoshita,  _Mass singularities of Feynman amplitudes_ ,  J. Math. Phys. 3, 650 (1962),  10.1063/1.1724268.
  * [118] T. D. Lee and M. Nauenberg,  _Degenerate Systems and Mass Singularities_ ,  Phys. Rev. 133, B1549 (1964),  10.1103/PhysRev.133.B1549.
  * [119] R. Gauld, J. Rojo, L. Rottoli and J. Talbert,  _Charm production in the forward region: constraints on the small-x gluon and backgrounds for neutrino astronomy_ ,  JHEP 11, 009 (2015),  10.1007/JHEP11(2015)009,  1506.08025.
  * [120] R. Gauld, J. Rojo, L. Rottoli, S. Sarkar and J. Talbert,  _The prompt atmospheric neutrino flux in the light of LHCb_ ,  JHEP 02, 130 (2016),  10.1007/JHEP02(2016)130,  1511.06346.
  * [121] R. Gauld and J. Rojo,  _Precision determination of the small- $x$ gluon from charm production at LHCb_,  Phys. Rev. Lett. 118(7), 072001 (2017),  10.1103/PhysRevLett.118.072001,  1610.09373.
  * [122] V. Bertone, R. Gauld and J. Rojo,  _Neutrino Telescopes as QCD Microscopes_ ,  JHEP 01, 217 (2019),  10.1007/JHEP01(2019)217,  1808.02034.
  * [123] R. Gauld,  _Precise predictions for multi-TeV and PeV energy neutrino scattering rates_ ,  Phys. Rev. D 100(9), 091301 (2019),  10.1103/PhysRevD.100.091301,  1905.03792.
  * [124] O. Zenaiev, M. V. Garzelli, K. Lipka, S. O. Moch, A. Cooper-Sarkar, F. Olness, A. Geiser and G. Sigl,  _Improved constraints on parton distributions using LHCb, ALICE and HERA heavy-flavour measurements and implications for the predictions for prompt atmospheric-neutrino fluxes_ ,  JHEP 04, 118 (2020),  10.1007/JHEP04(2020)118,  1911.13164.
  * [125] A. Garcia, R. Gauld, A. Heijboer and J. Rojo,  _Complete predictions for high-energy neutrino propagation in matter_ ,  JCAP 09, 025 (2020),  10.1088/1475-7516/2020/09/025,  2004.04756.
  * [126] O. Zenaiev _et al._ ,  _Impact of heavy-flavour production cross sections measured by the LHCb experiment on parton distribution functions at low x_ ,  Eur. Phys. J. C 75(8), 396 (2015),  10.1140/epjc/s10052-015-3618-z,  1503.04581.
  * [127] M. Cacciari, M. L. Mangano and P. Nason,  _Gluon PDF constraints from the ratio of forward heavy-quark production at the LHC at $\sqrt{S}=7$ and 13 TeV_,  Eur. Phys. J. C75(12), 610 (2015),  10.1140/epjc/s10052-015-3814-x,  1507.06197.
  * [128] R. D. Ball, V. Bertone, M. Bonvini, S. Marzani, J. Rojo and L. Rottoli,  _Parton distributions with small-x resummation: evidence for BFKL dynamics in HERA data_ ,  Eur. Phys. J. C 78(4), 321 (2018),  10.1140/epjc/s10052-018-5774-4,  1710.05935.
  * [129] H. Abdolmaleki _et al._ ,  _Impact of low- $x$ resummation on QCD analysis of HERA data_,  Eur. Phys. J. C 78(8), 621 (2018),  10.1140/epjc/s10052-018-6090-8,  1802.00064.
  * [130] M. Bonvini, S. Marzani and T. Peraro,  _Small- $x$ resummation from HELL_,  Eur. Phys. J. C 76(11), 597 (2016),  10.1140/epjc/s10052-016-4445-6,  1607.02153.
  * [131] M. Bonvini, S. Marzani and C. Muselli,  _Towards parton distribution functions with small- $x$ resummation: HELL 2.0_,  JHEP 12, 117 (2017),  10.1007/JHEP12(2017)117,  1708.07510.
  * [132] M. Bonvini and S. Marzani,  _Four-loop splitting functions at small $x$_,  JHEP 06, 145 (2018),  10.1007/JHEP06(2018)145,  1805.06460.
  * [133] M. Bonvini,  _Small- $x$ phenomenology at the LHC and beyond: HELL 3.0 and the case of the Higgs cross section_,  Eur. Phys. J. C 78(10), 834 (2018),  10.1140/epjc/s10052-018-6315-x,  1805.08785.
  * [134] S. Marzani,  _Combining $Q_{T}$ and small-$x$ resummations_,  Phys. Rev. D 93(5), 054047 (2016),  10.1103/PhysRevD.93.054047,  1511.06039.
  * [135] M. Aaboud _et al._ ,  _Measurement of the $W$-boson mass in pp collisions at $\sqrt{s}=7$ TeV with the ATLAS detector_,  Eur. Phys. J. C 78(2), 110 (2018),  10.1140/epjc/s10052-017-5475-4,  [Erratum: Eur.Phys.J.C 78, 898 (2018)],  1701.07240.