aainstitutetext: Université Paris-Saclay, CNRS, IJCLab, 91405 Orsay,
Francebbinstitutetext: National Centre for Nuclear Research (NCBJ), Pasteura
7, 02-093 Warsaw, Polandccinstitutetext: Institute for Theoretical Particle
Physics, KIT, 76128 Karlsruhe, Germany

# Matching next-to-leading-order and high-energy-resummed calculations of
heavy-quarkonium-hadroproduction cross sections

Jean-Philippe Lansberg 1,a,b Maxim Nefedov,111Corresponding author. a,c Melih
A. Ozcelik<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The energy dependence of the total hadroproduction cross section of pseudo-
scalar quarkonia is computed via matching Next-to-Leading Order (NLO)
Collinear-Factorisation (CF) results with resummed higher-order corrections,
proportional to $\alpha_{s}^{n}\ln^{n-1}(1/z)$, to the CF hard-scattering
coefficient, where $z=M^{2}/\hat{s}$ with $M$ and $\hat{s}$ being the
quarkonium mass and the partonic center-of-mass energy squared. The
resummation is performed using High-Energy Factorisation (HEF) in the Doubly-
Logarithmic (DL) approximation, which is a subset of the leading logarithmic
$\ln(1/z)$ approximation. Doing so, one remains strictly consistent with the
NLO and NNLO DGLAP evolution of the PDFs. By improving the treatment of the
small-$z$ asymptotics of the CF coefficient function, the resummation cures
the unphysical results of the NLO CF calculation. The matching is directly
performed in the $z$-space and, for the first time, by using the Inverse-Error
Weighting (InEW) matching procedure. As a by-product of the calculation, the
NNLO term of the CF hard-scattering coefficient proportional to
$\alpha_{s}^{2}\ln(1/z)$ is predicted from HEF.

††preprint: TTP-21-057

##  1 Introduction: heavy-quarkonium cross-section instability at high
energies

Historically, the discovery of charm quarks and charmonia has been one of the
most important milestones towards establishing QCD as the theory of strong
interaction. Nowadays, charmonia and bottomonia attract a lot of attention as
tools to study the proton structure, spin physics and/or to probe the quark-
gluon plasma, see e.g. recent reviews Chapon:2020heu ; Arbuzov:2020cqg ;
Lansberg:2019adr ; Brambilla:2014jmp ; Brambilla:2010cs . They have clean
experimental signatures and the fact that the hard scale of their production
process is limited from below by the heavy-quark mass justifies the
application of perturbation theory. However, the usage of quarkonia as tools
for precision studies is problematic because up to now there is no consensus
in the community on what are the main mechanisms of quarkonium production. The
Non-Relativistic QCD(NRQCD) factorisation approach Bodwin:1994jh  at NLO in
the $\alpha_{s}$ expansion and NNLO in the $v^{2}$ expansion, which seems to
be the most systematic theoretical tool available so far, is able to describe
the unpolarised $p_{T}$-differential hadroproduction cross sections of
charmonia Butenschoen:2010rq ; Butenschoen:2011yh ; Butenschoen:2012px ;
Butenschoen:2012qr ; Chao:2012iv ; Gong:2012ug ; Bodwin:2014gia  and
bottomonia Gong:2010bk ; Wang:2012is ; Gong:2013qka . However, this approach
has essentially failed all other phenomenological tests. It is incapable
Butenschoen:2010rq ; Butenschoen:2011yh ; Butenschoen:2012px ;
Butenschoen:2012qr  of describing the unpolarised photo- and hadroproduction
data together with the data on quarkonium polarisation observables using a
single set of Long-Distance Matrix Elements (LDMEs). This problem is known as
the “polarisation puzzle”. Moreover, the hadroproduction of $\eta_{c}$ at the
LHC has unexpectedly turned out to be dominated by the color-singlet state
$c\bar{c}[{}^{1}S_{0}^{[1]}]$, with no color-octet contribution needed to
describe data Butenschoen:2014dra . This seems inconsistent with Heavy-quark
spin symmetry relations between the LDMEs of $\eta_{c}$ and $J/\psi$. This
inconsistency is often referred to as the “heavy-quark-spin-symmetry puzzle”.

Alternative theoretical approaches, such as the (Improved) Colour-Evaporation
Model Barger:1979js ; Barger:1980mg ; Gavai:1994in ; Ma:2016exq  have their
own phenomenological problems, e.g. being incapable to describe the
$J/\psi$-pair hadroproduction Lansberg:2020rft , single $J/\psi$
hadroproduction at large $p_{T}$ Lansberg:2020rft  or $e^{+}e^{-}\to
J/\psi+c+\bar{c}$ cross section Lansberg:2019adr . Moreover, recent
theoretical developments, such as potential-NRQCD Brambilla:2021abf  or soft-
gluon factorisation Ma:2017xno ; Li:2019ncs ; Chen:2020yeg  stay within the
non-relativistic paradigm for the bound-state and try to either reduce number
of free parameters in NRQCD or more accurately treat the kinematic effects of
soft gluon emissions at hadronisation stage.

Given the situation described above, the logical way to proceed is to
scrutinise all possible sources of numerically large higher-order corrections
to the hard-scattering coefficient function in the NRQCD approach beyond
fixed-order NLO in $\alpha_{s}$ computations and find possible ways to resum
them to all orders in perturbation theory. For example, at high $p_{T}$, much
larger than the quarkonium mass $M$, such a resummation has already being
achieved via the fragmentation formalism which has been worked out up to the
next-to-leading power in $M/p_{T}$ Kang:2011mg ; Bodwin:2014gia . In the
present paper, we focus on a class of higher-order QCD corrections, which has
been largely ignored so far in the CF heavy-quarkonium-production literature
and which become important when the hadronic collision energy, $\sqrt{s}$, is
much larger than any other scale of our process, i.e. $M$ or $p_{T}$.

For simplicity, we will focus on the case of the inclusive hadroproduction of
the ${}^{1}S_{0}^{[1]}$, ${}^{3}P_{0}^{[1]}$ and ${}^{3}P_{2}^{[1]}$ NRQCD
Fock states of a heavy-quark pair of mass $M$. These states can be produced at
LO in CF via a $2\to 1$ process, namely by fusion of two on-shell gluons. The
NLO CF corrections to the total hadroproduction cross section at most involves
$2\to 2$ processes and can be computed in a closed form, which has been done
long time ago Kuhn:1992qw ; Petrelli:1997ge . In CF, the total hadroproduction
cross section of a state $m={}^{2S+1}L_{J}^{[1,8]}$ in a $h_{A}+h_{B}$
hadronic collision at a center-of-mass energy $\sqrt{s}$ can be computed as a
convolution:

$\sigma^{[m],h_{A}h_{B}}(\sqrt{s})=\sum\limits_{i,j=q,\bar{q},g}\hskip
2.84526pt\int\limits_{z_{\min}}^{1}\frac{dz}{z}{\cal
L}^{h_{A}h_{B}}_{ij}(z,\mu_{F})\hat{\sigma}^{[m]}_{ij}(z,\mu_{F},\mu_{R}),$
(1)

over the variable $z=M^{2}/\hat{s}$ which represents the squared fraction of
the energy of the partons initiating the hard subprocess used to produce the
observed final-state of mass $M$, with $z_{\min}=M^{2}/s$. The partonic
luminosity for partons $i,j=q,\bar{q},g$, entering eqn. (1) is:

${\cal L}^{h_{A}h_{B}}_{ij}(z,\mu_{F})=\int\limits_{-y_{\max}}^{+y_{\max}}dy\
\tilde{f}^{h_{A}}_{i}\left(\frac{M}{\sqrt{sz}}e^{y},\mu_{F}\right)\tilde{f}^{h_{B}}_{j}\left(\frac{M}{\sqrt{sz}}e^{-y},\mu_{F}\right),$
(2)

where $y_{\max}=\ln\sqrt{sz}/M$ and222In what follows, we will omit $h_{A}$
and $h_{B}$ when referring to PDFs.
$\tilde{f}^{h}_{i}(x,\mu_{F}^{2})=xf^{h}_{i}(x,\mu_{F}^{2})$ are the collinear
“momentum-density” PDFs of the parton $i$ in a hadron $h$, evaluated at the
factorisation scale $\mu_{F}$. The quantity $\hat{\sigma}^{[m]}_{ij}$ in eqn.
(1) is the hard-scattering coefficient for the corresponding partonic channel
$ij$ and state $m$.

The partonic luminosity is always decreasing when $z\to M^{2}/s$ corresponding
to $y_{\max}\to 0$. However, at $\mu_{F}\gg 1$ GeV this decrease is more rapid
than at $\mu_{F}\sim 1$ GeV, characteristic of quarkonium physics, due to
$\mu_{F}$-evolution of the $x$-dependence of PDFs in eqn. (2). It turns out
that, at high scales, the contribution of small values of $z$ to the integral
(1) is suppressed by the partonic luminosity. However, at smaller scales
$\mu_{F}\sim M\sim 1$ GeV most of the recent PDFs are relatively flat as
functions of $x$, and the region of small $z$ starts to significantly
contribute to the cross section at high energies $\sqrt{s}\gg M$, see the
detailed discussion in ref. Lansberg:2020ejc . For quarkonium-production cross
sections at NLO in CF, the effect of the small-$z$ region is dramatic, as it
was understood already in refs. Schuler:1994hy ; Mangano:1996kg .

For our forthcoming discussion, let us express the CF coefficient function
from eqn. (1) as follows:

$\hat{\sigma}^{[m]}_{ij}=\sigma_{0}^{[m]}\left[A_{0}^{[m]}\delta(1-z)+C_{ij}\frac{\alpha_{s}(\mu_{R})}{\pi}\left(A_{0}^{[m]}\ln\frac{M^{2}}{\mu_{F}^{2}}+A_{1}^{[m]}\right)+{\cal
O}(z\alpha_{s},\alpha_{s}^{2})\right],$ (3)

where the first term in square brackets contributes only for the states
$Q\bar{Q}[m]$ which can be produced in fusion of two on-shell gluons at ${\cal
O}(\alpha_{s}^{2})$. For such states, $A^{[m]}_{0}=1$. For all the other
$Q\bar{Q}[m]$ states, $A^{[m]}_{0}=0$ so that $\hat{\sigma}_{ij}^{[m]}$ starts
at ${\cal O}(\alpha_{s}^{3})$. The second term in the square brackets of the
eqn. (3) is the $z\to 0$ asymptotics of the ${\cal O}(\alpha_{s}^{3})$ part of
the coefficient function, while all the remaining terms are collected in the
term ${\cal O}(\alpha_{s}z,\alpha_{s}^{2})$. Overall factors
$\sigma_{0}^{[m]}$ are defined in eqn. (18) below.

The Fock-state dependent constants $A_{1}^{[m]}$ in eqn. (3) are collected,
for the cases when $A_{0}^{[m]}\neq 0$, in table 1 of ref. Lansberg:2020ejc
and they turn out to be negative. The colour-factors $C_{ij}$ in eqn. (3) are:
$C_{gg}=2C_{A}=2N_{c}$, $C_{qg}=C_{gq}=C_{F}=(N_{c}^{2}-1)/(2N_{c})$ and
$C_{q\bar{q}}=0$ (with $N_{c}=3$ being the number of colours). If the series
(3) is truncated at NLO in $\alpha_{s}$, then for scale choices satisfying
$\mu_{F}\geq M$, the contribution of the region $z\ll 1$ to the integral (1)
becomes purely negative and, at sufficiently high energies, it may outweigh
the positive NLO contribution from the $z\to 1$ limit and even the LO
contribution.

For charmonium production characterised by $M\simeq 3$ GeV, this already
happens at $\sqrt{s}$ as modest as several hundreds of GeV Mangano:1996kg ;
Lansberg:2020ejc  making NLO calculation completely unpredictive at LHC
energies. We find the problem described above very interesting because similar
instabilities plague the energy-dependence of rapidity-differential cross
sections of $J/\psi$ hadroproduction in NRQCD factorisation at NLO
Feng:2015cba  and could also affect the $p_{T}$-differential cross section at
$p_{T}\lesssim M\ll\sqrt{s}$.

As it was realised in ref. Lansberg:2020ejc , authored by two of us, these
negative cross sections in CF at NLO stem from an oversubtraction of the
collinear divergences inside the renormalised PDFs within the
$\overline{\text{MS}}$ scheme. In principle, such a subtraction should be
compensated by the evolution of the PDFs which progressively become steeper
when $\mu_{F}$ increases. However, the coefficients $A_{1}$, which are related
to the internal structure of the $gg\to Q\bar{Q}[m]$ form factor, are process-
dependent as opposed to the PDF evolution governed by the Dokshitzer-Gribov-
Lipatov-Altarelli-Parisi (DGLAP) DGLAP1 ; DGLAP2 ; DGLAP3  equations. As such,
the negative numbers cannot be systematically compensated by the PDF
evolution. This mismatch between the partonic cross section and the PDFs is
most dramatic for low scale processes for which the PDFs are the flattest.

The small-$z$ behaviour of the NLO partonic cross section of the type (3) is
characteristic of many hard processes, as it was shown for the first time in
ref. Ellis:1990hw  for the cases of total open heavy-flavour hadro- and
photoproduction as well as prompt-photon hadroproduction. In these cases, the
coefficients $A_{1}$ turn out to be positive, whereas they are negative for
the hadroproduction of ${}^{1}S_{0}$, ${}^{3}P_{0}$ and ${}^{3}P_{2}$-states
of heavy quark–anti-quark pairs as mentioned above, as well as for the
photoproduction of ${}^{3}S_{1}$ Serri:2021fhn . For Higgs-boson
hadroproduction, the sign of the coefficient $A_{1}$ depends on the ratio
$M_{H^{0}}/m_{t}$ (see ref. Lansberg:2020ejc ). These observations clearly
point at the fact that it is impossible to absorb the impact of $A_{1}$ into
the PDFs via the global scheme redefinition proposed in ref. Ellis:1990hw :
the optimal subtraction scheme for the open heavy-flavour hadroproduction will
be worse than $\overline{\text{MS}}$ for heavy quarkonia at high energies, and
vice-versa.

In ref. Lansberg:2020ejc , two of us have proposed a new scale prescription
for the factorisation scale to cure the mismatch between the partonic cross
section and the PDF evolution (for processes with $A_{0}=1$):

$\hat{\mu}_{F}=M\exp\left[{{A_{1}^{[m]}}/{2}}\right],$ (4)

which restores the positivity of quarkonium-production cross sections (see
refs. Lansberg:2020ejc ; Serri:2021fhn ). As it will be explained in section
2.5, this scale choice corresponds to an attempt to resum some higher-order
QCD corrections proportional to:

$\alpha_{s}^{n}\ln^{n-1}\frac{1}{z},$ (5)

at leading power in $z$ for $z\ll 1$ within $\hat{\sigma}_{ij}$. In what
follows, we will refer to the all order in $\alpha_{s}$ resummation of the
contributions of the type (5) as the Leading Logarithmic-$\ln(1/z)$
Approximation or LL($\ln(1/z)$) for short. The advantage of the
$\hat{\mu}_{F}$-prescription in comparison to the global scheme redefinition
is that it is process-dependent, i.e. every process is evaluated with its own
scale (4).

However, beyond ${\cal O}(\alpha_{s})$, the approximate resummation via the
$\hat{\mu}_{F}$-prescription does not correctly capture the structure of QCD
corrections of the type (5). For instance, as we will see, the
$\alpha_{s}^{2}\ln 1/z$ coefficient of the expanded expression of the NLO
cross section obtained with $\mu_{F}=\hat{\mu}_{F}$ does not correspond to the
one obtained from the resummation formalism. The systematic formalism for such
a resummation requires the use of the Balitsky-Fadin-Kuraev-Lipatov (BFKL)
evolution BFKL1 ; BFKL2 ; BFKL3  of the partonic amplitude in the rapidity
$Y\sim\ln(1/z)\gg 1$ with its non-trivial transverse-momentum dynamics to
correctly capture higher-order QCD corrections even in the LL($\ln(1/z)$)
approximation. Such a formalism is known as High Energy Factorisation (HEF)
and has been developed in refs. Catani:1990xk ; Catani:1990eg ; Collins:1991ty
; Catani:1994sq  in the LL($\ln(1/z)$) approximation. In the past, this
formalism has been successfully applied to study the high-energy structure of
hard-scattering coefficients of CF for many processes, such as heavy flavor
leptoproduction Catani:1992rn , Higgs Hautmann:2002tu ; Harlander:2009my ,
Drell-Yan lepton-pair Marzani:2008uh  and prompt-photon production
Diana:2009xv . In particular, in ref. Marzani:2008uh  the correctness of the
resummation predictions up to NNLO in $\alpha_{s}$ has been verified.
Recently, the NLL($\ln(1/z)$) resummation for lepton-hadron Deep-Inelastic
Scattering has been shown to significantly improve the quality of NNLO PDF
fits Ball:2017otu ; Abdolmaleki:2018jln .

With the present paper, we fill the gap in the existing HEF literature,
providing the first NLO CF + LL($\ln(1/z)$) matched calculation of the energy
dependence of the total heavy-quarkonium-hadroproduction cross section. The
aforementioned HEF calculations Catani:1992rn ; Hautmann:2002tu ;
Harlander:2009my ; Marzani:2008uh ; Diana:2009xv  use the representation of
the resummed cross section in Mellin space, which is convenient from the
computational point of view, but is difficult to match with NLO CF results.
For example, in the Ref. Catani:1992rn  the matching between LL HEF and NLO CF
calculations has been implemented via subtractive prescription in the Mellin
space ($N$) conjugate to $z$. With such a technique, however, it would be
quite difficult to come up with any alternative matching prescription, which
is necessary to estimate the uncertainty on the obtained cross section results
due to matching procedure. In contrast, we perform a direct $z$-space matching
between the HEF-resummed hard-scattering coefficient, which is valid at $z\ll
1$, and the full NLO CF correction, which contains numerically important
${\cal O}(z)$ power corrections. Along the same lines as in ref.
Echevarria:2018qyi , we propose a new matching procedure, described in section
3.2, which smoothly interpolates between $z\to 0$ and $z\to 1$ limits.

Let us also discuss now the relation of our work with other approaches to
quarkonium production at high energies existing in the small-$x$ physics
literature, such as calculations of the $J/\psi$-hadroproduction cross section
within Color-Glass-Condensate (CGC) formalism in refs. Ma:2014mri ;
Kang:2013hta . In these papers, the all-order in $\alpha_{s}$ resummation of
radiative corrections to the observable cross section $\sigma(\sqrt{s})$
enhanced by $\ln(s/M^{2})$ was the primary objective. The CGC formalism in
principle includes all such corrections in the LLA w.r.t. $\ln(s/M^{2})$ in
the limit $\sqrt{s}\gg M$, yet disregarding333It is interesting to note, that
the colour-singlet $Q\bar{Q}[{}^{3}S_{1}^{(1)}]$ state is produced in refs.
Ma:2014mri ; Kang:2013hta  via the subleading twist contribution, $g+(g+g)\to
Q\bar{Q}[{}^{3}S_{1}^{(1)}]$, requiring at least two off-shell gluons to be
taken from the “target” proton or nucleus. The leading-twist process $g+g\to
Q\bar{Q}[{}^{3}S_{1}^{(1)}]+g$ appears at NLO in $\alpha_{s}$ in the CGC
counting.  the “twist” expansion in powers of the hard scale $M^{2}$. In
contrast, in the present study, we are resumming large logarithms (5) in the
coefficient function of the CF $\hat{\sigma}(z)$, i.e. we are working in the
leading power approximation w.r.t. the hard scale $M^{2}$. As such, our
calculation includes at least partially the leading-twist part of the
logarithmic corrections resummed in Ma:2014mri ; Kang:2013hta . Indeed, at
relatively small scales $\mu_{F}\sim M\sim 1$ GeV, the partonic luminosity
function ${\cal L}_{ij}(z,\mu_{F})$ in eqn. (1) varies slowly for all values
of $z$ except those very close to the kinematic endpoint $z_{\min}=M^{2}/s$,
where it abruptly goes to zero, and therefore large logarithms (5) in
$\hat{\sigma}$ lead to the contributions proportional to:

$\int\limits_{z_{\min}}^{1}\frac{dz}{z}\
\alpha_{s}^{n}\ln^{n-1}\frac{1}{z}\sim\alpha_{s}^{n}\ln^{n}\frac{1}{z_{\min}}=\alpha_{s}^{n}\ln^{n}\frac{s}{M^{2}},$
(6)

which is another way to explain why the corrections of the type (5) are
especially important at $s\gg M^{2}$. The logarithmic contributions of the
type (6) can not be resummed completely in our calculation, because the HEF
formalism is based on the linear BFKL equation, while corrections of the type
(6) lead to the non-linear effects already in the LLA($\ln(s/M^{2})$). These
are taken into account, in the planar limit, by Balitsky-Kovchegov
Balitsky:1995ub ; Kovchegov:1999yj  or exactly in $N_{c}$, by
(Balitsky-)Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner
Balitsky:1998kc ; Jalilian-Marian:1997qno ; Jalilian-Marian:1997jhx ;
Jalilian-Marian:1997ubg ; Iancu:2001ad ; Iancu:2000hn  evolution equations.
However, we believe that for single-scale observables, such as the total cross
section $\sigma(\sqrt{s})$, in the leading twist approximation, all LLA
contributions of the type (6) can be included into the LLA HEF calculation
either via the perturbative resummation of logarithms (5) in $\hat{\sigma}$ or
phenomenologically, via the modification of the $x$ dependence of the
collinear PDFs. This point of view is supported by the fact that higher-twist
corrections originating from the modification of the ${\bf q}_{T}$-dependence
of the HEF resummation function in the region $|{\bf q}_{T}|\ll M$ are small,
see appendix C. Our conclusion is that the resummation of logarithms (5) as
well as its smooth matching with the NLO CF corrections at $z\sim 1$ are the
most numerically important effects for the single-scale observables such as
$\sigma(\sqrt{s})$ at high energy and in the leading-twist approximation.

The paper has the following structure. In section 2, we describe the
LL($\ln(1/z)$) HEF formalism, in particular the structure of the resummed
partonic cross section in $z$-space (section 2.1), the corresponding process-
dependent coefficient functions (section 2.2) and universal resummation
factors (section 2.3). In section 2.4, we expand the resummed cross section up
to NNLO in $\alpha_{s}$ to verify its consistency with NLO results Kuhn:1992qw
; Petrelli:1997ge  and provide predictions for future NNLO calculations of
quarkonium production cross section. In section 2.5, we explain the relation
of the HEF formalism with the $\hat{\mu}_{F}$ scale prescription. In section
3, we introduce two matching procedures in $z$-space, compare the
corresponding numerical results and discuss the uncertainties of our
calculation. Finally, in section 4, we summarise our conclusions and prospects
for future studies. The paper contains three appendices: in appendix A, we
explain two methods which we have developed to calculate the resummed cross
section in $z$-space in a numerically-stable way. In appendix B, we show the
numerical effects of the ${\cal O}(\alpha_{s}^{2})$ corrections to the
anomalous dimension (26) on the $\mu_{F}$ dependence of the resummed part of
the cross section and, in appendix C, we discuss the possible size of non-
perturbative effects, due to the intrinsic transverse momentum of gluons in
the proton, on the value of the total quarkonium production cross section in
our approximation.

##  2 High-energy factorisation in the leading-logarithmic approximation

###  2.1 The resummed cross section in $z$-space

Considering the $qg$ channel as an example, the corrections proportional to
$\alpha_{s}^{n}\ln^{n-1}1/z$ to the QCD hard subprocess for the production of
the final-state of interest, $Q\bar{Q}[m]$, with a four-momentum $p$ come from
processes with at most $n$ additional partons in the final state:

$g(p_{1})+q(p_{2})\to
g(k_{1})+g(k_{2})+\ldots+Q\bar{Q}[m](p)+\ldots+g(k_{n-1})+q(k_{n}),$ (7)

where the four-momentum labels are given in parentheses.

Figure 1: Typical diagrams contributing to the high-energy factorisation
amplitude in the $gq$-channel. The dashed lines denote Reggeised gluon
exchanges, while solid circles denote Lipatov’s vertices.

At leading-power in $z$ and in the LL($\ln(1/z)$) approximation, only the
final states with a strong ordering in the corresponding light-cone components
of momenta contribute, that is when:

$\displaystyle p_{1}^{+}$ $\displaystyle\simeq$ $\displaystyle k_{1}^{+}\gg
k_{2}^{+}\gg\ldots\gg p^{+}\gg\ldots\gg k_{n-1}^{+}\gg k_{n}^{+},$
$\displaystyle k_{1}^{-}\ll k_{2}^{-}\ll\ldots\ll p^{-}\ll\ldots\ll
k_{n-1}^{-}\ll k_{n}^{-}\simeq p_{2}^{-},$

where $k_{i}^{\pm}=k_{i}^{0}\pm k_{i}^{3}$, so that, for the four-momenta of
the incoming protons, one has $P_{1}^{+}=\sqrt{s}$, $P_{1}^{-}=0$ and
$P_{2}^{+}=0$, $P_{2}^{-}=\sqrt{s}$. This limit is referred to as the Multi-
Regge Kinematics (MRK). In this limit and at this accuracy, the QCD amplitudes
are known to factorise into a product of gauge-invariant blocks: the Lipatov’s
scattering, the production vertices BFKL1 ; BFKL2 ; BFKL3 , the $t$-channel
Reggeised gluon propagators and the effective vertex for production of the
$Q\bar{Q}[m]$ state, see the figure 1. For a review, the reader can consult
the monographs IFLQCD ; kovchegov_levin_2012  or the lecture notes
RevDelDuca95 . This factorisation at the amplitude level allows one to
factorise out the corresponding cross-section-level blocks: the universal
resummation factors ${\cal C}_{gi}$, with $i=g,q,\bar{q}$, and the process-
dependent coefficient function $H^{[m]}$.

The resummation factors are solutions of the BFKL equation, describing the
rapidity-evolution of the cascade of emissions from the rapidities of the most
backward/forward partons in the hard process to the rapidity of the observed
final state. The corresponding rapidity intervals (figure 1) can be estimated
in the MRK as:

$\displaystyle Y_{+}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\ln\frac{k_{1}^{+}}{k_{1}^{-}}-\ln\frac{p^{+}}{p^{-}}\right)\simeq\ln\left[\frac{p_{1}^{+}}{p_{+}}\frac{M_{T}}{|{\bf
k}_{T1}|}\right]=\ln\frac{1}{z_{+}}+\ln\frac{M_{T}}{|{\bf k}_{T1}|},$ (8)
$\displaystyle Y_{-}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\ln\frac{p^{+}}{p^{-}}-\ln\frac{k_{n}^{+}}{k_{n}^{-}}\right)\simeq\ln\left[\frac{p_{2}^{-}}{p_{-}}\frac{M_{T}}{|{\bf
k}_{Tn}|}\right]=\ln\frac{1}{z_{-}}+\ln\frac{M_{T}}{|{\bf k}_{Tn}|},$ (9)

where $M_{T}^{2}=M^{2}+{\bf p}_{T}^{2}$, $z_{+}=p_{+}/p_{1}^{+}$ and
$z_{-}=p_{-}/p_{2}^{-}$ are fractions of (+/-) light-cone components entering
the hard process, which are used up to produce the final state $Q\bar{Q}[m]$.
In terms of $z_{+}$ and $z_{-}$, the “global” $z$ variable introduced above is
equal to:

$z=\frac{M^{2}}{\hat{s}}=\frac{M^{2}}{p_{1}^{+}p_{2}^{-}}=\frac{M^{2}}{M_{T}^{2}}z_{+}z_{-}.$
(10)

The LL BFKL evolution of the resummation factors ${\cal C}_{gi}$ in rapidity
resums higher-order corrections to the cross section proportional to
$\alpha_{s}^{n}Y_{\pm}^{n-1}$. If one is only interested in the resummation of
the $\ln 1/z$-corrections in the LL($\ln(1/z)$) approximation, as we do in the
present paper, then one solves the BFKL equation for resummation factors,
neglecting transverse-momentum logarithms in eqns. (8) and (9), i.e. putting
$Y_{\pm}\simeq\ln 1/z_{\pm}$, as it has been done in the seminal papers
Catani:1990xk ; Catani:1990eg ; Collins:1991ty ; Catani:1994sq .

If the transverse-momentum logarithms in eqns. (8) and (9) are not neglected,
then in addition to $\ln 1/z$-corrections the HEF actually resums the LL
$\alpha_{s}^{n}\ln^{2n}(|{\bf p}_{T}|/M)$ and part of the NLL
$\alpha_{s}^{n}\ln^{n}(|{\bf p}_{T}|/M)$ corrections, which are traditionally
considered within the framework of Transverse-Momentum Dependent (TMD)
factorisation CollinsQCD . This overlap between TMD and High-Energy
factorisations has been explored in recent papers Nefedov:2021vvy ;
Hentschinski:2021lsh .

Postponing a more detailed discussion of the resummation factors until section
2.3, let us write down the factorisation formula for the rapidity-differential
cross section in HEF:

$\displaystyle\frac{d\sigma}{dy}$ $\displaystyle=$
$\displaystyle\sum\limits_{i,j=g,q,\bar{q}}\int\limits_{0}^{\infty}d{\bf
q}_{T1}^{2}d{\bf
q}_{T2}^{2}\int\limits_{x_{1}}^{1}\frac{dz_{+}}{z_{+}}\tilde{f}_{i}\left(\frac{x_{1}}{z_{+}},\mu_{F}\right){\cal
C}_{gi}(z_{+},{\bf q}_{T1}^{2},\mu_{F},\mu_{R})$ (11) $\displaystyle\times$
$\displaystyle\int\limits_{x_{2}}^{1}\frac{dz_{-}}{z_{-}}\tilde{f}_{j}\left(\frac{x_{2}}{z_{-}},\mu_{F}\right){\cal
C}_{gj}(z_{-},{\bf
q}_{T2}^{2},\mu_{F},\mu_{R})\int\limits_{0}^{2\pi}\frac{d\phi}{2}\frac{H^{[m]}({\bf
q}_{T1}^{2},{\bf q}_{T2}^{2},\phi)}{M_{T}^{4}},$

where $x_{1}=M_{T}e^{y}/\sqrt{s}$ and $x_{2}=M_{T}e^{-y}/\sqrt{s}$ with
$M_{T}^{2}=M^{2}+({\bf q}_{T1}+{\bf q}_{T2})^{2}$ and $\phi$ being the
azimuthal angle between ${\bf q}_{T1}$ and ${\bf q}_{T2}$. The HEF coefficient
function $H^{[m]}$ non-trivially depends on the transverse momenta ${\bf
q}_{T1,2}$ of the Reggeised gluons entering it (figure 1), and the procedure
for its computation is described in section 2.2. As it will become clear
below, this transverse-momentum dependence is the key to the factorisation of
higher-order terms enhanced by $\ln 1/z$ in the QCD perturbative series in
terms of the process-independent resummation functions ${\cal C}_{gj}$,
encapsulating logarithms $\ln 1/z_{\pm}$ and the process dependent
coefficients $H^{[m]}$, which are free from these large logarithms.

Integrating eqn. (11) over the rapidity and introducing the variables
$\eta=\ln(z_{+}/z_{-})/2$ and $z$ via eqn. (10), one can recast the total
hadronic cross section in HEF into the form of eqn. (1) with the partonic
coefficient:

$\displaystyle\hat{\sigma}_{ij}^{[m],{\rm\ HEF}}(z,\mu_{F},\mu_{R})$
$\displaystyle=$
$\displaystyle\int\limits_{-\infty}^{\infty}d\eta\int\limits_{0}^{\infty}d{\bf
q}_{T1}^{2}d{\bf q}_{T2}^{2}\ {\cal
C}_{gi}\left(\frac{M_{T}}{M}\sqrt{z}e^{\eta},{\bf
q}_{T1}^{2},\mu_{F},\mu_{R}\right)$ (12) $\displaystyle\times$
$\displaystyle{\cal C}_{gj}\left(\frac{M_{T}}{M}\sqrt{z}e^{-\eta},{\bf
q}_{T2}^{2},\mu_{F},\mu_{R}\right)\int\limits_{0}^{2\pi}\frac{d\phi}{2}\frac{H^{[m]}({\bf
q}_{T1}^{2},{\bf q}_{T2}^{2},\phi)}{M_{T}^{4}},$

where the integration over $\eta$ actually proceeds only over
$|\eta|\leq\eta_{\rm max}=\ln\left[M/(M_{T}\sqrt{z})\right]$, because the
arguments of the resummation functions $z_{\pm}$ should be smaller or equal
than one.

Finally let us emphasise that, for the $gg$ channel, the transverse-momentum
integrals in eqn. (12) can not be dealt with numerically unless the LO
contribution $\sigma_{0}^{[m]}\delta(1-z)$ is explicitly removed from eqn.
(12):

$\hat{\sigma}_{gg}^{[m],\ {\rm
HEF}}(z,\mu_{F},\mu_{R})=\sigma_{0}^{[m]}\delta(1-z)+\check{\sigma}_{gg}^{[m],\
{\rm HEF}}(z,\mu_{F},\mu_{R}),$ (13)

where $\check{\sigma}_{gg}^{[m],\ {\rm HEF}}$ is the LO-subtracted resummed
hard-scattering coefficient. The LO CF contribution is contained in the
$\hat{\sigma}_{gg}^{[m],\ {\rm HEF}}$ because the perturbative expansion for
resummation factors ${\cal C}_{gg}$ starts with the LO term
$\delta(1-z)\delta({\bf q}_{T}^{2})$, which corresponds to the case with no
additional emissions. The procedure (13) turns the integrals in eqn. (12) into
well-defined functions of $z$. We describe two methods which we have used to
isolate the LO contribution in the Appendix A. The method of the section A.1
relies on the form of resummation factor in the Doubly-Logarithmic
approximation, which we will introduce in eqn. (36) in section 2.3, while the
method in section A.2 is more general and is applicable to the exact
LL($\ln(1/z)$) resummation factor (24) and its generalisations.

###  2.2 HEF coefficient functions

For the following $2\to 1$ processes:

$R(q_{1})+R(q_{2})\to
c\bar{c}\left[{}^{1}S_{0}^{[1,8]},{}^{3}S_{1}^{[8]},{}^{3}P_{0,1,2}^{[1,8]}\right],$
(14)

the corresponding HEF coefficient functions are non-zero at LO in
$\alpha_{s}$. In eqn. (14), the incoming Reggeised gluons are denoted by the
symbol $R$ and their four-momenta can be parametrised as
$q_{1}^{\mu}=x_{1}P_{1}^{\mu}+q_{T1}^{\mu}$ and
$q_{2}^{\mu}=x_{2}P_{2}^{\mu}+q_{T2}^{\mu}$ (figure 1) so that
$q_{1,2}^{2}=-{\bf q}_{T1,2}^{2}<0$. The Feynman rules of Lipatov’s gauge-
invariant EFT for Multi-Regge processes in QCD Lipatov95  are used for the
computation of HEF coefficient functions for general QCD processes, as it was
done for quarkonia in Ref. Kniehl:2006vm . The use of the High-Energy EFT is
particularly important for the gauge invariance444With respect to gauge choice
of internal gluon propagator. of the computation of the coefficient function
for production of the ${}^{3}S_{1}^{[8]}$-state. Essentially, the Lipatov’s
vertex should be used instead of the usual three-gluon vertex.555Indeed, for
the states listed in eqn. (14), except for ${}^{3}S_{1}^{[8]}$, the EFT
computation is equivalent to the use of the following “nonsense” polarisation
prescription for the initial-state gluons:
$\varepsilon^{\mu}(q_{1,2})\to\frac{q_{T1,2}^{\mu}}{|{\bf q}_{T1,2}|},$ (15)
due to Slavnov-Taylor identities of QCD. However, for more complicated
processes containing final-state gluons coupling to initial-state ones, or
otherwise essentially non-Abelian, as e.g. double color-octet channels in
heavy-quarkonium pair production He:di-Jpsi , the full set of Feynman rules of
the EFT Lipatov95  have to be used instead of the prescription (15). The
coefficient functions for the processes (14) have been computed for the first
time in Hagler:2000dda ; Hagler:2000dd ; Kniehl:2006sk  using the prescription
(15) and later have been reproduced in the ref. Kniehl:2006vm  using Lipatov’s
EFT.

In the present paper we take the explicit expressions for HEF coefficient
functions, which can be directly plugged into our eqn. (12), from ref.
Kniehl:2006sk . The HEF coefficient functions satisfy the on-shell-limit
normalisation condition:

$\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}\lim\limits_{{\bf q}_{T1,2}^{2}\to
0}H^{[m]}({\bf q}^{2}_{T1},{\bf
q}_{T2}^{2},\phi)=\frac{1}{4(N_{c}^{2}-1)^{2}}\sum\limits_{\lambda_{1,2}=\pm}\left|{\cal
M}(g_{\lambda_{1}}g_{\lambda_{2}}\to Q\bar{Q}[m])\right|^{2},$ (16)

with the corresponding CF squared amplitude ${\cal M}$ averaged over colours
and helicities $\lambda_{1,2}$ of the initial-state gluons. It is convenient
to separate out the corresponding colour factors, the LDMEs of NRQCD
$\langle{\cal O}[m]\rangle$ and the CF spin-orbital factors from the
coefficient functions as follows:

$H^{[m]}({\bf q}_{T1}^{2},{\bf
q}_{T2}^{2},\phi)=\frac{M^{4}\sigma_{0}^{[m]}}{\pi}F^{[m]}({\bf
q}_{T1}^{2},{\bf q}_{T2}^{2},\phi),$ (17)

where the $\sigma_{0}^{[m]}$ factors are:

$\displaystyle\sigma_{0}^{[{}^{1}S_{0}^{[1]}]}$ $\displaystyle=$
$\displaystyle\frac{8}{C_{A}^{3}C_{F}}\pi^{3}\alpha_{s}^{2}\frac{\langle{\cal
O}[{}^{1}S_{0}^{[1])}]\rangle}{M^{5}},$ (18)
$\displaystyle\sigma_{0}^{[{}^{1}S_{0}^{[8]}]}$ $\displaystyle=$
$\displaystyle\frac{4(C_{A}^{2}-4)}{C_{A}^{3}C_{F}^{2}}\pi^{3}\alpha_{s}^{2}\frac{\langle{\cal
O}[^{1}S_{0}^{[8]}]\rangle}{M^{5}},$
$\displaystyle\sigma_{0}^{[{}^{3}S_{1}^{[8]}]}$ $\displaystyle=$
$\displaystyle\frac{8(C_{A}-2C_{F})}{3C_{F}^{2}}\pi^{3}\alpha_{s}^{2}\frac{\langle{\cal
O}[^{3}S_{1}^{[8]}]\rangle}{M^{5}},$
$\displaystyle\sigma_{0}^{[^{3}P_{0}^{[1]}]}$ $\displaystyle=$
$\displaystyle\frac{96}{C_{A}^{3}C_{F}}\pi^{3}\alpha_{s}^{2}\frac{\langle{\cal
O}[^{3}P_{0}^{[1]}]\rangle}{M^{7}},$
$\displaystyle\sigma_{0}^{[^{3}P_{1}^{[1]}]}$ $\displaystyle=$
$\displaystyle\frac{192}{C_{A}^{3}C_{F}}\pi^{3}\alpha_{s}^{2}\frac{\langle{\cal
O}[^{3}P_{1}^{[1]}]\rangle}{M^{7}},$
$\displaystyle\sigma_{0}^{[^{3}P_{2}^{[1]}]}$ $\displaystyle=$
$\displaystyle\frac{128}{5C_{A}^{3}C_{F}}\pi^{3}\alpha_{s}^{2}\frac{\langle{\cal
O}[^{3}P_{2}^{[1]}]\rangle}{M^{7}},$
$\displaystyle\sigma_{0}^{[^{3}P_{0}^{[8]}]}$ $\displaystyle=$
$\displaystyle\frac{48(C_{A}^{2}-4)}{C_{A}^{3}C_{F}^{2}}\pi^{3}\alpha_{s}^{2}\frac{\langle{\cal
O}[^{3}P_{0}^{[8]}]\rangle}{M^{7}},$
$\displaystyle\sigma_{0}^{[^{3}P_{1}^{[8]}]}$ $\displaystyle=$
$\displaystyle\frac{96(C_{A}^{2}-4)}{C_{A}^{3}C_{F}^{2}}\pi^{3}\alpha_{s}^{2}\frac{\langle{\cal
O}[^{3}P_{1}^{[8]}]\rangle}{M^{7}},$
$\displaystyle\sigma_{0}^{[^{3}P_{2}^{[8]}]}$ $\displaystyle=$
$\displaystyle\frac{64(C_{A}^{2}-4)}{5C_{A}^{3}C_{F}^{2}}\pi^{3}\alpha_{s}^{2}\frac{\langle{\cal
O}[^{3}P_{2}^{[8]}]\rangle}{M^{7}},$

where $C_{A}=N_{c}$, $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ and the averaging factors
accounting for the number of polarisation ($2J+1$) and colour-singlet
($2N_{c}$) or colour-octet ($N_{c}^{2}-1$) states of the $Q\bar{Q}[m]$-pair
have been included into the denominators of eqns. (18) according to the
standard definition of LDMEs. The functions $F^{[m]}$ in the eqn. (17) do not
depend on the colour state of the $Q\bar{Q}$ pair and depend only on its spin
and orbital momentum. For the $S$-states they are:

$\displaystyle F^{[^{1}S_{0}]}(t_{1},t_{2},\phi)$ $\displaystyle=$
$\displaystyle 2\frac{\left(M^{2}+{\bf
p}_{T}^{2}\right)^{2}}{(M^{2}+t_{1}+t_{2})^{2}}\sin^{2}{\phi},$ (19)
$\displaystyle F^{[^{3}S_{1}]}(t_{1},t_{2},\phi)$ $\displaystyle=$
$\displaystyle\frac{\left(M^{2}+{\bf
p}_{T}^{2}\right)\left[(t_{1}+t_{2})^{2}+M^{2}\left(t_{1}+t_{2}-2\sqrt{t_{1}t_{2}}\cos\phi\right)\right]}{M^{2}(M^{2}+t_{1}+t_{2})^{2}},$
(20)

with ${\bf p}_{T}^{2}=({\bf q}_{T1}+{\bf
q}_{T2})^{2}=t_{1}+t_{2}+2\sqrt{t_{1}t_{2}}\cos\phi$ with $t_{1,2}={\bf
q}_{T1,2}^{2}$. For the $P$-wave states, one has:

$\displaystyle F^{[^{3}P_{0}]}(t_{1},t_{2},\phi)$ $\displaystyle=$
$\displaystyle\frac{2}{9}\,\frac{\left(M^{2}+{\bf
p}_{T}^{2}\right)^{2}\left[(3M^{2}+t_{1}+t_{2})\cos{\phi}+2\sqrt{t_{1}t_{2}}\right]^{2}}{(M^{2}+t_{1}+t_{2})^{4}},$
$\displaystyle F^{[^{3}P_{1}]}(t_{1},t_{2},\phi)$ $\displaystyle=$
$\displaystyle\frac{2}{9}\,\frac{\left(M^{2}+{\bf
p}_{T}^{2}\right)^{2}\left[(t_{1}+t_{2})^{2}\sin^{2}{\phi}+M^{2}\left(t_{1}+t_{2}-2\sqrt{t_{1}t_{2}}\cos{\phi}\right)\right]}{(M^{2}+t_{1}+t_{2})^{4}},$
$\displaystyle F^{[^{3}P_{2}]}(t_{1},t_{2},\phi)$ $\displaystyle=$
$\displaystyle\frac{1}{3}\,\frac{\left(M^{2}+{\bf
p}_{T}^{2}\right)^{2}}{(M^{2}+t_{1}+t_{2})^{4}}\left\\{3M^{4}+3M^{2}(t_{1}+t_{2})+4t_{1}t_{2}\right.$
(21)
$\displaystyle{}+\left.(t_{1}+t_{2})^{2}\cos^{2}{\phi}+2\sqrt{t_{1}t_{2}}\left[3M^{2}+2(t_{1}+t_{2})\right]\cos{\phi}\right\\}.$

The on-shell limits (16) of the functions $F^{[m]}$ are equal to one or zero,
depending on whether the corresponding state is allowed to be produced in a
fusion of two on-shell gluons by the Landau-Yang theorem, i.e. they are equal
to constants $A_{0}$ first introduced in eqn. (3).

###  2.3 Resummation factor in the doubly-logarithmic approximation and beyond

Following the original refs. Catani:1990xk ; Catani:1990eg ; Collins:1991ty ;
Catani:1994sq , we introduce the Mellin transform with respect to the
$z$-dependence of the resummation factor as follows:

${\cal C}(N,{\bf q}_{T},\mu_{F},\mu_{R})=\int\limits_{0}^{1}dz\ z^{N-1}{\cal
C}(z,{\bf q}_{T},\mu_{F},\mu_{R}),$ (22)

so that the logarithms which we aim to resum are mapped to poles at $N=0$ in
Mellin space:

$\ln^{k-1}\frac{1}{z}\leftrightarrow\frac{(k-1)!}{N^{k}}.$ (23)

The emissions of additional partons at higher orders in $\alpha_{s}$ generate
collinear divergences in the resummation factor which have been removed in
ref. Collins:1991ty  using a transverse-momentum cut off and using dimensional
regularisation and the $\overline{\text{MS}}$-subtraction prescription in ref.
Catani:1994sq . After the subtraction of the collinear divergences, the
dependence on factorisation scale $\mu_{F}$ arises in the resummation factor.
The Mellin-space result Catani:1990xk ; Catani:1990eg ; Collins:1991ty ;
Catani:1994sq  for the subtracted resummation factor is:

${\cal C}_{gg}(N,{\bf
q}_{T}^{2},\mu_{F},\mu_{R})=R(\gamma_{gg}(N,\alpha_{s}))\frac{\gamma_{gg}(N,\alpha_{s})}{{\bf
q}_{T}^{2}}\left(\frac{{\bf
q}_{T}^{2}}{\mu_{F}^{2}}\right)^{\gamma_{gg}(N,\alpha_{s})},$ (24)

where $\hat{\alpha}_{s}=\alpha_{s}(\mu_{R})C_{A}/\pi$. The choice of the scale
at which the $\alpha_{s}$ is evaluated in the ${\cal C}_{gg}$ function is not
strictly dictated by the LL($\ln(1/z)$) resummation or NLO CF results.
However, since both ${\cal C}_{gg}$ functions and HEF coefficient function
belong to the coeffcient function of CF $\hat{\sigma}_{ij}^{[m]}$ (see eqn.
(12)), we opt to choose the same scale $\mu_{R}$ of the $\alpha_{s}$ in all
these quantities. The anomalous dimension $\gamma_{gg}$ is the solution of the
algebraic equation first derived in the ref. Jaroszewicz:1982gr :

$\frac{\hat{\alpha}_{s}}{N}\chi(\gamma_{gg}(N,\alpha_{s}))=1,$ (25)

where $\chi(\gamma)=2\psi_{0}(1)-\psi_{0}(\gamma)-\psi_{0}(1-\gamma)$ is he
Lipatov’s LO characteristic function and
$\psi_{n}(\gamma)=d^{n}\ln\Gamma(\gamma)/d\gamma^{n}$ is the Euler’s $\psi$
function. The first few terms of the perturbative solution of eqn. (25) are:

$\gamma_{gg}(N,\alpha_{s})=\gamma_{N}+2\zeta(3)\frac{\hat{\alpha}_{s}^{4}(\mu_{R})}{N^{4}}+2\zeta(5)\frac{\hat{\alpha}_{s}^{6}(\mu_{R})}{N^{6}}+{\cal
O}(\alpha_{s}^{7}),{\rm\ where\
}\gamma_{N}=\frac{\hat{\alpha}_{s}(\mu_{R})}{N}.$ (26)

Eqn. (24) and the anomalous dimension (26) resum the series of higher-order
corrections proportional to $\alpha_{s}^{n}/N^{n}$. These corrections by
virtue of the mapping (23) are equivalent to the LL($\ln(1/z)$)-approximation
(5). While the LL($\ln(1/z)$) anomalous dimension (26) is scheme-independent
in a wide class of $\overline{\text{MS}}$-like schemes Catani:1994sq , the
factor $R(\gamma)$ encapsulates the subtraction-scheme dependence. For the
$\overline{\text{MS}}$ scheme, its perturbative expansion starts at N3LO
Catani:1994sq ; Kirschner:2009qu :

$R(\gamma_{gg}(N,\alpha_{s}))=1+{\cal O}(\alpha_{s}^{3}).$ (27)

Let us now define the resummation factors for the quark-induced channels. As
depicted in figure 1, in the LL($\ln(1/z)$)-approximation, only Reggeised
gluons can propagate in the $t$-channel. Hence, the quarks can only
participate in so-called LO partonic impact factors. The corresponding
expression for the squared amplitude in the LL($\ln(1/z)$) approximation
differs from the similar squared amplitude with gluon substituted by the quark
only by an overall colour factor. This picture leads to the simple relation
between the quark and gluon-induced resummation factors in the LL($\ln(1/z)$)
approximation:

${\cal C}_{gq}(z,{\bf
q}_{T}^{2},\mu_{F},\mu_{R})=\frac{C_{F}}{C_{A}}\left[{\cal C}_{gg}(z,{\bf
q}_{T}^{2},\mu_{F},\mu_{R})-\delta(1-z)\delta({\bf q}_{T}^{2})\right],$ (28)

where the subtraction term in square brackets removes the contribution of the
direct $g\to R$ transition in the on-shell limit. The latter is absent in case
of quarks due to the fermion-number conservation.

To clarify the physical meaning of the anomalous dimension $\gamma_{gg}$ in
eqn. (24), let us observe that, since the Mellin transform turns
$z_{\pm}-$convolutions in eqn. (11) into products, the $\mu_{F}$ independence
of the cross section order by order in $\alpha_{s}$ is achieved provided that:

$\sum\limits_{i=q,\bar{q},g}\frac{\partial}{\partial\ln\mu_{F}^{2}}\left[\tilde{f}_{i}(N,\mu_{F}){\cal
C}_{gi}(N,{\bf q}_{T}^{2},\mu_{F},\mu_{R})\right]=0,$ (29)

at any ${\bf q}_{T}$. Substituting eqns. (24) and (28) one finds that eqn.
(29) is indeed satisfied due to the following form of DGLAP equations at
LL($\ln 1/z$) accuracy:

$\displaystyle\frac{\partial\tilde{f}_{g}(N,\mu_{F})}{\partial\ln\mu_{F}^{2}}$
$\displaystyle=$
$\displaystyle\gamma_{gg}(N,\alpha_{s})\left[\tilde{f}_{g}(N,\mu_{F})+\frac{C_{F}}{C_{A}}\sum\limits_{i=q,\bar{q}}\tilde{f}_{i}(N,\mu_{F})\right],$
(30)
$\displaystyle\frac{\partial\tilde{f}_{q}(N,\mu_{F})}{\partial\ln\mu_{F}^{2}}$
$\displaystyle=$ $\displaystyle 0.$ (31)

In other words, the $\gamma_{gg}(N,\alpha_{s})$ is the LL($\ln(1/z)$)
approximation to the anomalous dimension of the DGLAP equations, governing the
scale dependence of the gluon momentum density distribution, while
$C_{F}\gamma_{gg}/C_{A}$ determines the feed down from quarks to gluon in this
approximation. Indeed, the corresponding LO DGLAP splitting functions have the
following $z\to 0$ asymptotics:

$\displaystyle\frac{\alpha_{s}}{2\pi}zP_{gg}(z)$ $\displaystyle=$
$\displaystyle\hat{\alpha}_{s}+{\cal O}(z),$ (32)
$\displaystyle\frac{\alpha_{s}}{2\pi}zP_{gq}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{F}}{C_{A}}\hat{\alpha}_{s}+{\cal O}(z).$ (33)

Eqn. (32) is equivalent to the first term of the solution (26) via the mapping
(23), while eqn. (33) is related in the same way to the first term of the
expansion of $C_{F}\gamma_{gg}/C_{A}$ in $\alpha_{s}$. The coefficients in
front of $\alpha_{s}\ln(1/z)$ and $\alpha_{s}^{2}\ln^{2}(1/z)$ terms in the
NLO and NNLO contributions to $zP_{gg}(z)$ are zero in QCD, consistently with
eqn. (26). The first non-zero contributions at leading power in $z$ beyond LO
in $\alpha_{s}$ in QCD are ${\cal O}(\alpha_{s}^{2})$ and ${\cal
O}(\alpha_{s}^{3}\ln(1/z))$. They belong to the NLL($\ln(1/z)$) approximation
and their consistency with NLL($\ln(1/z)$) BFKL predictions has been verified
in ref. Vogt:2004mw , while the consistency of the LL($\ln(1/z)$) series (26)
with the DGLAP anomalous dimension in ${\cal N}=4$-supersymmetric Yang-Mills
theory in the large-$N_{c}$ limit has been checked in ref. Lukowski:2009ce  up
to five loops using integrability-based methods.

The connection with the DGLAP anomalous dimension, outlined above, is very
important for the phenomenological strategy of the usage of the LL($\ln(1/z)$)
resummation. Most of the common PDF sets666With the notable exception of the
sets obtained in refs. Ball:2017otu ; Abdolmaleki:2018jln . use the fixed-
order approximation for the DGLAP splittings/anomalous dimensions, i.e. no
$\ln(1/z)$-resummation is performed in the evolution. In particular, it means
that existing LO, NLO and NNLO PDFs contain information only about the
$\gamma_{N}$ term in the series (26). Including further terms of this series
into the CF coefficient function via the resummation will only increase the
mismatch between the $\mu_{F}$ dependence of PDFs and the resummation factors
and, hence, blow up our $\mu_{F}$-scale uncertainty at high energy.

The goal of the present paper is to use HEF to cure the unphysical behaviour
of the CF coefficient function at small $z$ in a way consistent with the NLO
DGLAP evolution of PDFs. Therefore, we will use the following Doubly-
Logarithmic approximation (DLA) for resummation factor, which is obtained by
neglecting all terms of the series (26) except the first one:

${\cal C}_{gg}^{\rm DLA}(N,{\bf q}_{T},\mu_{F},\mu_{R})=\frac{\gamma_{N}}{{\bf
q}_{T}^{2}}\left(\frac{{\bf q}_{T}^{2}}{\mu_{F}^{2}}\right)^{\gamma_{N}}.$
(34)

In the DLA, only terms $\propto\left(\frac{\alpha_{s}}{N}\ln\frac{{\bf
q}_{T}^{2}}{\mu_{F}^{2}}\right)^{n}\leftrightarrow\alpha_{s}^{n}\ln^{n-1}\frac{1}{z}\ln^{n}\frac{{\bf
q}_{T}^{2}}{\mu_{F}^{2}}$ are included in the resummation factor. One also
takes $R(\gamma)=1$ in the eqn. (34) for the same reason: the corresponding
scheme-dependence starts at N3LO, eqn. (27), and is not taken into account in
the NLO or even NNLO PDFs. The inclusion of this factor to the resummation
function ${\cal C}$ would only lead to an unphysical subtraction-scheme
mismatch between the PDFs and the CF coefficient function.

The inverse Mellin transform:

${\cal C}(z,{\bf
q}_{T}^{2},\mu_{F},\mu_{R})=\int\limits_{-i\infty}^{+i\infty}\frac{dN}{2\pi
i}z^{-N}{\cal C}(N,{\bf q}_{T}^{2},\mu_{F},\mu_{R}),$ (35)

of the resummation factor in the DLA can be computed straightforwardly, e.g.
order-by-order in $\alpha_{s}$ with the use of relation (23), leading to the
following result:

${\cal C}_{gg}^{\text{DLA}}(z,{\bf
q}_{T}^{2},\mu_{F},\mu_{R})=\frac{\hat{\alpha}_{s}(\mu_{R})}{{\bf
q}_{T}^{2}}\left\\{\begin{matrix}J_{0}\left(2\sqrt{\hat{\alpha}_{s}(\mu_{R})\ln\left(\frac{1}{z}\right)\ln\left(\frac{\mu_{F}^{2}}{{\bf
q}_{T}^{2}}\right)}\right)&\text{if }{\bf q}_{T}^{2}<\mu_{F}^{2},\\\
I_{0}\left(2\sqrt{\hat{\alpha}_{s}(\mu_{R})\ln\left(\frac{1}{z}\right)\ln\left(\frac{{\bf
q}_{T}^{2}}{\mu_{F}^{2}}\right)}\right)&\text{if }{\bf
q}_{T}^{2}>\mu_{F}^{2},\end{matrix}\right.$ (36)

where $J_{0}/I_{0}$ are Bessel functions of first/second kind. Eqn. (36) is
known in the HEF community777See e.g. Ref. Kniehl:2011hc  for its application
to the inclusive jet production at the LHC. as the (Collins-Ellis-)Blümlein
formula, and first appears in ref. Blumlein:1995eu  where the solution of the
evolution equation written in ref. Collins:1991ty  has been studied.

The DLA resummation function (34), or equivalently (36), has a very important
normalisation property:

$\displaystyle\int\limits_{0}^{\mu_{F}^{2}}d{\bf q}_{T}^{2}\ {\cal
C}_{gg}^{\rm DLA}(N,{\bf
q}_{T}^{2},\mu_{F},\mu_{R})=1\leftrightarrow\int\limits_{0}^{\mu_{F}^{2}}d{\bf
q}_{T}^{2}\ {\cal C}_{gg}^{\rm DLA}(z,{\bf
q}_{T}^{2},\mu_{F},\mu_{R})=\delta(1-z),$ (37)

which will be used extensively in the rest of the paper. The similar
normalisation condition for ${\cal C}_{gq}$ in the DLA follows form eqn. (28):

$\int\limits_{0}^{\mu_{F}^{2}}d{\bf q}_{T}^{2}\ {\cal C}^{\rm DLA}_{gq}(z,{\bf
q}_{T}^{2},\mu_{F},\mu_{R})=0.$ (38)

###  2.4 Small-$z$ hard-scattering coefficient at NLO and NNLO from
resummation

As a consistency check, the LL($\ln(1/z)$) resummation described above should
reproduce the small-$z$ behaviour (3) of the known NLO result for the
coefficient function Kuhn:1992qw ; Petrelli:1997ge ; Lansberg:2020ejc . To
this end, one has to expand the resummed CF coefficient functions (12) at
least up to NLO in $\alpha_{s}$. In this section, we perform such an expansion
up to NNLO in order to provide predictions for the $\alpha_{s}^{2}\ln(1/z)$
NNLO terms in the CF coefficient functions. The results of this section will
also be employed in sec. 3.2 to construct the weight functions needed for our
matching procedure. The DLA resummation factor (34) is appropriate for the
expansion up to NNLO because the exact LL($\ln(1/z)$)-resummation factor (24)
differs from it only starting at ${\cal O}(\alpha_{s}^{3})$.

To perform such an expansion, it is convenient to work in Mellin space,
plugging in the inverse Mellin transform (35) of eqns. (34) into eqn. (12).
Passing to the dimensionless transverse momenta ${\bf n}_{T1}={\bf q}_{T1}/M$
and ${\bf n}_{T2}={\bf q}_{T2}/M$, one can rewrite the resummed CF coefficient
functions as:

$\hat{\sigma}^{[m],\text{\
HEF}}_{ij}=\sigma_{0}^{[m]}\int\limits_{-\infty}^{\infty}d\eta\int\limits_{-i\infty}^{+i\infty}\frac{dN_{+}dN_{-}}{(2\pi
i)^{2}}\
z^{-\frac{N_{+}+N_{-}}{2}}e^{\eta(N_{+}-N_{-})}G^{[m]}_{ij}\left(\gamma_{N_{+}},\gamma_{N_{-}},\frac{N_{+}+N_{-}}{2},\frac{M^{2}}{\mu_{F}^{2}}\right),$
(39)

where

$G^{[m]}_{gg}\left(\gamma_{1},\gamma_{2},\nu,\frac{M^{2}}{\mu_{F}^{2}}\right)=\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}+\gamma_{2}}\int\limits_{0}^{\infty}\frac{\gamma_{1}\gamma_{2}d{\bf
n}_{T1}^{2}d{\bf n}_{T2}^{2}}{({\bf n}_{T1}^{2})^{1-\gamma_{1}}({\bf
n}_{T2}^{2})^{1-\gamma_{2}}}\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}f^{[m]}({\bf
n}_{T1}^{2},{\bf n}_{T2}^{2},\phi,\nu),$ (40)

for the gluon channel, while for quark-induced channels, one takes into
account eqn. (28) to obtain

$\displaystyle G^{[m]}_{qg}$ $\displaystyle=$
$\displaystyle\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{2}}\frac{C_{F}}{C_{A}}\int\limits_{0}^{\infty}\frac{\gamma_{2}d{\bf
n}_{T2}^{2}}{({\bf n}_{T2}^{2})^{1-\gamma_{2}}}\int\limits_{0}^{\infty}d{\bf
n}_{T1}^{2}\left[\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}}\gamma_{1}({\bf
n}_{T1}^{2})^{-1+\gamma_{1}}-\delta({\bf n}_{T1}^{2})\right]$ (41)
$\displaystyle\times$
$\displaystyle\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}f^{[m]}({\bf
n}_{T1}^{2},{\bf n}_{T2}^{2},\phi),$

and

$\displaystyle G^{[m]}_{q\bar{q}}$ $\displaystyle=$
$\displaystyle\left(\frac{C_{F}}{C_{A}}\right)^{2}\int\limits_{0}^{\infty}d{\bf
n}_{T1}^{2}d{\bf
n}_{T2}^{2}\left[\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}}\gamma_{1}({\bf
n}_{T1}^{2})^{-1+\gamma_{1}}-\delta({\bf n}_{T1}^{2})\right]$ (42)
$\displaystyle\times$
$\displaystyle\left[\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{2}}\gamma_{2}({\bf
n}_{T2}^{2})^{-1+\gamma_{2}}-\delta({\bf
n}_{T2}^{2})\right]\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}f^{[m]}({\bf
n}_{T1}^{2},{\bf n}_{T2}^{2},\phi),$

where, in eqns. (40), (41) and (42), the dimensionless function

$f^{[m]}({\bf n}_{T1}^{2},{\bf n}_{T2}^{2},\phi,\nu)=\frac{\pi
H^{[m]}(M^{2}{\bf n}_{T1}^{2},M^{2}{\bf
n}_{T2}^{2},\phi)}{\sigma_{0}^{[m]}M^{4}(1+({\bf n}_{T1}+{\bf
n}_{T2})^{2})^{2+\nu}}=\frac{F^{[m]}(M^{2}{\bf n}_{T1}^{2},M^{2}{\bf
n}_{T2}^{2},\phi)}{(1+({\bf n}_{T1}+{\bf n}_{T2})^{2})^{2+\nu}},$ (43)

encapsulates the HEF coefficient function $H^{[m]}$ (or $F^{[m]}$ of eqn.
(17)) and the kinematic factor $1/M_{T}^{4}$ from eqn. (12).

The only quantity depending on $\alpha_{s}$ in eqn. (39) is
$\gamma_{N}=\hat{\alpha}_{s}/N$. One thus only has to Taylor-expand the
function $G_{ij}^{[m]}$ with respect to its arguments $\gamma_{N_{+}}$ and
$\gamma_{N_{-}}$. The corresponding poles in $N_{+}$ and $N_{-}$ map to the
following functions of $z$ via eqn. (39):

$\displaystyle 1$ $\displaystyle\to$ $\displaystyle\delta(1-z),$
$\displaystyle\frac{1}{N_{+}}\text{\ and\ }\frac{1}{N_{-}}$ $\displaystyle\to$
$\displaystyle\theta(1-z),$ $\displaystyle\frac{1}{N_{+}^{2}},\
\frac{1}{N_{-}^{2}},\text{\ and\ }\frac{1}{N_{+}N_{-}}$ $\displaystyle\to$
$\displaystyle\theta(1-z)\ln\frac{1}{z}.$ (44)

The expansion of the ${\bf n}_{T1,2}^{2}$ integrand of the function
$G_{ij}^{[m]}$ in $\gamma_{1}$ and $\gamma_{2}$ however has to be done in a
distributional sense because otherwise the order-by-order integrals will just
diverge at ${\bf n}_{T1,2}^{2}\to 0$. To this end, one isolates the ${\bf
n}_{T1,2}^{2}\to 0$ behaviour of the $f^{[m]}$ function888Which can be
understood as a test function, on which functionals (40), (41) and (42)
depend. in the integrands of eqns. (40), (41) and (42) as follows:

$\displaystyle f^{[m]}({\bf n}_{T1}^{2},{\bf n}_{T2}^{2},\phi,\nu)=$
$\displaystyle{\left\\{f^{[m]}({\bf n}_{T1}^{2},{\bf
n}_{T2}^{2},\phi,\nu)-f^{[m]}(0,{\bf n}_{T2}^{2},\phi,\nu)\theta(1-{\bf
n}_{T1}^{2})\right.}$ $\displaystyle{\left.-f^{[m]}({\bf
n}_{T1}^{2},0,\phi,\nu)\theta(1-{\bf
n}_{T2}^{2})+f^{[m]}(0,0,\phi,\nu)\theta(1-{\bf n}_{T1}^{2})\theta(1-{\bf
n}_{T2}^{2})\right\\}}$ $\displaystyle+{\left[f^{[m]}({\bf
n}_{T1}^{2},0,\phi,\nu)-f^{[m]}(0,0,\phi,\nu)\theta(1-{\bf
n}_{T1}^{2})\right]\theta(1-{\bf n}_{T2}^{2})}$
$\displaystyle+{\left[f^{[m]}(0,{\bf
n}_{T2}^{2},\phi,\nu)-f^{[m]}(0,0,\phi,\nu)\theta(1-{\bf
n}_{T2}^{2})\right]\theta(1-{\bf n}_{T1}^{2})}$
$\displaystyle{+f^{[m]}(0,0,\phi,\nu)\theta(1-{\bf n}_{T1}^{2})\theta(1-{\bf
n}_{T2}^{2}).}$ (45)

The expression in curly brackets in eqn. (45) tends to zero when ${\bf
n}^{2}_{T1}\to 0$ and/or ${\bf n}_{T2}^{2}\to 0$, so the factors $({\bf
n}_{T1}^{2})^{-1+\gamma_{1}}$ and $({\bf n}_{T2}^{2})^{-1+\gamma_{2}}$ in
front of it in eqns. (40) – (42) can be safely Taylor-expanded in
$\gamma_{1,2}$. The expression in the fourth line of eqn. (45) tends to zero
when ${\bf n}_{T1}^{2}\to 0$ while its dependence on ${\bf n}_{T2}^{2}$ is
just a step function. As a consequence, when this expression is substituted to
eqns. (40) – (42), the dependence on ${\bf n}_{T2}^{2}$ can be integrated out,
while the factor $({\bf n}_{T1}^{2})^{-1+\gamma_{1}}$ can be Taylor-expanded
in $\gamma_{1}$. For the term in the fifth line of eqn. (45), the situation is
symmetric, up to the replacement ${\bf n}_{T1}\leftrightarrow{\bf n}_{T2}$.
Finally, with the last term in eqn. (45), the integrations over ${\bf
n}_{T1}^{2}$ and ${\bf n}_{T2}^{2}$ can be performed straightforwardly. The
resulting expansions for the functions $G_{ij}^{[m]}$ in $\gamma_{1}$ and
$\gamma_{2}$ at $\nu=0$ are:

$\displaystyle G^{[m]}_{gg}$ $\displaystyle=$
$\displaystyle\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}+\gamma_{2}}\left[A^{[m]}_{0}+(\gamma_{1}+\gamma_{2})A^{[m]}_{1}+(\gamma_{1}^{2}+\gamma_{2}^{2})A^{[m]}_{2}+\gamma_{1}\gamma_{2}B^{[m]}_{2}+{\cal
O}(\gamma^{3})\right].$ (46) $\displaystyle G^{[m]}_{qg}$ $\displaystyle=$
$\displaystyle\frac{C_{F}}{C_{A}}\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{2}}\left\\{\left[\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}}-1\right]A^{[m]}_{0}+\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}}(\gamma_{1}A^{[m]}_{1}+\gamma_{1}^{2}A^{[m]}_{2})\right.$
$\displaystyle+$
$\displaystyle\left.\left[\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}}-1\right](\gamma_{2}A^{[m]}_{1}+\gamma_{2}^{2}A^{[m]}_{2})+\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}}\gamma_{1}\gamma_{2}B^{[m]}_{2}+{\cal
O}(\gamma^{3})\right\\},$ $\displaystyle G^{[m]}_{q\bar{q}}$ $\displaystyle=$
$\displaystyle\left(\frac{C_{F}}{C_{A}}\right)^{2}\left\\{\left[\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}}-1\right]\left[\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{2}}-1\right]A^{[m]}_{0}\right.$
$\displaystyle+$
$\displaystyle\left[\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{2}}-1\right]\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}}(\gamma_{1}A^{[m]}_{1}+\gamma_{1}^{2}A^{[m]}_{2})$
$\displaystyle+$
$\displaystyle\left.\left[\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}}-1\right]\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{2}}(\gamma_{2}A^{[m]}_{1}+\gamma_{2}^{2}A^{[m]}_{2})+\left(\frac{M^{2}}{\mu_{F}^{2}}\right)^{\gamma_{1}+\gamma_{2}}\gamma_{1}\gamma_{2}B^{[m]}_{2}+{\cal
O}(\gamma^{3})\right\\},$

where the Taylor expansion of the factors $(M^{2}/\mu_{F}^{2})^{\gamma_{1,2}}$
up to the second order have to be done as well. The coefficients
$A_{0,1,2}^{[m]}$ entering this expansion are:

$\displaystyle A_{0}^{[m]}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}f^{[m]}(0,0,\phi,0),$
(49) $\displaystyle A_{1}^{[m]}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\infty}\frac{d{\bf n}_{T1}^{2}}{{\bf
n}_{T1}^{2}}\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}\left[f^{[m]}({\bf
n}_{T1}^{2},0,\phi,0)-f^{[m]}(0,0,\phi,0)\theta(1-{\bf n}_{T1}^{2})\right],$
(50) $\displaystyle A_{2}^{[m]}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\infty}\frac{d{\bf n}_{T1}^{2}}{{\bf
n}_{T1}^{2}}\ln{\bf
n}_{T1}^{2}\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}\left[f^{[m]}({\bf
n}_{T1}^{2},0,\phi,0)-f^{[m]}(0,0,\phi,0)\theta(1-{\bf n}_{T1}^{2})\right],$
(51)

while, for the coefficient $B_{2}^{[m]}$, one has

$\displaystyle B^{[m]}_{2}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\infty}\frac{d{\bf n}_{T1}^{2}d{\bf
n}_{T2}^{2}}{{\bf n}_{T1}^{2}{\bf
n}_{T2}^{2}}\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}\left\\{f^{[m]}({\bf
n}_{T1}^{2},{\bf n}_{T2}^{2},\phi,0)-f^{[m]}(0,{\bf
n}_{T2}^{2},\phi,0)\theta(1-{\bf n}_{T1}^{2})\right.$ (52)
$\displaystyle\left.-f^{[m]}({\bf n}_{T1}^{2},0,\phi,0)\theta(1-{\bf
n}_{T2}^{2})+f^{[m]}(0,0,\phi,0)\theta(1-{\bf n}_{T1}^{2})\theta(1-{\bf
n}_{T2}^{2})\right\\}.$

Substituting the expansion (46) to eqn. (39) due to the mappings (44), one
obtains the following result for the $z$-dependent CF coefficient function in
$gg$ channel:

$\displaystyle\hat{\sigma}^{[m],\text{\ HEF}}_{gg}$ $\displaystyle=$
$\displaystyle\sigma^{[m]}_{0}\left\\{A^{[m]}_{0}\delta(1-z)+\frac{\alpha_{s}}{\pi}2C_{A}\left[A^{[m]}_{1}+A^{[m]}_{0}\ln\frac{M^{2}}{\mu_{F}^{2}}\right]\right.$
$\displaystyle+$
$\displaystyle\left.\left(\frac{\alpha_{s}}{\pi}\right)^{2}C_{A}^{2}\ln\frac{1}{z}\left[2A^{[m]}_{2}+B^{[m]}_{2}+4A^{[m]}_{1}\ln\frac{M^{2}}{\mu_{F}^{2}}+2A^{[m]}_{0}\ln^{2}\frac{M^{2}}{\mu_{F}^{2}}\right]+{\cal
O}(\alpha_{s}^{3})\right\\},$

where $\sigma^{[m]}_{0}$ has been defined in eqn. (18). The coefficients
$A_{0,1,2}^{[m]}$ and $B_{2}^{[m]}$ can be computed for various spin-orbital
states $m={}^{2S+1}L_{J}$ by substituting the HEF coefficient functions from
eqn. (21) to eqns. (49), (50), (51) and (52). The results of this computation
are listed in the table 1.

State | $A_{0}^{[m]}$ | $A_{1}^{[m]}$ | $A_{2}^{[m]}$ | $B_{2}^{[m]}$
---|---|---|---|---
${}^{1}S_{0}$ | $1$ | $-1$ | $\frac{\pi^{2}}{6}$ | $\frac{\pi^{2}}{6}$
${}^{3}S_{1}$ | $0$ | $1$ | $0$ | $\frac{\pi^{2}}{6}$
${}^{3}P_{0}$ | $1$ | $-\frac{43}{27}$ | $\frac{\pi^{2}}{6}+\frac{2}{3}$ | $\frac{\pi^{2}}{6}+\frac{40}{27}$
${}^{3}P_{1}$ | $0$ | $\frac{5}{54}$ | $-\frac{1}{9}$ | $-\frac{2}{9}$
${}^{3}P_{2}$ | $1$ | $-\frac{53}{36}$ | $\frac{\pi^{2}}{6}+\frac{1}{2}$ | $\frac{\pi^{2}}{6}+\frac{11}{9}$
Table 1: Coefficients entering into eqns. (2.4) – (55) for several states
$m={}^{2S+1}L_{J}$ of a $Q\bar{Q}$ pair. The highlighted numbers coincide with
the corresponding results from the table 1 in ref. Lansberg:2020ejc .

Similarly, from the expansions (2.4) and (2.4), one obtains the following
predictions for the small-$z$ behaviour of the quark-induced channels at NLO
and NNLO:

$\displaystyle\hat{\sigma}^{[m],\text{\ HEF}}_{qg}$ $\displaystyle=$
$\displaystyle\sigma_{0}^{[m]}\left\\{\frac{\alpha_{s}}{\pi}C_{F}\left[A^{[m]}_{1}+A^{[m]}_{0}\ln\frac{M^{2}}{\mu_{F}^{2}}\right]+\left(\frac{\alpha_{s}}{\pi}\right)^{2}C_{A}C_{F}\ln\frac{1}{z}\right.$
(54) $\displaystyle\times$
$\displaystyle\left.\left[A^{[m]}_{2}+B^{[m]}_{2}+3A^{[m]}_{1}\ln\frac{M^{2}}{\mu_{F}^{2}}+\frac{3}{2}A^{[m]}_{0}\ln^{2}\frac{M^{2}}{\mu_{F}^{2}}\right]\right\\},$
$\displaystyle\hat{\sigma}^{[m],\text{\ HEF}}_{q\bar{q}}$ $\displaystyle=$
$\displaystyle\sigma_{0}^{[m]}\left(\frac{\alpha_{s}}{\pi}\right)^{2}C_{F}^{2}\ln\frac{1}{z}\left[B^{[m]}_{2}+2A^{[m]}_{1}\ln\frac{M^{2}}{\mu_{F}^{2}}+A^{[m]}_{0}\ln^{2}\frac{M^{2}}{\mu_{F}^{2}}\right].$
(55)

Where the same coefficients $A_{0,1,2}^{[m]}$ and $B_{2}^{[m]}$ appear as in
the $gg$ case. The NLO parts of eqns. (2.4) and (54) coincide with the $z\to
0$ asymptotics of the full NLO results, obtained in refs. Mangano:1996kg ;
Lansberg:2020ejc , which is yet another non-trivial cross-check of the HEF
formalism.

Figure 2:  Dimensionless integrand function (43) averaged over $\phi$, at
${\bf n}_{T2}^{2}=0$, as a function of ${\bf n}_{T1}^{2}={\bf
q}_{T1}^{2}/(4m_{Q}^{2})$, normalised on its value at ${\bf n}_{T1}^{2}=0$.
The three dashed curves correspond to the coefficient functions for production
of $Q\bar{Q}[m]$-states with $m={}^{1}S_{0},$ ${}^{3}P_{0}$ and ${}^{3}P_{2}$
(section 2.2). The solid curve depicts the coefficient function for the
unbound $Q\bar{Q}$-pair production Ellis:1990hw , integrated over the phase
space of the $Q\bar{Q}$ pair. The short-dashed line represents the
$\overline{\text{MS}}$ subtraction term $\theta(1-{\bf n}_{T1}^{2})$.

Let us emphasise that the specific transverse-momentum dependence of the HEF
coefficient functions is crucial for the consistency of the HEF results with
the high-energy limit of the QCD scattering amplitudes. This is very different
from the case of TMD factorisation, which is applicable only at ${\bf
p}_{T}^{2}\ll M^{2}$. The NLO coefficient $A_{1}^{[m]}$ is determined by the
behaviour of the coefficient function for ${\bf p}_{T}^{2}\simeq{\bf
q}_{T1,2}^{2}\sim M^{2}$ as follows. In the integrand of eqn. (50), the step-
function $\theta(1-{\bf n}_{T1}^{2})$ is subtracted from the dimensionless
function $f^{[m]}({\bf n}_{T1}^{2},0,\phi,0)$, which describes the coefficient
function in the case when one gluon is off-shell with momenta ${\bf q}_{T,1}$
and the second one is on-shell ${\bf q}_{T,2}^{2}=0$. In the case of
quarkonia, the latter function starts to decrease already for ${\bf
n}_{T1}^{2}={\bf q}_{T1}^{2}/M^{2}<1$, so the subtraction leads to a large
negative contribution, as illustrated by the plot in the figure 2. For
comparison, the coefficient function for open heavy flavour production,
integrated over the phase-space of the final-state heavy quarks, given by
eqns. (4.8) – (4.10) in ref. Ellis:1990hw , is also plotted in figure 2. From
this figure, one immediately realises that the bound state is easily broken by
the transverse-momentum imbalance between the incoming off-shell gluons, which
leads to a quickly decreasing HEF coefficient function. This dependence is
quite different from the “expectation” of the $\overline{\text{MS}}$
scheme999See e.g. section 3.4 of ref. CollinsQCD  for the relation of the
$\overline{\text{MS}}$ subtraction of $1/\epsilon$ singularity with the
transverse momentum cut off., according to which the coefficient function
should behave roughly like $\theta(1-{\bf q}_{T1,2}^{2}/M^{2})$ thus leading
to the consequent oversubtraction of the collinear behaviour. For the case of
open heavy flavour, integrated over its invariant mass, the situation is
opposite and the coefficient function has a substantial tail at ${\bf
q}_{T1,2}^{2}>M^{2}$, which leads to large positive NLO corrections at high
energy.

###  2.5 Resummation by a $\mu_{F}$ choice

The connection between the HEF resummation and the $\mu_{F}$-scale
optimisation approach (the $\hat{\mu}_{F}$ prescription) proposed in ref.
Lansberg:2020ejc  is most easily illustrated by eqn. (40). Substituting the
scale choice (4), one obtains:

$\displaystyle
G^{[m]}_{gg}\left(\gamma_{N_{+}},\gamma_{N_{-}},\nu,\frac{M^{2}}{\hat{\mu}_{F}^{2}}\right)$
$\displaystyle=$
$\displaystyle\exp\left[-A_{1}^{[m]}(\gamma_{N_{+}}+\gamma_{N_{-}})\right]$
$\displaystyle\times$
$\displaystyle\int\limits_{0}^{\infty}\frac{\gamma_{N_{+}}\gamma_{N_{-}}d{\bf
n}_{T1}^{2}d{\bf n}_{T2}^{2}}{({\bf n}_{T1}^{2})^{1-\gamma_{N_{+}}}({\bf
n}_{T2}^{2})^{1-\gamma_{N_{-}}}}\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}f^{[m]}({\bf
n}_{T1}^{2},{\bf n}_{T2}^{2},\phi,\nu).$

The exponent in the first line of eqn. (2.5), which arose from the scale
choice (4), resums a series of corrections proportional to
$\gamma_{N_{\pm}}^{n}=\hat{\alpha}_{s}^{n}/N_{\pm}^{n}$, which belong to the
LL($\ln(1/z)$)-approximation. This resummation would be equivalent to that
performed by the HEF only if this exponent is able to cancel the
$\gamma_{N_{\pm}}$ dependence of the integral in the second line of eqn.
(2.5). If this was the case then all LL($\ln(1/z)$) corrections would be
removed from the CF coefficient function and absorbed into the scale evolution
of PDFs. However, such a perfect cancellation is possible only if the
dimensionless function $f^{[m]}$ complies to the relation:

$\int\limits_{0}^{2\pi}\frac{d\phi}{2\pi}f^{[m]}({\bf n}_{T1}^{2},{\bf
n}_{T2}^{2},\phi,\nu)=\theta\left(e^{A_{1}^{[m]}}-{\bf
n}_{T1}^{2}\right)\theta\left(e^{A_{1}^{[m]}}-{\bf n}_{T2}^{2}\right),$ (57)

which when substituted to the eqn. (2.5) leads to
$G_{gg}^{[m]}\left(\gamma_{N_{+}},\gamma_{N_{-}},\nu,{M^{2}}/{\hat{\mu}_{F}^{2}}\right)=1$.
Neither the coefficient functions for quarkonium production listed in the
section 2.2, nor any other coefficient functions for physical processes known
to the authors provide a sharp cut-off in the transverse momentum as in eqn.
(57). Instead, they always smoothly depend on transverse momenta, as e.g. in
the figure 2. Therefore, the HEF is not equivalent to the $\hat{\mu}_{F}$
prescription. By construction, the $\hat{\mu}_{F}$ prescription takes into
account the coefficient $A_{1}^{[m]}$ and, hence, it is correct up to NLO in
$\alpha_{s}$. However, already at NNLO in $\alpha_{s}$, the results of HEF and
the $\hat{\mu}_{F}$ prescription differ. This is clear because the NNLO
coefficients obtained from HEF (table 1) contain the transcendental number
$\pi^{2}/6$, which can not arise from the expansion of the exponent in the
first line of eqn. (2.5).

##  3 NLO+DLA matching in $z$-space

Before discussing various approaches to match the DLA HEF and NLO CF
contributions into a single prediction for the total cross section, let us set
the baseline by showing how the NLO CF cross section with the “canonical”
scale choice $\mu_{F}=\mu_{R}=M$ depends on the energy and how large the
associated scale uncertainty is. The plots of figure 3 show the LO and NLO CF
predictions for the total cross section of the production of
${}^{1}S_{0}^{[1]}$ states with masses of 3 and 9.4 GeV. The scale
uncertainties, depicted as solid or shaded bands in all these plots, has been
estimated using the five-point scale-variation procedure. The latter procedure
consists in evaluating the cross section for each value of the energy with the
scale choice $\mu_{F}=2^{\zeta_{1}}M$ and $\mu_{R}=2^{\zeta_{2}}M$ for
$(\zeta_{1},\zeta_{2})\in\\{(0,\pm 1),(\pm 1,0)\\}$ and taking, as an
uncertainty estimate, the largest positive or negative deviation of the
obtained value form the cross section with the default scale choice.

Figure 3: Comparison of the energy dependence of the hadroproduction cross
sections of a $Q\bar{Q}[{}^{1}S_{0}^{(1)}]$ state with $M=3$ GeV (left panel)
and $M=9.4$ GeV (right panel) in CF at LO (grey curve) and NLO(blue curve).
The central member of the CT18NLO PDF set Hou:2019efy  has been used. The
dashed line depicts the central prediction of the NLO calculation using the
$\hat{\mu}_{F}$ prescription (Equation 4) Lansberg:2020ejc . The LDME
$\langle{\cal O}[{}^{1}S_{0}^{(1)}]\rangle$ was set to 1 GeV3 in both cases
for illustration purposes.

From the left panel of figure 3, one can see that, for $c\bar{c}$ states, the
instability of the cross section at high energy is dramatic. Above
$\sqrt{s}=1$ TeV, such a computation does not show any more predictive power
and, above $\sqrt{s}=100$ GeV, the NLO $K$ factor, defined as the ratio of the
NLO to the LO cross section, decreases with energy, which is very different
from the behaviour of the matched result as we will see. For bottomonia (see
the plot in the right panel of figure 3), the situation is much better: the
scale uncertainty of the NLO calculation is significantly smaller than the LO
one all the way up to $\sqrt{s}\sim 100$ TeV. However, the decrease of the NLO
$K$ factor is evident even at this scales. The dependence of this behaviour on
the PDF choice has been analysed in detail in ref. Lansberg:2020ejc .

A specific factorisation scale choice (4) has been proposed in ref.
Lansberg:2020ejc  as a possible resolution of the problem of high-energy
instability of quarkonium hadroproduction cross sections. As we have shown in
section 2.5, this prescription is equivalent to the resummation of some
LL($\ln 1/z$)-terms, which is correct at NLO but fails at higher orders. In
practice, for $\eta_{Q}$ hadroproduction, this procedure leads to
significantly smaller cross section at high energies, shown by dashed lines in
figure 3, than that predicted by LO CF at the default scale. Indeed the PDFs
are less evolved thus smaller. In this regards, we stress that, for other
processes like $H^{0}$ production, $\hat{\mu}_{F}$ can be larger than the mass
of the produced particle Lansberg:2020ejc  and will lead to larger cross
sections with further evolved gluon PDFs.

The DL HEF resummation proposed in section 2.3 is compatible with the
factorisation scheme and the evolution of the usual NLO and NNLO PDFs and
correctly reproduces the LL($\ln 1/z$) terms in the CF coefficient function up
to NNLO in $\alpha_{s}$, as it was shown in sections 2.3 and 2.4. However this
resummation is applicable only in the region of $z\ll 1$. Power corrections
${\cal O}(z)$ to the CF coefficient function are missing in the HEF
calculation with any NkLL accuracy, while they are equally importantas the
logarithmic terms at $z\sim 1$. For the value of the total cross section, the
whole range of $z$ from 0 to 1 contributes due to the integration over $z$ in
eqn. (1). In the following subsections, we will propose two approaches to
match these DL HEF and NLO CF results in $z$ space and will compare the
corresponding numerical results.

###  3.1 The subtractive-matching prescription

The simplest matching prescription between NLO CF and DLA HEF calculations
consists in the subtraction of the $z\to 0$ asymptotics of the NLO CF
coefficient function, to avoid double-counting it with the HEF contribution:

$\sigma^{[m]}_{\rm NLO+HEF}=\sigma^{[m]}_{\text{LO
CF}}+\sum\limits_{i,j=q,\bar{q},g}\int\limits_{z_{\min}}^{1}\frac{dz}{z}\
\left[\check{\sigma}^{[m],ij}_{\text{HEF}}(z)+\hat{\sigma}^{[m],ij}_{\text{NLO
CF}}(z)-\hat{\sigma}_{\text{NLO CF}}^{[m],ij}(0)\right]{\cal L}_{ij}(z),$ (58)

where $\check{\sigma}^{[m],gg}_{\text{HEF}}$ is defined in eqn. (13) while,
for other partonic channels,
$\check{\sigma}^{[m],ij}_{\text{HEF}}=\hat{\sigma}^{[m],ij}_{\text{HEF}}$ and
$\hat{\sigma}^{[m],ij}_{\text{NLO CF}}(z)$ includes only the ${\cal
O}(\alpha_{s})$ CF terms from the original papers Kuhn:1992qw ;
Petrelli:1997ge .

The corresponding numerical results for $Q\bar{Q}[{}^{1}S_{0}^{(1)}]$ states
with masses of 3 and 9.4 GeV are shown in the figure 4. One can see that the
high-energy instability of the NLO cross section has gone away and the ratio
of the matched cross section to the LO one becomes almost constant at high
energy, which signal that energy dependence is now mostly driven by the PDFs,
not by the hard-scattering coefficient.

Figure 4: Comparison of the energy dependence of the hadroproduction cross
sections of a $Q\bar{Q}[{}^{1}S_{0}^{[1]}]$ state with $M=3$ GeV (left panel)
and $M=9.4$ GeV (right panel) for the matched NLO CF + DLA HEF calculations
using the subtractive (orange curve and right-shaded band) and InEW
matchings(red curve and left-shaded band). The dashed line depicts the central
prediction of the NLO calculation using the $\hat{\mu}_{F}$ prescription
(Equation 4) Lansberg:2020ejc . The LDME $\langle{\cal
O}[{}^{1}S_{0}^{[1]}]\rangle$ was set to 1 GeV3.

###  3.2 The inverse-error-weighting matching

The InEW matching method has been proposed in ref. Echevarria:2018qyi  to
match the $|{\bf p}_{Tl^{+}l^{-}}|$ spectrum of Drell-Yan lepton pairs
predicted by the TMD factorisation for $|{\bf p}_{Tl^{+}l^{-}}|\ll
M_{l^{+}l^{-}}$ with the NLO CF calculation of the same spectrum, which is
free from large logarithms of $|{\bf p}_{Tl^{+}l^{-}}|/M_{l^{+}l^{-}}$ for
$|{\bf p}_{Tl^{+}l^{-}}|\gtrsim M_{l^{+}l^{-}}$. Following the same ideas, we
introduce smooth weight functions $w_{\rm HEF}(z)$ and $w_{\rm CF}(z)$ into
our cross section formula as follows:

$\displaystyle\sigma^{[m]}_{\rm NLO+HEF}=\sigma^{[m]}_{\text{LO CF}}+$ (59)
$\displaystyle+\sum\limits_{i,j=q,\bar{q},g}\int\limits_{z_{\min}}^{1}dz\ $
$\displaystyle\left\\{\left[\check{\sigma}^{[m],ij}_{\text{HEF}}(z)\frac{{\cal
L}_{ij}(z)}{z}\right]w_{\rm
HEF}^{ij}(z)+\left[\hat{\sigma}^{[m],ij}_{\text{NLO CF}}(z)\frac{{\cal
L}_{ij}(z)}{z}\right]w_{\rm CF}^{ij}(z)\right\\},$

where we emphasise that the weights are introduced on the level of the $z$
integrand of the cross section, so that the distributions $1/(1-z)_{+}$ and
$\ln(1-z)/(1-z)_{+}$, contained in the $\hat{\sigma}^{[m]}_{\text{NLO
CF},ij}$, have already “acted” on the function ${\cal L}_{ij}(z)/z$. The
weights should suppress one of the contributions in the region where it is
unreliable, so naturally they are chosen to be proportional to the inverse
squared errors of each contribution and to add up to one:

$w^{ij}_{\rm HEF}(z)=\frac{[\Delta\sigma_{\rm
HEF}^{ij}(z)]^{-2}}{[\Delta\sigma_{\rm HEF}^{ij}(z)]^{-2}+[\Delta\sigma_{\rm
CF}^{ij}(z)]^{-2}},\ w^{ij}_{\rm CF}(z)=1-w^{ij}_{\rm HEF}(z),$ (60)

where, by $\Delta\sigma^{ij}_{\rm HEF/CF}$, we refer to an estimate of the
theoretical uncertainty of the corresponding contributions to the $z$
integrand of the total cross section, which should correctly capture the fact
that the HEF or the CF contributions respectively become unreliable in the
$z\to 1$ or $z\to 0$ limits. The integrand uncertainty due to the matching
procedure then follows from:

$\Delta\sigma^{ij}_{\rm Match.}(z)=\left([\Delta\sigma_{\rm
HEF}^{ij}(z)]^{-2}+[\Delta\sigma_{\rm CF}^{ij}(z)]^{-2}\right)^{-1/2},$ (61)

which, upon integration over $z$, gives the total cross section uncertainty
due to matching procedure.

The errors entering eqns. (60) and (61) can be estimated as follows. The NLO
CF cross section contains no information about the small-$z$ behaviour of
$\hat{\sigma}^{[m],ij}(z)$ beyond ${\cal O}(\alpha_{s})$, while the LL HEF
calculation provides the $\alpha_{s}^{2}\ln 1/z$ and higher-order LLA terms.
Hence one can take the ${\cal O}(\alpha_{s}^{2})$ expansion of the HEF cross
section, obtained in Sec. 2.4 as an estimate for the error of the NLO CF
integrand. We should also take into account the fact, that the constant ${\cal
O}(\alpha_{s}^{2})$ term of the $z\to 0$ asymptotics of the CF coefficient
function can not be obtained from the LL HEF calculation, because it belongs
to the NLL approximation of the HEF. We estimate the coefficient in front of
this term as equal to $1$ and treat it as an independent source of
uncertainty. Adding the constant term to the $\alpha_{s}^{2}\log 1/z$
contribution in quadrature we obtain the following uncertainty estimate for
the NLO CF contribution:

$\Delta\sigma_{\rm CF}^{ij}(z)=\alpha_{s}^{2}\sigma_{0}^{[m]}\frac{{\cal
L}_{ij}(z)}{z}\sqrt{\left[\frac{C_{\rm
LL}^{ij}}{\pi}\ln\frac{1}{z}\right]^{2}+1},$ (62)

where

$\displaystyle C^{gg}_{LL}$ $\displaystyle=$ $\displaystyle
C_{A}^{2}\left[2A_{2}+B_{2}+4A_{1}\ln\frac{M^{2}}{\mu_{F}^{2}}+2A_{0}\ln^{2}\frac{M^{2}}{\mu_{F}^{2}}\right],$
(63) $\displaystyle C^{qg}_{LL}$ $\displaystyle=$ $\displaystyle
C_{A}C_{F}\left[A_{2}+B_{2}+3A_{1}\ln\frac{M^{2}}{\mu_{F}^{2}}+\frac{3}{2}A_{0}\ln^{2}\frac{M^{2}}{\mu_{F}^{2}}\right],$
(64) $\displaystyle C^{q\bar{q}}_{LL}$ $\displaystyle=$ $\displaystyle
C_{F}^{2}\left[B_{2}+2A_{1}\ln\frac{M^{2}}{\mu_{F}^{2}}+A_{0}\ln^{2}\frac{M^{2}}{\mu_{F}^{2}}\right],$
(65)

and the coefficients $A_{0}$, $A_{1,2}$, $B_{2}$ can be found in table 1.

Conversely, the HEF result is valid only at $z\ll 1$ and contains no
information about the ${\cal O}(z)$ power corrections to the CF coefficient
function. On the other hand, the NLO CF coefficient function is an exact
function of $z$ at NLO in $\alpha_{s}$. Hence one can take the NLO CF hard-
scattering coefficient with its $z\to 0$ asymptotics subtracted, as an
estimate of the missing ${\cal O}(z)$ power corrections in HEF, i.e. of the
error of the latter. We also take into account the unknown $z\to 1$
asymptotics of the NNLO CF coefficient function as ${\cal O}(\alpha_{s}^{2})$
term and add this estimate to the total error of the HEF contribution in
quadrature as follows:

$\Delta\sigma_{\rm
HEF}^{ij}(z)=\sqrt{\left[\left(\hat{\sigma}^{[m],ij}_{\text{NLO
CF}}(z)-\hat{\sigma}_{\text{NLO CF}}^{[m],ij}(0)\right)\frac{{\cal
L}_{ij}(z)}{z}\right]^{2}+\left[\alpha_{s}^{2}\sigma^{[m]}_{0}\frac{{\cal
L}_{ij}(z)}{z}\right]^{2}}.$ (66)

Using formulas (62) and (66), one can construct weights $w_{\rm CF}$ and
$w_{\rm HEF}$ with expected behaviour: e.g. $w^{(gg)}_{\rm CF}(z)$ increases
between 0 and 1 when $z\to 1$, figure 5. The latter property is ensured by the
behaviour of errors (62) and (66). The CF error (62) increases towards $z\to
0$, while HEF error (66) increases towards $z=1$. Among three weights shown in
the figure 5, only the weight for $q\bar{q}$ channel does not reach $1$ when
$z\to 1$. This behaviour is expectable, because the NLO coefficient function,
and hence the error (66) has no $\ln(1-z)$ enhancement at $z\to 1$ for this
channel, compared to other channels. But this feature of the weight function
of the $q\bar{q}$ channel is not important numerically, because the $q\bar{q}$
contribution amounts to ${\cal O}(1\%)$ of the total cross section at high
energy.

Figure 5: Plots of the InEW weights of the CF contribution for different
channels as function of $z$ with $\mu_{F}=\mu_{R}=M$.

In figures 6, 7 and 8 several plots of the $z$ integrand of eqn. (59) are
shown as red solid lines, using weights obtained above. One can see that the
InEW procedure provides a smooth interpolation between the HEF curves at $z\ll
1$ and the CF curves at $z$ closer to 1.

The corresponding curves for the subtractive matching are shown for comparison
by the red dotted lines. From figures 6, 7 and 8, the problem of subtractive
matching is evident: it subtracts the $z\to 0$ asymptotics of the NLO CF
coefficient function at all values of $z$, thus introducing the unphysical
modification of the $z\to 1$ behaviour. The InEW matching, on the other hand,
smoothly connects the asymptotic regions.

Figure 6: Matching plots for the $gg$-channel contribution to the
${}^{1}S_{0}^{[1]}$-state hadroproduction with $M=3$ GeV. The solid curve
depicts the $z$-integrand of eqn. (59), the dashed curve the HEF contribution
(without the InEW weight), the dash-dotted curve the NLO CF contribution
(without the InEW weight), and the dotted red line the integrand of eqn. (58),
i.e. the result of the subtractive matching prescription for comparison. The
plots of the InEW weights are shown in the bottom inset, while the matching
uncertainty (61) is shown as the solid band. The LDME $\langle{\cal
O}[{}^{1}S_{0}^{[1]}]\rangle$ was set to 1 GeV3.

Figure 7: Matching plots for $qg$-channel to the ${}^{1}S_{0}^{[1]}$-state
hadroproduction with $M=3$ GeV. The integrand of eqn. (58), i.e. subtractive
matching prescription is shown for comparison by the dotted red line. The
plots of InEW weights are shown in the bottom inset, while the matching
uncertainty (61) is shown as the solid band. The LDME $\langle{\cal
O}[{}^{1}S_{0}^{[1]}]\rangle$ was set to 1 GeV3.

Figure 8: Matching plots for $q\bar{q}$-channel to the
${}^{1}S_{0}^{[1]}$-state hadroproduction with $M=3$ GeV. The integrand of
eqn. (58), i.e. subtractive matching prescription is shown for comparison by
the dotted red line. The plots of InEW weights are shown in the bottom inset,
while the matching uncertainty (61) is shown as the solid band. The LDME
$\langle{\cal O}[{}^{1}S_{0}^{[1]}]\rangle$ was set to 1 GeV3.

The results for the matched cross section in the InEW procedure are shown by
red solid curves on figure 4. The matching uncertainty, estimated by eqn. (61)
is shown separately as the green solid band in the figure 4 and it is
significantly smaller than the residual scale uncertainty, shown as the shaded
band. The scale uncertainty band is slightly smaller in the InEW matching
case.

Finally, figure 9 illustrates how the matched cross section varies depending
on the choice of the collinear PDF. Since our computations are rather computer
intensive, only the central members of the CT18NLO Hou:2019efy ,
MSHT20nlo_as118 Bailey:2020ooq , NNPDF31_nlo_as_0118 NNPDF:2017mvq  and
NNPDF31sx_nlonllx_as_0118 Ball:2017otu  PDF sets, as implemented in the LHAPDF
library Buckley:2014ana  has been used for the comparison. The first three
PDFs evolve according to the usual NLO DGLAP splitting functions and thus fit
seamlessly into our NLO+DLA matching scheme. The last PDF set Ball:2017otu
contains the NLL($\ln(1/z)$) resummation and is included in figure 9 to show
that the corresponding cross section, obtained in the NLO+DLA
$\ln(1/z)$-resummation scheme, is consistent with the other results within
uncertainties.

On one hand, figure 9 shows that there is a significant spread in the
charmonium production cross sections at high energy, due to the PDF
uncertainties at small values of $x$. At $\sqrt{s}=13$ TeV, the central values
of NLO+LL($\ln(1/z)$) predictions obtained with different PDFs vary as much as
by factor 1.5 up and down from our default CT18NLO-based estimations. This
demonstrates a potential usefulness of this observable to constrain PDFs at
small $x$. On the other hand, the theoretical uncertainties of our
predictions, dominated by the scale-variation uncertainties shown in figure 9,
are still larger than the spread between central curves, so to reliably
discriminate between them, the scale uncertainties have to be significantly
reduced which we expected to be realised once a NLL computation is available.

Figure 9: Dependence of the matched cross section on the choice of collinear
PDF. The LDME $\langle{\cal O}[{}^{1}S_{0}^{(1)}]\rangle$ was set to 1 GeV3.

##  4 Conclusions and outlook

In the present paper, we have performed an exploratory study of the effects of
High-Energy resummation on heavy-quarkonium hadroproduction cross sections. We
have introduced the DLA for the resummation contributions, which is consistent
with the fixed-order DGLAP evolution of PDFs and have shown that matching the
NLO CF and the DL resummed HEF results solves the problem of the negative NLO
CF cross sections at high energy. This opens up the possibility to use this
observable to constrain PDFs at small-$x$, however the scale uncertainties of
our predictions have to be significantly reduced. To this end, advancing such
calculation to a complete NLO+NLL($\ln(1/z)$) accuracy is required, which is a
work in progress. We also emphasize, that performing the matching between NLO
CF and HEF contributions in $z$-space (as opposed to matching in Mellin-
conjugate $N$-space, as in Ref. Catani:1992rn ) is a simple and flexible
procedure, which allows one to estimate the uncertainty of the cross section
due to the matching. When the NLL($\ln(1/z)$) calculation in HEF will become
available, it will be easy to match it with the NLO or possibly even NNLO CF
results using the $z$-space technique.

###### Acknowledgements.

We thank Renaud Boussarie, Francesco Hautmann and Samuel Wallon for useful
discussions as well as Kate Lynch and Yelyzaveta Yedelkina for comments on the
manuscript. This project has received funding from the European Union’s
Horizon 2020 research and innovation programme under grant agreement
STRONG–2020 No 824093 in order to contribute to the EU Virtual Access
NLOAccess (VA1-WG10) & to the Joint Research activity “Fixed Target
Experiments at the LHC” (JRA2). This project has also received funding from
the Agence Nationale de la Recherche (ANR) via the grant ANR-20-CE31-0015
(“PrecisOnium”) and via the IDEX Paris-Saclay “Investissements d’Avenir”
(ANR-11-IDEX-0003-01) through the GLUODYNAMICS project funded by the “P2IO
LabEx (ANR-10-LABX-0038)”. This work was also partly supported by the French
CNRS via the GDR QCD, via the IN2P3 project GLUE@NLO, via the Franco-Polish
EIA (Gluegraph). M.A.O.’s work was partly supported by the ERC grant 637019
“MathAm”.

##  Appendix A Numerics in $z$-space

###  A.1 Method 1: regularised resummation factor

This method is based on the following regularised form of the resummation
factor, where the small-${\bf q}_{T}$ oscillations are replaced by a constant
in ${\bf q}_{T}$-behaviour for ${\bf q}_{T}^{2}<{\bf q}_{T0}^{2}$ in such a
way that the normalisation condition (37) still holds exactly:

$\widetilde{\cal C}^{\rm\ DLA}_{gg}(z,{\bf
q}_{T}^{2},\mu_{F}^{2})=\left\\{\begin{matrix}\frac{1}{{\bf
q}_{T0}^{2}}\left[\delta(1-z)-{\cal F}\left(\frac{{\bf
q}^{2}_{T0}}{\mu_{F}^{2}},z\right)\right],&\text{ if }{\bf q}_{T}^{2}<{\bf
q}_{T0}^{2},\\\ {\cal C}_{gg}^{\rm DLA}(z,{\bf q}_{T}^{2},\mu_{F}^{2}),&\text{
if }{\bf q}_{T}^{2}>{\bf q}_{T0}^{2},\end{matrix}\right.$ (67)

where

${\cal F}\left(\frac{{\bf
q}_{T0}^{2}}{\mu_{F}^{2}},z\right)=\int\limits_{0}^{L_{\max}}dL\
J_{0}\left(2\sqrt{L\ln\frac{1}{z}}\right)=L_{\max}\
{}_{0}\tilde{F}_{1}\left(2,-L_{\max}\ln\frac{1}{z}\right),$ (68)

with $L_{\max}=\hat{\alpha}_{s}\ln({{\bf q}_{T0}^{2}}/{\mu_{F}^{2}})$ and
${}_{0}\tilde{F}_{1}$, which is the regularised confluent hypergeometric
function. In the ${\bf q}_{T0}^{2}\to 0$ limit, the function ${\cal F}$ turns
into $\delta(1-z)$ but this limit is approached very slowly. So, for all
reasonable values of $|{\bf q}_{T0}|$, even as small as $10^{-3}$ GeV, the
function (68) is still rather smoothly dependent on $z$ and does not cause any
numerical problems. For the quark-induced channel, the regularised resummation
factor:

$\widetilde{\cal C}^{\rm\ DLA}_{gq}(z,{\bf
q}_{T}^{2},\mu_{F}^{2})=\frac{C_{F}}{C_{A}}\left\\{\begin{matrix}-\frac{1}{{\bf
q}_{T0}^{2}}{\cal F}\left(\frac{{\bf
q}^{2}_{T0}}{\mu_{F}^{2}},z\right),&\text{ if }{\bf q}_{T}^{2}<{\bf
q}_{T0}^{2},\\\ {\cal C}^{\rm DLA}_{gg}(z,{\bf q}_{T}^{2},\mu_{F}^{2}),&\text{
if }{\bf q}_{T}^{2}>{\bf q}_{T0}^{2},\end{matrix}\right.$ (69)

satisfies the normalisation condition (38).

After substituting the regularised resummation factors (67) or (69) into eqn.
(12), one splits the ${\bf q}_{T1,2}^{2}$-integrations into regions with ${\bf
q}_{T1,2}^{2}<{\bf q}_{T0}^{2}$ and ${\bf q}_{T1,2}^{2}>{\bf q}_{T0}^{2}$. If,
for example ${\bf q}_{T1}^{2}<{\bf q}_{T0}^{2}$, then this transverse momentum
can be neglected in all the rest of the integrand in eqn. (12) because the
dependence of the HEF coefficient function as well as the dependence of the
factor $1/M_{T}^{4}$ on ${\bf q}_{T1}^{2}$ is smooth. This approximation
allows one to trivially integrate out ${\bf q}_{T1}^{2}$ in the region ${\bf
q}_{T1}^{2}<{\bf q}_{T0}^{2}$ introducing the error which scales as ${\cal
O}({\bf q}_{T0}^{2}/M^{2})$ and can be made negligible by choosing a
reasonably small value of ${\bf q}_{T0}^{2}$. In the region ${\bf
q}_{T1}^{2}>{\bf q}_{T0}^{2}$ no approximations are made. The ${\bf
q}_{T2}^{2}$-integration can be decomposed in the same manner and the region
where both ${\bf q}_{T1}^{2}<{\bf q}_{T0}^{2}$ and ${\bf q}_{T2}^{2}<{\bf
q}_{T0}^{2}$ produces the LO CF term in eqn. (13), while the LO-subtracted
partonic cross section in the method under consideration is decomposed as
follows:

$\check{\sigma}_{gg}^{[m],\ {\rm HEF}}=\check{\sigma}_{\rm
2UI}^{[m]}+\check{\sigma}_{\rm 2F}^{[m]}+2\left[\check{\sigma}_{\rm
1UI-C}^{[m]}-\check{\sigma}_{\rm 1UI-CF}^{[m]}-\check{\sigma}_{\rm
F}^{[m]}\right],$ (70)

where the separate contributions depend on the value of the cut ${\bf q}_{T0}$
but this dependence should cancel in their sum for sufficiently small ${\bf
q}_{T0}^{2}\ll M^{2}$. For the quark-induced channels, the resummed partonic
coefficient can be calculated in terms of the same contributions as:

$\displaystyle\hat{\sigma}_{gq}^{[m],\ {\rm HEF}}$ $\displaystyle=$
$\displaystyle\frac{C_{F}}{C_{A}}\left[\check{\sigma}_{\rm
2UI}^{[m]}+\check{\sigma}_{\rm 2F}^{[m]}+\check{\sigma}_{\rm
1UI-C}^{[m]}-2\check{\sigma}_{\rm 1UI-CF}^{[m]}-\check{\sigma}_{\rm
F}^{[m]}\right],$ (71) $\displaystyle\hat{\sigma}_{q\bar{q}}^{[m],\ {\rm
HEF}}$ $\displaystyle=$
$\displaystyle\left(\frac{C_{F}}{C_{A}}\right)^{2}\left[\check{\sigma}_{\rm
2UI}^{[m]}+\check{\sigma}_{\rm 2F}^{[m]}-2\check{\sigma}_{\rm 1UI-
CF}^{[m]}\right].$ (72)

The doubly-unintegrated contribution, $\check{\sigma}_{\rm 2UI}^{[m]}$,
entering eqns. (70), (71) and (72) is just eqn. (12) with both ${\bf
q}_{T1,2}^{2}>{\bf q}_{T0}^{2}$ and the resummation factor (36). Other
contributions are calculated as follows:

$\displaystyle\check{\sigma}^{[m]}_{\rm 1UI-CF}$ $\displaystyle=$
$\displaystyle\int\limits_{-\infty}^{+\infty}d\eta\int\limits_{{\bf
q}_{T0}^{2}}^{\infty}d{\bf q}_{T1}^{2}{\cal C}^{\rm
DLA}_{gg}\left(\sqrt{z}\frac{M_{T1}}{M}e^{\eta},{\bf
q}_{T1}^{2},\mu_{F},\mu_{R}\right)$ (73) $\displaystyle\times$
$\displaystyle{\cal F}\left(\frac{{\bf
q}_{T0}^{2}}{\mu_{F}^{2}},\sqrt{z}\frac{M_{T1}}{M}e^{-\eta}\right)\int\limits_{0}^{2\pi}\frac{d\phi}{2}\
\frac{H({\bf q}^{2}_{T1},0,\phi)}{M_{T1}^{4}},$
$\displaystyle\check{\sigma}^{[m]}_{\rm 1UI-C}$ $\displaystyle=$
$\displaystyle\int\limits_{{\bf q}_{T0}^{2}}^{\infty}d{\bf q}_{T1}^{2}\ {\cal
C}^{\rm DLA}_{gg}\left(z\frac{M_{T1}^{2}}{M^{2}},{\bf
q}_{T1}^{2},\mu_{F},\mu_{R}\right)\int\limits_{0}^{2\pi}\frac{d\phi}{2}\
\frac{H({\bf q}^{2}_{T1},0,\phi)}{M_{T1}^{4}},$ (74)
$\displaystyle\check{\sigma}^{[m]}_{\rm 2F}$ $\displaystyle=$
$\displaystyle\int\limits_{-\infty}^{\infty}d\eta\ {\cal F}\left(\frac{{\bf
q}_{T0}^{2}}{\mu_{F}^{2}},\sqrt{z}e^{\eta}\right){\cal F}\left(\frac{{\bf
q}_{T0}^{2}}{\mu_{F}^{2}},\sqrt{z}e^{-\eta}\right)\int\limits_{0}^{2\pi}\frac{d\phi}{2}\
\frac{H(0,0,\phi)}{M^{4}},$ (75) $\displaystyle\check{\sigma}^{[m]}_{\rm F}$
$\displaystyle=$ $\displaystyle{\cal F}\left(\frac{{\bf
q}_{T0}^{2}}{\mu_{F}^{2}},z\right)\int\limits_{0}^{2\pi}\frac{d\phi}{2}\
\frac{H(0,0,\phi)}{M^{4}},$ (76)

where $M_{T1}^{2}=M^{2}+{\bf q}_{T1}^{2}$.

###  A.2 Method 2: small-${\bf q}_{T}$ subtraction

The approach described in this section has been inspired by the treatment of
the same problem in the ref. Collins:1991ty . The method relies on the same
kind of decomposition of the transverse-momentum dependence of the integrand
as in eqn. (45), but now we perform it to isolate the LO contribution, using
the normalisation condition (37). To save space, let us introduce the short-
hand notation for the integrand of eqn. (12) :

$\displaystyle{\cal J}(\xi_{1},\xi_{2})$ $\displaystyle=$ $\displaystyle{\cal
C}_{gg}\left(\sqrt{z}\frac{M_{T}(\xi_{1},\xi_{2})}{M}e^{\eta},{\bf
q}_{T1}^{2},\mu_{F},\mu_{R}\right){\cal
C}_{gg}\left(\sqrt{z}\frac{M_{T}(\xi_{1},\xi_{2})}{M}e^{-\eta},{\bf
q}_{T2}^{2},\mu_{F},\mu_{R}\right)$ (77) $\displaystyle\times$
$\displaystyle\frac{H(\xi_{1}^{2}{\bf q}^{2}_{T1},\xi_{2}^{2}{\bf
q}^{2}_{T2},\phi)}{M_{T}^{4}(\xi_{1},\xi_{2})},$

where $M_{T}^{2}(\xi_{1},\xi_{2})=M^{2}+(\xi_{1}{\bf q}_{T1}+\xi_{2}{\bf
q}_{T2})^{2}$ and variables $\xi_{1,2}=\\{0,1\\}$ are introduced for us to be
able to turn on and off the “smooth part” of the dependence of the integrand
on transverse momenta, while the dependence of ${\cal J}$ on ${\bf
q}^{2}_{T1,2}$, $\phi$ and $\eta$ is implicit. Similarly to eqn. (45) we
decompose

${\cal J}(1,1)={\cal J}_{1}+{\cal J}_{2}+{\cal J}_{3}+{\cal
J}(0,0)\theta(\mu_{F}^{2}-{\bf q}_{T1}^{2})\theta(\mu_{F}^{2}-{\bf
q}_{T2}^{2}),$

where the last term of this decomposition, when substituted to eqn. (12) gives
us $\sigma_{\rm LO}^{[m]}\delta(1-z)$ due to normalisation condition (37),
while the terms ${\cal J}_{1,2,3}$ are:

$\displaystyle{\cal J}_{1}$ $\displaystyle=$ $\displaystyle{\cal J}(1,1)-{\cal
J}(0,1)\theta(\mu_{F}^{2}-{\bf q}_{T1}^{2})$ $\displaystyle-$
$\displaystyle{\cal J}(1,0)\theta(\mu_{F}^{2}-{\bf q}_{T2}^{2})+{\cal
J}(0,0)\theta(\mu_{F}^{2}-{\bf q}_{T1}^{2})\theta(\mu_{F}^{2}-{\bf
q}_{T2}^{2}),$ $\displaystyle{\cal J}_{2}$ $\displaystyle=$
$\displaystyle\left[{\cal J}(1,0)-{\cal J}(0,0)\theta(\mu_{F}^{2}-{\bf
q}_{T1}^{2})\right]\theta(\mu_{F}^{2}-{\bf q}_{T2}^{2}),$ $\displaystyle{\cal
J}_{3}$ $\displaystyle=$ $\displaystyle\left[{\cal J}(0,1)-{\cal
J}(0,0)\theta(\mu_{F}^{2}-{\bf q}_{T2}^{2})\right]\theta(\mu_{F}^{2}-{\bf
q}_{T1}^{2}).$

Using these expressions, the doubly-unintegrated contribution
$\check{\sigma}_{\rm 2UI-S}^{[m]}$ in this method is defined as eqn. (12) with
integrand ${\cal J}_{1}$ and without any low-${\bf q}_{T1,2}$ cuts. The
single-unintegrated contribution is obtained by substituting ${\cal J}_{2}$
into eqn. (12) and integrating out ${\bf q}_{T2}^{2}$ via eqn. (37):

$\displaystyle\check{\sigma}^{[m]}_{\rm 1UI-S}=\int\limits_{0}^{\infty}d{\bf
q}_{T1}^{2}\int\limits_{0}^{2\pi}\frac{d\phi}{2}\ \Bigg{[}$
$\displaystyle{\cal C}_{gg}\left(z\frac{M_{T1}^{2}}{M^{2}},{\bf
q}_{T1}^{2},\mu_{F}^{2},\mu_{R}^{2}\right)\frac{H({\bf
q}_{T1}^{2},0,\phi)}{M_{T1}^{4}}$ $\displaystyle-$ $\displaystyle{\cal
C}_{gg}\left(z,{\bf
q}_{T1}^{2},\mu_{F}^{2},\mu_{R}^{2}\right)\frac{H(0,0,\phi)}{M^{4}}\theta(\mu_{F}^{2}-{\bf
q}_{T1}^{2})\Bigg{]}.$ (78)

The LO subtracted resummed cross section for the $gg$ channel is calculated
as:

$\check{\sigma}_{gg}^{[m],\ {\rm HEF}}=\check{\sigma}_{\rm
2UI-S}^{[m]}+2\check{\sigma}^{[m]}_{\rm 1UI-S},$ (79)

while, for the quark-induced channels, one has:

$\displaystyle\hat{\sigma}_{gq}^{[m],\ {\rm HEF}}$ $\displaystyle=$
$\displaystyle\frac{C_{F}}{C_{A}}\left[\check{\sigma}_{\rm
2UI-S}^{[m]}+\check{\sigma}^{[m]}_{\rm 1UI-S}\right],$ (80)
$\displaystyle\hat{\sigma}_{q\bar{q}}^{[m],\ {\rm HEF}}$ $\displaystyle=$
$\displaystyle\left(\frac{C_{F}}{C_{A}}\right)^{2}\check{\sigma}_{\rm
2UI-S}^{[m]}.$ (81)

For the DLA resummation factor, the numerical results obtained in both methods
agree within the integration accuracy. However, the subtraction method allows
one to work with resummation factor given in the numerical form. The only
requirement is that it should satisfy the normalisation condition (37), which
is true for example for the resummation factor (24) with full LL($\ln(1/z)$)
anomalous dimension, obtained as numerical solution of eqn. (25), but with
$R(\gamma)=1$, as discussed in Appendix B, or for the Gaussian-smeared
resummation factor discussed in Appendix C.

##  Appendix B Effects of the anomalous dimension beyond LO

In this appendix, we study the effects of terms beyond LO in $\alpha_{s}$ in
the anomalous dimension (26) on the resummed cross section, while still
omitting the scheme-transformation factor $R(\gamma)$ from the eqn. (24). We
have obtained the numerical solution of eqn. (25), which in terms of the
variable $\rho=N/\hat{\alpha}_{s}$ has the form:

$\chi(\gamma_{gg}(\rho))=\rho.$ (82)

This equation has already been studied numerically in refs. Ellis:1995gv ;
Blumlein:1995eu . The perturbative branch of the solution has a cut-
discontinuity depicted in the figure 10, and inverse Mellin transform integral
(35) can be written as an integral over the contour in $\rho$-plane which
should avoid crossing this cut. Our numerical results for the real and
imaginary parts of the solution agree with results of ref. Blumlein:1995eu
and are also shown on figure 10.

Figure 10: Left panel: branch cut discontinuity of the function
$\gamma_{gg}(\rho)$ and the integration contour of the inverse Mellin
transform in the $\rho$-plane. Right panel: real and imaginary parts of
$\gamma_{gg}(\rho)$ as functions of $\mathop{\rm Im}\rho$ at fixed
$\mathop{\rm Re}\rho=1,0.5,0,-0.5,-1,-1.2,-1.42$. To be compared with the
plots in the figure 1 of the ref. Blumlein:1995eu .

Plots of the resummation factor resulting from the numerical inverse Mellin
transform are shown on figure 11. At moderately small values of $z$, it is
close to the DLA, but at smaller values of $z$ it deviates from DLA and the
most numerically important region of ${\bf q}_{T}^{2}\sim 1$ GeV2 features a
significant peak in the negative direction.

On figure 12, the $\sqrt{s}$-dependence of the HEF-resummed part of the total
cross section of production of the $Q\bar{Q}[{}^{1}S_{0}^{(1)}]$ state in the
$gg$ channel is plotted. The central lines as well as uncertainty bands
resulting from the variation of $\mu_{F}$ by a factor of 2 above and below the
default value $\mu_{F}=M$ are presented. Two approximations for the
resummation factor are compared on figure 12: the DLA, introduced in the
section 2.3, and the approximation described above in this appendix, which
includes effects beyond LO in $\gamma_{gg}(N,\alpha_{s})$. The LO CF central
curve and the $\mu_{F}$ uncertainty is also shown on figure 12 for comparison.
One can see that, while DLA significantly reduces the scale uncertainty in
comparison to the LO CF result, the inclusion of higher-order corrections to
$\gamma_{gg}(N,\alpha_{s})$ blows up the scale uncertainty at high energies so
much that, for $\sqrt{s}>4$ TeV (left panel of the figure 12) or $\sqrt{s}>7$
TeV (right panel of the figure 12), the resummed part of cross section can
become negative.

Figure 11: Resummation factor ${\cal C}_{gg}$ as a function of $|{\bf q}_{T}|$
for several values of $z$. The DLA result is plotted by the dashed curve and
numerical result with LL($\ln(1/z)$) anomalous dimension, obtained as a
solution of eqn. (25) is plotted by the solid curve.

Figure 12: Energy dependence of the HEF resummed part of the cross section (LO
CF included) in the $gg$ channel: the DLA result and it’s $\mu_{F}$-variation
are plotted as solid line and left-shaded band, dashed line and right-shaded
$\mu_{F}$-variation band correspond to the calculation using the full
numerical solution of eqn. (25) for $\gamma_{gg}(N,\alpha_{s})$. For
comparison, the LO CF cross section together with its $\mu_{F}$-variation is
shown by the dash-dotted line and the corresponding solid band. The gluon PDFs
used are NNPDF31_nlo_as_0118 NNPDF:2017mvq (left panel) and
NNPDF31sx_nlonllx_as_0118 Ball:2017otu (right panel). The LDME $\langle{\cal
O}[{}^{1}S_{0}^{(1)}]\rangle$ was set to 1 GeV3.

This catastrophic factorisation-scale uncertainty is a consequence of the
mismatch between the $\mu_{F}$-dependence of the resummation factor and one of
the collinear PDFs used. For the left plot on figure 12, we have employed the
central eigenset of the NNPDF31_nlo _as_0118 NNPDF:2017mvq  PDFs and, in the
right panel, the NLL($\ln(1/z)$)-resummed version of the NNPDF31 PDFs has been
used Ball:2017otu . Figure 12 serves as a numerical illustration of our
statement in section 2.3 that one should not increase the accuracy of the
resummation factors in HEF beyond DLA if collinear PDFs result from a fixed-
order DGLAP evolution, as in the left plot of figure 12. In the case of the
NLL($\ln(1/z)$)-resummed PDFs (right panel of figure 12), the addition of the
full LL($\ln(1/z)$) correction to $\gamma_{gg}$ is also not justified because
in this case the formal perturbative accuracy of the PDFs is greater that that
of the resummation factor and the NLL resummed version $\gamma_{gg}$ is
significantly Altarelli:1999vw ; Ball:2017otu  different from the full
LL($\ln(1/z)$)-resummed one, so even in this case we have a mismatch of the
$\mu_{F}$-dependence between the PDFs and the resummation factor.

##  Appendix C Higher-twist effects

It may be argued that the oscillatory behaviour of the resummation factors at
${\bf q}_{T}^{2}\ll 1$ GeV2, depicted on figure 11, is un-physical and that in
this domain the resummation factor is dominated by non-perturbative effects.
In this context, one would like to understand how much the total cross section
depends on the detailed behaviour of the resummation factor at small values of
${\bf q}_{T}^{2}$. As a toy model to study this question, we have performed a
convolution of the DLA resummation factor (34) with Gaussian of the width
$\sigma_{T}$ in transverse momentum which can be done analytically with the
Mellin-space result (34):

$\displaystyle{\cal C}_{gg}^{{\rm DLA},\sigma}(N,{\bf
q}_{T}^{2},\mu_{F},\mu_{R})$ $\displaystyle=$
$\displaystyle\int\frac{d^{2}{\bf
k}_{T}}{\pi\sigma_{T}^{2}}\exp\left[-\frac{{\bf
k}_{T}^{2}}{\sigma_{T}^{2}}\right]{\cal C}^{\rm DLA}_{gg}(N,({\bf q}_{T}+{\bf
k}_{T})^{2},\mu_{F},\mu_{R})$ (83) $\displaystyle=$
$\displaystyle\frac{1}{\sigma_{T}^{2}}\left(\frac{\sigma_{T}^{2}}{\mu_{F}^{2}}\right)^{\gamma_{N}}e^{-\frac{{\bf
q}_{T}^{2}}{\sigma_{T}^{2}}}{}_{1}F_{1}\left(\gamma_{N},1,\frac{{\bf
q}_{T}^{2}}{\sigma_{T}^{2}}\right),$

and then converted it to the $z$-space numerically. Eqn. (83) still satisfies
the normalisation condition (37) up to corrections suppressed as
$e^{-\mu_{F}^{2}/\sigma_{T}^{2}}$, so the method from the section A.2 can be
used to calculate the cross section if $\sigma_{T}\lesssim 1$ GeV. Several
plots of the $|{\bf q}_{T}|$-dependence of the smeared resummation factor
(83), transformed numerically to the $z$-space are presented in the left panel
of figure 13.

In the right panel of figure 13, the ratio of the matched cross section
obtained with the Gaussian-smeared resummation factor to the cross section
without ${\bf k}_{T}$-smearing is plotted as a function of $\sqrt{s}$ for
$\sigma_{T}=0.5$ and $1$ GeV. The results with $\sigma_{T}=0.5$ and $1$ GeV
differ from the matched cross section calculated without a Gaussian by at most
ten percent for $\mu_{F}=0.5M$ and $2M$, while the results at the default
scale are even closer.

From the detailed consideration of the small-${\bf q}_{T}$ subtraction method,
described in section A.2, one can see that the effect of the ${\bf
k}_{T}$-smearing on the total cross section is suppressed as ${\cal
O}(\sigma_{T}^{2}/M^{2})\sim O(\Lambda_{\rm QCD}^{2}/M^{2})$ if
$\sigma_{T}\sim\Lambda_{QCD}$. In other words, these effect follow a
subleading power in our hard scale $M$; it is a higher-twist effect.

The width of the “intrinsic” ${\bf k}_{T}$-distribution is expected to be
increasing with energy, due to increasing saturation scale (see e.g. ref.
Albacete:2014fwa  and references therein), so that, at high enough energies,
the twist expansion will completely break down due to the effect discussed
above. However, as one can see from figure 13 this effects are unlikely to be
the main source of uncertainty in our calculation at LHC energies and even
beyond.

Figure 13: Left panel: the ${\bf q}_{T}$-dependence of the Gaussian smeared
DLA resummation factor for two values of $\sigma_{T}$ (dashed curves). The DLA
resummation factor without smearing (36) is shown by solid curve for
comparison. Right panel: the energy-dependence of the ratio of the matched
cross section obtained in the DLA with Gaussian smearing with $\sigma_{T}=0.5$
and 1 GeV to the result without smearing and its $\mu_{F}$-scale variation.

## References

  * (1) E. Chapon et al., _Perspectives for quarkonium studies at the high-luminosity LHC_ , 2012.14161.
  * (2) A. Arbuzov et al., _On the physics potential to study the gluon content of proton and deuteron at NICA SPD_ , _Prog. Part. Nucl. Phys._ 119 (2021) 103858 [2011.15005].
  * (3) J.-P. Lansberg, _New Observables in Inclusive Production of Quarkonia_ , _Phys. Rept._ 889 (2020) 1 [1903.09185].
  * (4) N. Brambilla et al., _QCD and Strongly Coupled Gauge Theories: Challenges and Perspectives_ , _Eur. Phys. J. C_ 74 (2014) 2981 [1404.3723].
  * (5) N. Brambilla et al., _Heavy Quarkonium: Progress, Puzzles, and Opportunities_ , _Eur. Phys. J. C_ 71 (2011) 1534 [1010.5827].
  * (6) G.T. Bodwin, E. Braaten and G.P. Lepage, _Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium_ , _Phys. Rev. D_ 51 (1995) 1125 [hep-ph/9407339].
  * (7) M. Butenschoen and B.A. Kniehl, _Reconciling $J/\psi$ production at HERA, RHIC, Tevatron, and LHC with NRQCD factorization at next-to-leading order_, _Phys. Rev. Lett._ 106 (2011) 022003 [1009.5662].
  * (8) M. Butenschoen and B.A. Kniehl, _World data of J/psi production consolidate NRQCD factorization at NLO_ , _Phys. Rev. D_ 84 (2011) 051501 [1105.0820].
  * (9) M. Butenschoen and B.A. Kniehl, _J/psi polarization at Tevatron and LHC: Nonrelativistic-QCD factorization at the crossroads_ , _Phys. Rev. Lett._ 108 (2012) 172002 [1201.1872].
  * (10) M. Butenschoen and B.A. Kniehl, _Next-to-leading-order tests of NRQCD factorization with $J/\psi$ yield and polarization_, _Mod. Phys. Lett. A_ 28 (2013) 1350027 [1212.2037].
  * (11) K.-T. Chao, Y.-Q. Ma, H.-S. Shao, K. Wang and Y.-J. Zhang, _$J/\psi$ Polarization at Hadron Colliders in Nonrelativistic QCD_, _Phys. Rev. Lett._ 108 (2012) 242004 [1201.2675].
  * (12) B. Gong, L.-P. Wan, J.-X. Wang and H.-F. Zhang, _Polarization for Prompt J/ $\psi$ and $\psi$(2s) Production at the Tevatron and LHC_, _Phys. Rev. Lett._ 110 (2013) 042002 [1205.6682].
  * (13) G.T. Bodwin, H.S. Chung, U.-R. Kim and J. Lee, _Fragmentation contributions to $J/\psi$ production at the Tevatron and the LHC_, _Phys. Rev. Lett._ 113 (2014) 022001 [1403.3612].
  * (14) B. Gong, J.-X. Wang and H.-F. Zhang, _QCD corrections to $\Upsilon$ production via color-octet states at the Tevatron and LHC_, _Phys. Rev. D_ 83 (2011) 114021 [1009.3839].
  * (15) K. Wang, Y.-Q. Ma and K.-T. Chao, _$\Upsilon(1S)$ prompt production at the Tevatron and LHC in nonrelativistic QCD_, _Phys. Rev. D_ 85 (2012) 114003 [1202.6012].
  * (16) B. Gong, L.-P. Wan, J.-X. Wang and H.-F. Zhang, _Complete next-to-leading-order study on the yield and polarization of $\Upsilon(1S,2S,3S)$ at the Tevatron and LHC_, _Phys. Rev. Lett._ 112 (2014) 032001 [1305.0748].
  * (17) M. Butenschoen, Z.-G. He and B.A. Kniehl, _$\eta_{c}$ production at the LHC challenges nonrelativistic-QCD factorization_, _Phys. Rev. Lett._ 114 (2015) 092004 [1411.5287].
  * (18) V.D. Barger, W.-Y. Keung and R.J.N. Phillips, _On psi and Upsilon Production via Gluons_ , _Phys. Lett. B_ 91 (1980) 253.
  * (19) V.D. Barger, W.-Y. Keung and R.J.N. Phillips, _Hadroproduction of $\psi$ and $\Upsilon$_, _Z. Phys. C_ 6 (1980) 169.
  * (20) R. Gavai, D. Kharzeev, H. Satz, G.A. Schuler, K. Sridhar and R. Vogt, _Quarkonium production in hadronic collisions_ , _Int. J. Mod. Phys. A_ 10 (1995) 3043 [hep-ph/9502270].
  * (21) Y.-Q. Ma and R. Vogt, _Quarkonium Production in an Improved Color Evaporation Model_ , _Phys. Rev. D_ 94 (2016) 114029 [1609.06042].
  * (22) J.-P. Lansberg, H.-S. Shao, N. Yamanaka, Y.-J. Zhang and C. Noûs, _Complete NLO QCD study of single- and double-quarkonium hadroproduction in the colour-evaporation model at the Tevatron and the LHC_ , _Phys. Lett. B_ 807 (2020) 135559 [2004.14345].
  * (23) N. Brambilla, H.S. Chung and A. Vairo, _Inclusive production of heavy quarkonia in pNRQCD_ , _JHEP_ 09 (2021) 032 [2106.09417].
  * (24) Y.-Q. Ma and K.-T. Chao, _New factorization theory for heavy quarkonium production and decay_ , _Phys. Rev. D_ 100 (2019) 094007 [1703.08402].
  * (25) R. Li, Y. Feng and Y.-Q. Ma, _Exclusive quarkonium production or decay in soft gluon factorization_ , _JHEP_ 05 (2020) 009 [1911.05886].
  * (26) A.-P. Chen and Y.-Q. Ma, _Theory for quarkonium: from NRQCD factorization to soft gluon factorization_ , _Chin. Phys. C_ 45 (2021) 013118 [2005.08786].
  * (27) Z.-B. Kang, J.-W. Qiu and G. Sterman, _Heavy quarkonium production and polarization_ , _Phys. Rev. Lett._ 108 (2012) 102002 [1109.1520].
  * (28) J.H. Kuhn and E. Mirkes, _QCD corrections to toponium production at hadron colliders_ , _Phys. Rev. D_ 48 (1993) 179 [hep-ph/9301204].
  * (29) A. Petrelli, M. Cacciari, M. Greco, F. Maltoni and M.L. Mangano, _NLO production and decay of quarkonium_ , _Nucl. Phys. B_ 514 (1998) 245 [hep-ph/9707223].
  * (30) J.-P. Lansberg and M.A. Ozcelik, _Curing the unphysical behaviour of NLO quarkonium production at the LHC and its relevance to constrain the gluon PDF at low scales_ , _Eur. Phys. J. C_ 81 (2021) 497 [2012.00702].
  * (31) G.A. Schuler, _Quarkonium production and decays_ , 1994.
  * (32) M.L. Mangano and A. Petrelli, _NLO quarkonium production in hadronic collisions_ , _Int. J. Mod. Phys. A_ 12 (1997) 3887 [hep-ph/9610364].
  * (33) Y. Feng, J.-P. Lansberg and J.-X. Wang, _Energy dependence of direct-quarkonium production in $pp$ collisions from fixed-target to LHC energies: complete one-loop analysis_, _Eur. Phys. J. C_ 75 (2015) 313 [1504.00317].
  * (34) V.N. Gribov and L.N. Lipatov, _$e^{+}e^{-}$ -annihilation and deep-inelastic $ep$-scattering in perturbation theory_, _Sov. J. Nucl. Phys._ 15 (1972) 438.
  * (35) Y.L. Dokshitzer, _Calculation of structure functions of deep-inelastic scattering and $e^{+}e^{-}$-annihilation in perturbation theory of quantum chromodynamics_, _Sov. Phys. JETP_ 46 (1977) 641.
  * (36) G. Altarelli and G. Parisi, _Asymptotic freedom in parton language_ , _Nucl. Phys._ B126 (1977) 298.
  * (37) R.K. Ellis and D.A. Ross, _The Coupling of the QCD Pomeron in Various Semihard Processes_ , _Nucl. Phys. B_ 345 (1990) 79.
  * (38) A.C. Serri, Y. Feng, C. Flore, J.-P. Lansberg, M.A. Ozcelik, H.-S. Shao et al., _Revisiting NLO QCD corrections to total inclusive J/psi and Upsilon photoproduction cross sections in lepton-proton collisions_ , 2112.05060.
  * (39) E.A. Kuraev, L.N. Lipatov and V.S. Fadin, _Multi - Reggeon processes in the Yang-Mills theory_ , _Sov. Phys. JETP_ 44 (1976) 443.
  * (40) E.A. Kuraev, L.N. Lipatov and V.S. Fadin, _The Pomeranchuk singularity in non-Abelian gauge theories_ , _Sov. Phys. JETP_ 45 (1977) 199.
  * (41) Y.Y. Balitsky and L.N. Lipatov, _The Pomeranchuk singularity in Quantum Chromodynamics_ , _Sov. J. Nucl. Phys._ 28 (1978) 822.
  * (42) S. Catani, M. Ciafaloni and F. Hautmann, _GLUON CONTRIBUTIONS TO SMALL x HEAVY FLAVOR PRODUCTION_ , _Phys. Lett. B_ 242 (1990) 97.
  * (43) S. Catani, M. Ciafaloni and F. Hautmann, _High-energy factorization and small x heavy flavor production_ , _Nucl. Phys. B_ 366 (1991) 135.
  * (44) J.C. Collins and R.K. Ellis, _Heavy quark production in very high-energy hadron collisions_ , _Nucl. Phys._ B360 (1991) 3.
  * (45) S. Catani and F. Hautmann, _High-energy factorization and small x deep inelastic scattering beyond leading order_ , _Nucl. Phys._ B427 (1994) 475 [hep-ph/9405388].
  * (46) S. Catani, M. Ciafaloni and F. Hautmann, _Leptoproduction of heavy flavor at high energies_ , _Nucl. Phys. B Proc. Suppl._ 29 (1992) 182.
  * (47) F. Hautmann, _Heavy top limit and double logarithmic contributions to Higgs production at m(H)**2 / s much less than 1_ , _Phys. Lett._ B535 (2002) 159 [hep-ph/0203140].
  * (48) R.V. Harlander, H. Mantler, S. Marzani and K.J. Ozeren, _Higgs production in gluon fusion at next-to-next-to-leading order QCD for finite top mass_ , _Eur. Phys. J. C_ 66 (2010) 359 [0912.2104].
  * (49) S. Marzani and R.D. Ball, _High Energy Resummation of Drell-Yan Processes_ , _Nucl. Phys. B_ 814 (2009) 246 [0812.3602].
  * (50) G. Diana, _High-energy resummation in direct photon production_ , _Nucl. Phys. B_ 824 (2010) 154 [0906.4159].
  * (51) 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._ C78 (2018) 321 [1710.05935].
  * (52) xFitter Developers’ Team collaboration, _Impact of low- $x$ resummation on QCD analysis of HERA data_, _Eur. Phys. J._ C78 (2018) 621 [1802.00064].
  * (53) M.G. Echevarria, T. Kasemets, J.-P. Lansberg, C. Pisano and A. Signori, _Matching factorization theorems with an inverse-error weighting_ , _Phys. Lett. B_ 781 (2018) 161 [1801.01480].
  * (54) Y.-Q. Ma and R. Venugopalan, _Comprehensive Description of J/ $\psi$ Production in Proton-Proton Collisions at Collider Energies_, _Phys. Rev. Lett._ 113 (2014) 192301 [1408.4075].
  * (55) Z.-B. Kang, Y.-Q. Ma and R. Venugopalan, _Quarkonium production in high energy proton-nucleus collisions: CGC meets NRQCD_ , _JHEP_ 01 (2014) 056 [1309.7337].
  * (56) I. Balitsky, _Operator expansion for high-energy scattering_ , _Nucl. Phys. B_ 463 (1996) 99 [hep-ph/9509348].
  * (57) Y.V. Kovchegov, _Small x F(2) structure function of a nucleus including multiple pomeron exchanges_ , _Phys. Rev. D_ 60 (1999) 034008 [hep-ph/9901281].
  * (58) I. Balitsky, _Factorization for high-energy scattering_ , _Phys. Rev. Lett._ 81 (1998) 2024 [hep-ph/9807434].
  * (59) J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, _The BFKL equation from the Wilson renormalization group_ , _Nucl. Phys. B_ 504 (1997) 415 [hep-ph/9701284].
  * (60) J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, _The Wilson renormalization group for low x physics: Towards the high density regime_ , _Phys. Rev. D_ 59 (1998) 014014 [hep-ph/9706377].
  * (61) J. Jalilian-Marian, A. Kovner and H. Weigert, _The Wilson renormalization group for low x physics: Gluon evolution at finite parton density_ , _Phys. Rev. D_ 59 (1998) 014015 [hep-ph/9709432].
  * (62) E. Iancu, A. Leonidov and L.D. McLerran, _The Renormalization group equation for the color glass condensate_ , _Phys. Lett. B_ 510 (2001) 133 [hep-ph/0102009].
  * (63) E. Iancu, A. Leonidov and L.D. McLerran, _Nonlinear gluon evolution in the color glass condensate. 1._ , _Nucl. Phys. A_ 692 (2001) 583 [hep-ph/0011241].
  * (64) B.L. Ioffe, V.S. Fadin and L.N. Lipatov, _Quantum Chromodynamics Perturbative and Nonperturbative Aspects_ , Cambridge University Press, Cambridge (2010).
  * (65) Y.V. Kovchegov and E. Levin, _Quantum Chromodynamics at High Energy_ , Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology, Cambridge University Press (2012), 10.1017/CBO9781139022187.
  * (66) V. Del Duca, _An introduction to the perturbative QCD pomeron and to jet physics at large rapidities_ , _Scientifica Acta_ 10 (1995) 91 [hep-ph/9503226].
  * (67) J.C. Collins, _Foundations of perturbative QCD_ , Cambridge University Press, Cambridge (2011).
  * (68) M. Nefedov, _Sudakov resummation from the BFKL evolution_ , _Phys. Rev. D_ 104 (2021) 054039 [2105.13915].
  * (69) M. Hentschinski, _Transverse momentum dependent gluon distribution within high energy factorization at next-to-leading order_ , _Phys. Rev. D_ 104 (2021) 054014 [2107.06203].
  * (70) L.N. Lipatov, _Gauge invariant effective action for high-energy processes in QCD_ , _Nucl. Phys._ B452 (1995) 369.
  * (71) B.A. Kniehl, V.A. Saleev and D.V. Vasin, _Bottomonium production in the Regge limit of QCD_ , _Phys. Rev._ D74 (2006) 014024 [hep-ph/0607254].
  * (72) Z.-G. He, B.A. Kniehl, M.A. Nefedov and V.A. Saleev, _Double Prompt $J/\psi$ Hadroproduction in the Parton Reggeization Approach with High-Energy Resummation_, _Phys. Rev. Lett._ 123 (2019) 162002 [1906.08979].
  * (73) P. Hagler, R. Kirschner, A. Schafer, L. Szymanowski and O. Teryaev, _Heavy quark production as sensitive test for an improved description of high-energy hadron collisions_ , _Phys. Rev. D_ 62 (2000) 071502 [hep-ph/0002077].
  * (74) P. Hagler, R. Kirschner, A. Schafer, L. Szymanowski and O.V. Teryaev, _Towards a solution of the charmonium production controversy: $k^{-}$ perpendicular factorization versus color octet mechanism_, _Phys. Rev. Lett._ 86 (2001) 1446 [hep-ph/0004263].
  * (75) B.A. Kniehl, D.V. Vasin and V.A. Saleev, _Charmonium production at high energy in the $k_{T}$ -factorization approach_, _Phys. Rev. D_ 73 (2006) 074022 [hep-ph/0602179].
  * (76) T. Jaroszewicz, _Gluonic Regge Singularities and Anomalous Dimensions in QCD_ , _Phys. Lett._ 116B (1982) 291.
  * (77) R. Kirschner and M. Segond, _Small x resummation in collinear factorisation_ , _Eur. Phys. J. C_ 68 (2010) 425 [0910.5443].
  * (78) A. Vogt, S. Moch and J.A.M. Vermaseren, _The Three-loop splitting functions in QCD: The Singlet case_ , _Nucl. Phys. B_ 691 (2004) 129 [hep-ph/0404111].
  * (79) T. Lukowski, A. Rej and V.N. Velizhanin, _Five-Loop Anomalous Dimension of Twist-Two Operators_ , _Nucl. Phys. B_ 831 (2010) 105 [0912.1624].
  * (80) B.A. Kniehl, V.A. Saleev, A.V. Shipilova and E.V. Yatsenko, _Single jet and prompt-photon inclusive production with multi-Regge kinematics: From Tevatron to LHC_ , _Phys. Rev. D_ 84 (2011) 074017 [1107.1462].
  * (81) J. Blumlein, _On the k(T) dependent gluon density of the proton_ , in _Deep inelastic scattering and QCD. Proceedings, Workshop, Paris, France, April 24-28, 1995_ , pp. 265–268, 1995, http://www-library.desy.de/cgi-bin/showprep.pl?desy95-121 [hep-ph/9506403].
  * (82) 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].
  * (83) 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].
  * (84) NNPDF collaboration, _Parton distributions from high-precision collider data_ , _Eur. Phys. J. C_ 77 (2017) 663 [1706.00428].
  * (85) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _Eur. Phys. J._ C75 (2015) 132 [1412.7420].
  * (86) R.K. Ellis, F. Hautmann and B.R. Webber, _QCD scaling violation at small x_ , _Phys. Lett. B_ 348 (1995) 582 [hep-ph/9501307].
  * (87) G. Altarelli, R.D. Ball and S. Forte, _Resummation of singlet parton evolution at small x_ , _Nucl. Phys._ B575 (2000) 313 [hep-ph/9911273].
  * (88) J.L. Albacete and C. Marquet, _Gluon saturation and initial conditions for relativistic heavy ion collisions_ , _Prog. Part. Nucl. Phys._ 76 (2014) 1 [1401.4866].