11institutetext: School of Computing, National University of Singapore
22institutetext: Institute of Data Science, National University of Singapore
33institutetext: Institute for Infocomm Research, A*STAR, Singapore
33email: adam.goodge@u.nus.edu,{dcsbhk<EMAIL_ADDRESS><EMAIL_ADDRESS>
# ARES: Locally Adaptive Reconstruction-based Anomaly Scoring

Adam Goodge 1133 Bryan Hooi 1122 See Kiong Ng 1122 Wee Siong Ng 33

###### Abstract

How can we detect anomalies: that is, samples that significantly differ from a
given set of high-dimensional data, such as images or sensor data? This is a
practical problem with numerous applications and is also relevant to the goal
of making learning algorithms more robust to unexpected inputs. Autoencoders
are a popular approach, partly due to their simplicity and their ability to
perform dimension reduction. However, the anomaly scoring function is not
adaptive to the natural variation in reconstruction error across the range of
normal samples, which hinders their ability to detect real anomalies. In this
paper, we empirically demonstrate the importance of local adaptivity for
anomaly scoring in experiments with real data. We then propose our novel
Adaptive Reconstruction Error-based Scoring approach, which adapts its scoring
based on the local behaviour of reconstruction error over the latent space. We
show that this improves anomaly detection performance over relevant baselines
in a wide variety of benchmark datasets.

###### Keywords:

Anomaly detection machine learning unsupervised learning.

##  1 Introduction

The detection of anomalous data is important in a wide variety of
applications, such as detecting fraudulent financial transactions or malignant
tumours. Recently, deep learning methods have enabled significant improvements
in performance on ever larger and higher-dimensional datasets. Despite this,
anomaly detection remains a challenging task, most notably due to the
difficulty of obtaining accurate labels of anomalous data. For this reason,
supervised classification methods are unsuitable for anomaly detection.
Instead, unsupervised methods are used to learn the distribution of an all-
normal training set. Anomalies are then detected amongst unseen samples by
measuring their closeness to the normal-data distribution.

Autoencoders are an extremely popular approach to learn the behaviour of
normal data. The reconstruction error of a sample is directly used as its
anomaly score; anomalies are assumed to have a higher reconstruction error
than normal samples due to their difference in distribution. However, this
anomaly score fails to account for the fact that reconstruction error can vary
greatly even amongst different types of normal samples. For example, consider
a sensor system in a factory with various activities on weekdays but zero
activity on weekends. The weekday samples are diverse and complex, therefore
we could expect the reconstruction error of these samples to vary greatly.
Meanwhile, even a small amount of activity in a weekend sample would be
anomalous, even if the effect on the reconstruction error is minimal. The
anomaly detector would likely detect false positives (high reconstruction
error) among weekday samples and false negatives among weekend samples.
Although encoded within the input data attributes, the context of each
individual sample (i.e. day of the week) is neglected at the detection stage
by the standard scoring approach. Instead, all samples are assessed according
to the same error threshold or standard. This invokes the implicit assumption
that the reconstruction errors of all samples are identically distributed,
regardless of any individual characteristics and context which could
potentially influence the reconstruction error significantly.

In this work, we aim to address this problem by proposing Adaptive
Reconstruction Error-based Scoring (ARES). Our locally-adaptive scoring method
is able to automatically account for any contextual information which affects
the reconstruction error, resulting in more accurate anomaly detection. We use
a flexible, neighbourhood-based approach to define the context, based on the
location of its latent representation learnt by the model.

Our scoring approach is applied at test time, so it can be retrofitted to pre-
trained models of any size and architecture or used to complement existing
anomaly detection techniques. Our score is simple and efficient to compute,
requiring very little additional computational time, and the code is available
online111https://github.com/agoodge/ARES. In summary, we make the following
contributions:

  1. 1.

We empirically show with real, multi-class data the variation in
reconstruction error among different samples in the normal set and why this
justifies the need for local adaptivity in anomaly scoring.

  2. 2.

We propose a novel anomaly scoring method which adapts to this variation by
evaluating anomaly status based on the local context of a given sample in the
latent space.

  3. 3.

We evaluate our method against a wide range of baselines with various
benchmark datasets. We also study the effect of different components and
formulations of our scoring method in an ablation study.

##  2 Background

###  2.1 Autoencoders

Autoencoders are neural networks that output a reconstruction of the input
data [26]. They are comprised of two components: an encoder and decoder. The
encoder compresses data from the input level into lower-dimensional latent
representations through its hidden layers. The encoder output is typically
known as the bottleneck, from which the reconstruction of the original data is
found through the decoder hidden layers. The network is trained to minimise
the reconstruction error over the training set:

$\displaystyle\min_{\theta,\phi}\|\mathbf{x}-(f_{\theta}\circ
g_{\phi})(\mathbf{x})\|_{2}^{2},$ (1)

where $g_{\phi}$ is the encoder, $f_{\theta}$ the decoder. The assumption is
that data lies on a lower-dimensional manifold within the high-dimensional
input space. The autoencoder learns to reconstruct data on this manifold by
performing dimension reduction. By training the model on only normal data, the
reconstruction error of a normal sample should be low as it close to the
learnt manifold on which it has been reconstructed by the autoencoder, whereas
anomalies are far away and are reconstructed with higher error.

###  2.2 Local Outlier Factor

The Local Outlier Factor (LOF) method is a neighbourhood-based approach to
anomaly detection; it measures the density of a given sample relative to its
$k$ nearest neighbours’ [8]. Anomalies are assumed to be in sparse regions far
away from the one or more high-density clusters of normal data. A lower
density therefore suggests the sample is anomalous. The original method uses
the ‘reachability distance’; defined for a point $A$ from point $B$ as:

$\max\\{k\text{-distance}(B),d(A,B)\\}$ (2)

where $d(\cdot,\cdot)$ is a chosen distance metric and $k$-distance$(B)$ is
the distance of $B$ to its $k^{th}$ nearest neighbour. The local reachability
density and subsequently the local outlier factor of $A$ based on its set of
neighbours $\mathsf{N}_{k}(A)$, is:

$\text{lrd}_{k}(A):=\left(\frac{\sum_{B\in\mathsf{N}_{k}(A)}\text{reachability-
distance}_{k}(A,B)}{|\mathsf{N}_{k}(A)|}\right)^{-1}.$ (3)

$\text{LOF}_{k}(A)=\frac{\sum_{B\in\mathsf{N}_{k}(A)}\text{lrd}_{k}(B)}{|\mathsf{N}_{k}(A)|\cdot\text{lrd}_{k}(A)},$
(4)

##  3 Related Work

Anomalies are often assumed to occupy sparse regions far away from high-
density clusters of normal data. As such, a great variety of methods detect
anomalies by measuring the distances to their nearest neighbours, most notably
KNN [35], which directly uses the distance of a point to its $k^{th}$ nearest
neighbour as the anomaly score. LOF [8] measures the density of a point
relative to the density of points in its local neighbourhood, as seen in
Section 2. This is beneficial as a given density may be anomalous in one
region but normal in another region. This local view accommodates the natural
variation in density and therefore allows for the detection of more meaningful
anomalies. More recently, deep models such as graph neural networks have been
used to learn anomaly scores from neighbour-distances in a data-driven way
[15]. Naturally, this relies on an effective distance metric. This is often
itself a difficult problem; even the most established metrics have been shown
to lose significance in high-dimensional spaces [6] due to the ‘curse of
dimensionality’.

Reconstruction-based methods rely on deep models, such as autoencoders, to
learn to reconstruct samples from the normal set accurately. Samples are then
flagged as anomalous if their reconstruction error higher than some pre-
defined threshold, as it is assumed that the model will reconstruct samples
outside of the normal set with a higher reconstruction error [26, 13, 4, 9,
14].

More recent works have used autoencoders or other models are deep feature
extractors, with the learnt features then used in another anomaly detection
module downstream, such as Gaussian mixture models [7], DBSCAN [3], KNN [5]
and auto-regressive models [1]. [11] calibrates different types of
autoencoders against the effects of varying hardness within normal samples.

Other methods measure the distance of a sample to a normal set-enclosing
hypersphere [30, 25]. or its likelihood under a learnt model [12, 24, 34].
Generative adversarial networks have also been proposed for anomaly detection,
mostly relying either on the discriminator score or the accuracy of the
generator for a given sample to determine anomalousness [2, 33].

In all of these methods, the normal set is often restricted to a subset of the
available data; e.g. just one class from a multi-class dataset in experiments.
In practice, normal data could belong to multiple classes or modes of
behaviour which all need to be modelled adequately. As such, developing
anomaly detection methods that can model a diverse range of normal behaviours
and adapt their scoring appropriately are important.

##  4 Methodology

###  4.1 Problem Definition

We assume to have $m$ normal training samples
$\mathbf{x}_{1}^{\text{(train)}},\cdots,\mathbf{x}_{m}^{\text{(train)}}\in\mathbb{R}^{d}$
and $n$ testing samples,
$\mathbf{x}_{1}^{\text{(test)}},\cdots,\mathbf{x}_{n}^{\text{(test)}}\in\mathbb{R}^{d}$,
each of which may be normal or anomalous. For each test sample $\mathbf{x}$,
our algorithm should indicate how anomalous it is through computing its
anomaly score $s(\mathbf{x})$. The goal is for anomalies to be given higher
anomaly scores than normal points. In this work, the fundamental question is:

Given an autoencoder with encoder $g_{\phi}$ and decoder $f_{\theta}$, how can
we use the latent encoding $\mathbf{z}=g_{\phi}(\mathbf{x}$) and
reconstruction $\hat{\mathbf{x}}=(f_{\theta}\circ g_{\phi})(\mathbf{x})$ of a
sample $\mathbf{x}$ to score its anomalousness?

Our approach is not limited to standard autoencoders, and can be applied to
other models that involve a latent encoding $\mathbf{z}$ and a reconstruction
$\hat{\mathbf{x}}$.

In practice, the anomaly score is compared with a user-defined anomaly
threshold; samples which exceed this threshold are flagged as anomalies.
Different approaches can be employed to set this threshold, such as extreme
value theory [28]. Our focus is instead on the approach to anomaly scoring
itself, which allows for any choice of thresholding scheme to be used
alongside it.

###  4.2 Statistical Interpretation of Reconstruction-based Anomaly Detection

In the standard approach, with the residual as
$\bm{\varepsilon}:=\mathbf{x}-\hat{\mathbf{x}}$, the anomaly score is:

$R(\bm{\varepsilon}):=\|\bm{\varepsilon}\|_{2}^{2}=\|\mathbf{x}-\hat{\mathbf{x}}\|_{2}^{2}.$
(5)

This can be seen as a negative log-likelihood-based score that assumes
$\bm{\varepsilon}$ follows a Gaussian distribution with zero mean and unit
variance: $\bm{\varepsilon}\sim\mathcal{N}(0,\mathbf{I})$:

$-\log
P(\bm{\varepsilon})=\frac{d}{2}\log(2\pi)+\frac{R(\bm{\varepsilon})}{2}$ (6)

The negative log-likelihood is an intuitive anomaly score because anomalous
samples should have lower likelihood, thus higher negative log-likelihood.

The key conclusion, ignoring the constant additive and multiplicative factors,
is that $R(\bm{\varepsilon})$ can be equivalently seen as a negative log-
likelihood-based score, based on a model with the implicit assumption of
constant mean as well as _constant variance residuals across all samples_ ,
known in the statistical literature as homoscedasticity. Most crucially, this
assumption applies to all $\mathbf{x}$ regardless of its individual latent
representation, $\mathbf{z}$, which implies that the reconstruction error of
$\mathbf{x}$ does not depend on the location of $\mathbf{z}$ in the latent
space. We question the validity of this assumption in real datasets.

In our earlier example of a factory sensor system, the variance of a sample
$\mathbf{x}$ depends greatly on its characteristics (i.e. day of the week),
with high variance during weekdays and low variance on weekends. This is
instead a prime example of heteroscedasticity. Furthermore, since the latent
representation $\mathbf{z}$ typically encodes these important characteristics
of $\mathbf{x}$, there is likely to be a clear dependence between $\mathbf{z}$
and the reconstruction error $R(\bm{\varepsilon})$, which will be examined in
Section 4.3. This relationship could be exploited to improve the modelling of
reconstruction errors and thus detection accuracy by adapting the anomaly
score to this relationship at the sample-level.

###  4.3 Motivation & Empirical Results

Figure 1: t-SNE plot of the latent encodings with colour determined by
reconstruction error of the associated sample.

In this section, we train an autoencoder on all 10 classes of MNIST to
empirically examine the variation in the reconstruction error of real data to
scrutinise the homoscedastic assumption. In doing so, we exemplify the
shortcomings of the standard approach which implicitly assumes all
reconstruction errors from a fixed Gaussian distribution. The autoencoder
architecture and training procedure is the same as those described in Section
5. We make the following observations:

  1. 1.

Inter-class variation: In Figure 1, the t-SNE [21] projections of the latent
representations of training points are coloured according to their
reconstruction error. We see that the autoencoder learns to separate each
class approximately into distinct clusters. In this context, the class label
can be seen as a variable characteristic between different samples within the
normal set (similar to the weekday vs. weekend example).

Figure 2: Probability density function of reconstruction errors associated
with different classes of training samples, estimated via kernel density
estimation.

  2. 2.

Intra-class variation Also in Figure 1, we can see that there is significant
variation even within any single cluster. In other words, there are noticeably
distinct regions of high (and low) reconstruction errors within each
individual cluster. This shows that there are additional characteristics that
influence whether a given normal sample has a high or low reconstruction error
besides its class label.

We observe that there is significant variation in reconstruction error between
the different clusters. Most notably, samples in the leftmost cluster
(corresponding to class 1) has significantly lower reconstruction errors than
most others. Indeed, as shown in Figure 2, the distribution of reconstruction
errors associated with class 1 samples is very different to those of the other
classes.

This shows that there is significant variation in reconstruction errors
between classes, and that it is inappropriate to assume that the
reconstruction errors of all samples can be modelled by a single, fixed
Gaussian distribution. For example, a reconstruction error of $0.06$ would be
high for a class 1 or class 9 sample, but low for a class 2 or class 8 sample.

Figure 3: Average reconstruction error of the training set nearest neighbours
of a point versus its own reconstruction error for normal test samples.

  3. 3.

Neighbourhood correlation: In Figure 3, we plot the reconstruction error of
test samples against the average reconstruction error of its nearest training
set neighbours in the latent space (neighbourhood error). We see that the
reconstruction error of a test point increases as its neighbourhood error
increases. Furthermore, the variance in test errors increases for larger
neighbourhood errors: a clear heteroscedastic relationship. This information
is useful in determining anomalousness: a test error of $0.1$ would be
anomalously high if its neighbourhood error is $0.02$, but normal if it is
$0.06$.

Given these observations, we propose that anomalies could be more accurately
detected by incorporating contextual information into the anomaly scoring
beyond reconstruction error alone. Furthermore, analysing the neighbourhood of
a given sample provides this contextual information to help determine its
anomalousness. In section 4.4, we propose our Adaptive Reconstruction Error-
based Scoring method to achieve this.

###  4.4 Adaptive Reconstruction Error-based Scoring

In Section 4.2, we saw that the standard approach assumes that all residuals
come from a fixed Gaussian with constant mean and unit variance:
$\bm{\varepsilon}\sim\mathcal{N}(0,\mathbf{I})$. In Section 4.3, we saw that
this assumption is inappropriate. In this section, we detail ARES, a novel
anomaly scoring methodology which aims to address this flaw by adapting the
scoring for each samples local context in the latent space.

The normal level of reconstruction error varies for samples in different
regions of the latent space, meaning the latent encoding of a sample holds
important information regarding anomalousness. As such, ARES is inspired by
the joint-likelihood of a samples residual $\bm{\varepsilon}$ with its latent
encoding $\mathbf{z}$, defined as follows:

$-\log P(\bm{\varepsilon},\mathbf{z})=-\log
P(\bm{\varepsilon}|\mathbf{z})-\log P(\mathbf{z}).$ (7)

The first term, $-\log P(\bm{\varepsilon}|\mathbf{z})$, is the conditional
(negative log-) likelihood of the points residual conditioned on its latent
encoding. The second term, $-\log P(\mathbf{z})$, is the likelihood of
observing the latent encoding from the normal set.

We now detail our approach to interpret these terms into tractable, efficient
scores which we name the local reconstruction score and local density score
respectively. These scores are combined to give the overall ARES anomaly
score.

###  4.5 Local Reconstruction Score

The local reconstruction score is based on an estimate of how likely a given
residual $\bm{\varepsilon}$ (and consequent reconstruction error) is to come
from the corresponding sample $\mathbf{x}$, based on its latent encoding
$\mathbf{z}$:

$r(\mathbf{x})=-\log P(\bm{\varepsilon}|\mathbf{z})$ (8)

This likelihood cannot be calculated directly for any individual $\mathbf{z}$.
Instead, we consider the $k$ nearest neighbours of $\mathbf{z}$ in the latent
space, denoted $\mathsf{N}_{k}(\mathbf{z})$, to be a sample population which
$z$ belongs to. This is intuitive as data points with similar characteristics
to x in the input space are more likely to be encoded nearer to $\mathbf{z}$
in the latent space. The full local reconstruction score algorithm is shown in
Algorithm 1.

After training the autoencoder, we fix the model weights and store in memory
all training samples’ reconstruction errors and latent encodings taken from
the bottleneck layer. For a test sample $\mathbf{x}$, we find its own latent
encoding $\mathbf{z}$ and reconstruction $\hat{\mathbf{x}}$. We then find
$\mathsf{N}_{k}(\mathbf{z})$, the set of $k$ nearest neighbours to
$\mathbf{z}$ amongst the training set encodings. We only find nearest
neighbours among the training samples as they are all assumed to be normal and
therefore should have reconstruction errors within the normal range.

We want to obtain an estimate of the reconstruction error that could be
expected of $\mathbf{x}$, assuming its a normal point, based on the
reconstruction errors of its (normal) neighbours. We can then compare this
expected value, conditioned on its unique location in the latent space, with
the points true reconstruction error to determine its anomalousness.

In a fully probabilistic approach, we could do this by measuring the
likelihood of the test points reconstruction error under a probability
distribution, e.g. a Gaussian, fit to the neighbours’ reconstruction errors.
However, it is unnecessarily restrictive to assume any closed-form probability
distribution to adequately model this population for the unique neighbourhood
of each and every test sample. Instead, we opt for a non-parametric approach;
we measure the difference between the test points reconstruction error and the
median reconstruction error of its neighbouring samples. The larger the
difference between them, the more outlying the test point is in comparison to
its local neighbours, therefore the more likely it is to be anomalous. In
practice, using the median was found to perform better than the mean as it is
more robust to extrema. With this, we obtain the local reconstruction score
as:

$r(\mathbf{x})=\|\mathbf{x}-\mathbf{\hat{x}}\|_{2}^{2}-\underset{\mathbf{n}\in\mathsf{N}_{k}(\mathbf{z})}{\mathsf{median}}(\|\mathbf{n}-\hat{\mathbf{n}}\|_{2}^{2}).$
(9)

Nearest neighbour search in the latent space is preferable over the input
space as dimension reduction helps to alleviates the curse of dimensionality
(see Section 3), resulting in more semantically-meaningful neighbors.
Secondly, from a practical perspective, the neighbour search is less time-
consuming and computationally intensive in lower dimensional spaces.

Algorithm 1 Local Reconstruction Score

Input: Autoencoder $A(\cdot)=(f_{\theta}\circ g_{\phi})(\cdot)$, training set
$\mathbf{X}_{\text{train}}$, test sample $\mathbf{x}_{\text{test}}$
Parameters: neighbour count $k$
Output: Local reconstruction score $r(\mathbf{x}_{\text{test}})$

1: Train autoencoder $A$ on training set according to:
$\min_{\theta,\phi}\|\mathbf{X}_{\text{train}}-\mathbf{\hat{X}}_{\text{train}}\|_{2}^{2}$
where $\mathbf{\hat{x}}^{i}_{\text{train}}=(f_{\theta}\circ
g_{\phi})(\mathbf{x}_{\text{train}}^{i})$ for
$\mathbf{x}_{\text{train}}^{i}\in\mathbf{X}_{\text{train}}$.

2: Find latent encoding of $\mathbf{x}_{\text{test}}$:
$\mathbf{z}_{\text{test}}:=g_{\phi}(\mathbf{x}_{\text{test}})$

3: Find the set of $k$ nearest neighbours to $\mathbf{z}_{\text{test}}$ among
latent encodings of the training data:
$\mathsf{N}_{k}(\mathbf{z}_{\text{test}}):=\\{\mathbf{z}_{\text{train}}^{1},...,\mathbf{z}_{\text{train}}^{k}\\}$
where
$\mathbf{z}_{\text{train}}^{i}=g_{\phi}(\mathbf{x}^{i}_{\text{train}})\in$ for
$i=\\{1,...,k\\}$.

4: return
$r(\mathbf{x}_{\text{test}})=\|\mathbf{x}_{\text{test}}-\mathbf{\hat{x}}_{\text{test}}\|_{2}^{2}-\underset{\mathbf{n}\in\mathsf{N}_{k}(\mathbf{z})}{\mathsf{median}}(\|\mathbf{x}^{i}_{\text{train}}-\mathbf{\hat{x}}^{i}_{\text{train}}\|_{2}^{2})$

###  4.6 Local Density Score

The local reconstruction score corresponds to the conditional term of the
joint likelihood. We now introduce the local density score, which corresponds
to the likelihood of observing the given encoding in the latent space. This is
a density estimation task, which concerns the relative distance of
$\mathbf{z}$ to its nearest neighbours, unlike the local reconstruction score
which focuses on the reconstruction error of neighbours. Anomalies are assumed
to exist in sparse regions, where normal samples are unlikely to be found in
significant numbers. Any multivariate distribution $P$, with trainable
parameters $\Theta$, could be used to estimate this density:

$\displaystyle d(\mathbf{x}):=-\log P(\mathbf{z};\Theta)$ (10)

Note that it is common to ignore constant factor shifts in the anomaly score.
Thus, the distribution $P$ need not be normalized; even unnormalized density
estimation techniques can be used as scoring functions. We note that LOF is an
example of an unnormalized score which is similarly locally adaptive like the
local reconstruction score, so LOF is used for the local density score in our
main experiments. Other methods are also tested and their performance is shown
in the ablation study.

The overall anomaly score for sample $\mathbf{x}$ is:

$s(\mathbf{x}):=r(\mathbf{x})+\alpha d(\mathbf{x}),$ (11)

where $r(\mathbf{x})$ is its local reconstruction score and $d(\mathbf{x})$
the local density score. These two scores are unnormalized, so we use a
scaling factor $\alpha$ to balance the relative magnitudes of the two scores.
We heuristically set it equal to $0.5$ for all datasets and settings for
simplicity, as this was found to balance the two scores sufficiently fairly in
most cases. We choose not to treat $\alpha$ as a hyper-parameter to be tuned
to optimise performance, although different values could give better
performance for different datasets. The effect of changing $\alpha$ is shown
in the supplementary material.

#### Computational Runtime:

The average runtimes of experiments with the MNIST dataset can be found in the
supplementary material. We see that, despite taking longer than the standard
reconstruction error approach, the additional computational runtime of ARES
anomaly scoring is insignificant in relation to the model training time.
Anomaly scoring with ARES is just $1.2\%$ of the overall time taken (including
training) in the one-class case (0.017 minutes for scoring versus 1.474
minutes for training), and $3.5\%$ in the multi-class case. The additional
run-time is a result of the $k$ nearest neighbour search. An exact search
algorithms would be $\mathcal{O}(nm)$ with $n$ train and $m$ test samples,
however approximate methods can achieve near-exact accuracy much more
efficiently.

##  5 Experiments

In our experiments, we aim to answer the following research questions:

RQ1 (Accuracy): Does ARES perform better than existing anomaly detection
methods?
RQ2 (Ablation Study): How do different components and design choices of ARES
contribute to its performance?

###  5.1 Datasets and Experimental Setup

Dataset | #Dim | #Classes | #Samples | Description
---|---|---|---|---
SNSR [31] | $48$ | Multi-class | $58,509$ | Electric current signals
MNIST [19] | $784$ | Multi-class | $70,000$ |  $0$-$9$ digit images
FMNIST [32] | $784$ | Multi-class | $70,000$ | Fashion article images
OTTO [22] | $93$ | Multi-class | $61,878$ | E-commerce types
MI-F [29] | $58$ | Single-class | $25,286$ | CNC milling defects
MI-V [29] | $58$ | Single-class | $23,125$ | CNC milling defects
EOPT [16] | $20$ | Single-class | $90,515$ | Storage system failures
Table 1: Name and descriptions of the datasets used in experiments, including
the number of samples and dimensions.

Table 1 shows the datasets we use in experiments. The single-class datasets
consist of ground truth normal-vs-anomaly labels, as opposed to the multiple
class labels in multi-class datasets. In the latter, we distinguish between
the ’one-class normality’ setup, in which one class label is used as the
normal class and all other classes are anomalous. Alternatively, in the
’multi-class normality’ setup, one class is anomalous and all other classes
are normal. For a dataset with $N$ classes, there are $N$ possible
arrangements of normal and anomaly classes. We train separate models for each
arrangement and find their average score for the final result. We use the
Area-Under-Curve (AUC) metric to measure performance as it does not require an
anomaly score threshold to be set.

As anomaly scores are calculated for each test sample independently of each
other, the proportion of anomalies in the test set has no impact on the
anomaly score assigned to any given sample. Therefore, we are able to use a
normal:anomaly ratio of 50:50 in our experiments for the sake of simplicity
and an unbiased AUC metric. Besides the normal sampels in the test set, the
remaining normal samples are split $80$:$20$ into training and validation sets
for all models. Full implementation details can be found in the supplementary
material.

###  5.2 Baselines

We test the performance of ARES against a range of baselines. We use the sci-
kit learn implementations of LOF (in the input space) [8], IForest [20], PCA
[27] and OC-SVM [10]. Publicly available codes are used for DAGMM [7], RAPP-
SAP and RAPP-NAP [17], and we use Pytorch to build the autoencoder (AE and
ARES) and variational autoencoder (VAE). All experiments were conducted in
Windows OS using an Nvidia GeForce RTX 2080 Ti GPU.

We do not tune hyper-parameters relating to the model architectures or
training procedures for any method. The effect of variation in hyper-
parameters is studied in the ablation study in Section 5.4 and the
supplementary material instead. We set the number of neighbours $k=10$ for
both the local density and local reconstruction score in our main experiments.

###  5.3 RQ1 (Accuracy):

Table 2 shows the average AUC scores as a percentage (i.e. multiplied by 100).
In the one-class normality setting, ARES significantly improves performance
over the baselines in all multi-class datasets besides FMNIST. In the single-
class datasets, this improvement is even greater, e.g. $+8\%$ lift for MI-F
and EOPT. Compared with AE, we see that local adaptivity helps to detect true
anomalies by correcting for the natural variation in reconstruction error in
the latent space. This effect is even more pronounced in the multi-class
normality setting, where ARES gives the best performance on all datasets.
Here, the normal set is much more diverse and therefore local adaptivity is
even more important. We show the standard deviations and additional
significance test scores in the supplementary material, which also show
statistically significant ($p<0.01$) improvement.

Dataset | LOF | IForest | PCA | OC-SVM | SAP | NAP | DAGMM | VAE | AE | ARES
---|---|---|---|---|---|---|---|---|---|---
One-class Normality
SNSR | $97.98$ | $89.16$ | $92.01$ | $95.85$ | $98.79$ | $98.74$ | $88.08$ | $89.49$ | $98.30$ |  $\mathbf{98.83}$**
MNIST | $96.85$ | $85.44$ | $95.68$ | $90.35$ | $95.35$ | $97.25$ | $89.60$ | $91.73$ | $96.96$ |  $\mathbf{97.89}$**
FMNIST | $91.35$ | $91.39$ | $90.13$ | $90.74$ | $89.66$ | $\mathbf{93.08}$ | $87.97$ | $77.59$ | $92.33$ | $91.63$
OTTO | $84.76$ | $70.34$ | $80.09$ | $81.43$ | $81.61$ | $82.77$ | $68.02$ | $82.39$ | $85.26$ |  $\mathbf{87.86}$**
MI-F | $59.79$ | $81.53$ | $55.07$ | $76.69$ | $81.78$ | $80.61$ | $82.23$ | $76.93$ | $71.19$ |  $\mathbf{89.52}$
MI-V | $83.97$ | $84.35$ | $87.32$ | $83.58$ | $88.24$ | $89.35$ | $75.45$ | $89.03$ | $90.75$ |  $\mathbf{93.94}$
EOPT | $55.01$ | $61.61$ | $54.72$ | $59.66$ | $59.87$ | $61.69$ | $60.63$ | $68.08$ | $59.85$ |  $\mathbf{68.43}$
Multi-class Normality
SNSR | $60.74$ | $52.70$ | $52.94$ | $52.79$ | $57.52$ | $58.32$ | $54.77$ | $61.36$ | $57.28$ |  $\mathbf{69.78}$**
MNIST | $77.40$ | $56.49$ | $70.41$ | $58.56$ | $84.76$ | $86.38$ | $54.24$ | $84.12$ | $80.04$ |  $\mathbf{93.25}$**
FMNIST | $71.50$ | $64.75$ | $66.92$ | $60.27$ | $68.50$ | $72.09$ | $57.56$ | $71.18$ | $71.03$ |  $\mathbf{72.49}$
OTTO | $63.01$ | $54.14$ | $58.27$ | $62.96$ | $57.47$ | $63.44$ | $58.96$ | $61.88$ | $59.59$ |  $\mathbf{63.54}$
Table 2: Mean AUC scores for each datasets and normality setting. The best
scores are highlighted in bold and we mark the most significant improvements
over AE ($p<0.01$) with **. Further tests are in the supplementary material.

The neighbours of a normal sample with high reconstruction error tend to be
have high reconstruction errors themselves. By basing the anomaly scoring on
their relative difference, ARES uses this to better detect truly anomalous
samples. Furthermore, ARES also uses the local density of the point, which
depends purely on its distance to the training samples in the latent space.
This is important as there may be some anomalies with such low reconstruction
error that comparison with neighbours does not alone indicate anomalousness
(for example the weekend samples mentioned earlier). These samples could be
expected to be occupy very sparse regions of the latent space due to their
significant deviation from the normal set, which means they can be better
detected through the local density score. By combining these two scores, we
are able to detect a wider range of anomalous data than either could
individually.

###  5.4 RQ2 (Ablation Study):

#### Density Estimation Method

Table 3 shows the performance of ARES with other density estimation methods.
KNN is the distance to the $k^{th}$ nearest neighbour in the latent space with
$k=20$. GD is the distance of a point to the closest of $N$ Gaussian
distributions fit to the latent encodings of samples for each of the classes
in the training set ($N$ = 1 in the one-class normality case). NF is the
likelihood under a RealNVP normalizing flow [23].

Dataset |  LOF |  KNN |  GD |  NF
---|---|---|---|---
One-class Normality
SNSR |  $\mathbf{98.83}$ |  $98.66$ |  $95.68$ |  $98.29$
MNIST |  $\mathbf{97.89}$ |  $95.24$ |  $87.29$ |  $97.30$
FMNIST |  $91.63$ |  $\mathbf{92.94}$ |  $90.26$ |  $92.14$
OTTO |  $\mathbf{87.76}$ |  $83.14$ |  $81.27$ |  $86.74$
MI-F |  $\mathbf{89.52}$ |  $76.12$ |  $80.41$ |  $84.47$
MI-V |  $\mathbf{93.94}$ |  $91.60$ |  $79.61$ |  $92.89$
EOPT |  $\mathbf{68.43}$ |  $67.63$ |  $63.35$ |  $61.07$
Multi-class Normality
SNSR |  $\mathbf{69.78}$ |  $67.42$ |  $61.86$ |  $60.63$
MNIST |  $\mathbf{93.25}$ |  $92.92$ |  $91.14$ |  $86.25$
FMNIST |  $72.49$ |  $72.77$ |  $\mathbf{73.02}$ |  $70.48$
OTTO |  $63.54$ |  $61.71$ |  $\mathbf{66.30}$ |  $63.67$
Table 3: Mean AUC scores for each choice of local density score. The best
scores are highlighted in bold.

LOF performs best overall, closely followed by KNN. GD performs poorly in one-
class experiments, however it is better in multi-class normality experiments
and even the best for FMNIST and OTTO. This could be as the use of multiple
distributions provides more flexbility. NF is noticeably worse; previous
studies have found that normalizing flows are not well-suited to detect out-
of-distribution data [18].

In the supplementary material, we further test these density estimation
methods by varying their hyper-parameters and find LOF to still come out best.
Further supplementary experiments show that the local density score generally
performs better than the local reconstruction score in the multi-class
normality case and vice versa in the one-class case. Combining them, as in
ARES, generally gives the best performance overall across different latent
embedding sizes.

#### Robustness to Training Contamination

In practice, it is likely that a small proportion of anomalies ‘contaminate’
the training set. In this section, we study the effect of different levels of
contamination on ARES, defined as $n\%$ of the total number of samples in the
training set. Table 4 shows that increasing $n$ worsens performance overall.
With more anomalies in the training set, it is more likely that anomalies are
found in the nearest neighbour set of more test points, which skews both the
average neighbourhood reconstruction error as well as their density in the
latent space and degrades performance.

Neighbourhood Size
---
$n$ | $10$ | $50$ | $100$ | $200$ | $500$
One-class Normality
$0$ | $\mathbf{97.89}$ | $97.85$ | $97.80$ | $97.73$ | $97.56$
$0.5$ | $95.77$ | $\mathbf{96.08}$ | $96.06$ | $95.95$ | $95.72$
$1$ | $82.50$ | $94.45$ | $\mathbf{95.97}$ | $95.81$ | $95.42$
$2$ | $90.65$ | $92.88$ | $93.29$ | $\mathbf{93.33}$ | $92.99$
$3$ | $88.41$ | $91.35$ | $91.92$ | $\mathbf{92.07}$ | $91.86$
$5$ | $84.91$ | $87.85$ | $89.20$ | $90.26$ | $\mathbf{90.65}$
$10$ | $79.69$ | $81.51$ | $83.16$ | $84.88$ | $\mathbf{86.82}$
Multi-class Normality
$0$ | $\mathbf{93.25}$ | $92.71$ | $92.08$ | $91.31$ | $89.68$
$0.5$ | $68.38$ | $76.07$ | $78.36$ | $80.75$ | $\mathbf{81.42}$
$1$ | $52.38$ | $58.42$ | $67.62$ | $65.75$ | $\mathbf{68.91}$
$2$ | $59.50$ | $63.13$ | $65.71$ | $67.98$ | $\mathbf{69.68}$
$3$ | $57.09$ | $59.00$ | $60.78$ | $62.89$ | $\mathbf{65.35}$
$5$ | $50.86$ | $51.48$ | $52.40$ | $53.47$ | $\mathbf{55.78}$
$10$ | $\mathbf{52.50}$ | $50.80$ | $50.87$ | $51.14$ | $52.11$
Table 4: Mean AUC scores in MNIST with training set anomaly contamination
($n$%) and different neighbourhood sizes. The best scores are highlighted in
bold.

We find that ARES is more robust to training set contamination with a higher
setting of the number of neighbours ($k$). In the one-class setup, we see that
ARES performs better with higher values of $k$ as the proportion of anomalies
increases. By using more neighbours, the effect of any individual anomalies on
the overall neighbourhood error is reduced, which helps to maintain better
performance. This effect is even more stronger in the multi-class normality
setup. The highest neighbour count of $k=500$ gives the best performance in
all cases except for $n=0\%$ and $10\%$. As the multi-class training sets are
much larger than their one-class normality counterparts, $n\%$ corresponds to
a much larger number of anomalous contaminants, which explains their greater
effect for a given $k$.

##  6 Conclusion

Autoencoders are extremely popular deep learning models used for anomaly
detection through their reconstruction error. We have shown that the
assumption made by the standard reconstruction error score, that
reconstruction errors are identically distributed for all normal samples, is
unsuitable for real datasets. We empirically show that there is a
heteroscedastic relationship between latent space characteristics and
reconstruction error, which demonstrates why adaptivity to local latent
information is important for anomaly scoring. As such, we have developed a
novel approach to anomaly scoring which adaptively evaluates the anomalousness
of a samples reconstruction error, as well as its density in the latent space,
relative to those of its nearest neighbours. We show that our approach results
in significant performance improvements over the standard approach, as well as
other prominent baselines, across a range of real datasets.

## Acknowledgements

This work was supported in part by NUS ODPRT Grant R252-000-A81-133.

## References

  * [1] Abati, D., Porrello, A., Calderara, S., Cucchiara, R.: Latent space autoregression for novelty detection. In: ICCV. pp. 481–490 (2019)
  * [2] Akcay, S., Atapour-Abarghouei, A., Breckon, T.P.: Ganomaly: Semi-supervised anomaly detection via adversarial training. In: ACCV. pp. 622–637. Springer (2018)
  * [3] Amarbayasgalan, T., Jargalsaikhan, B., Ryu, K.H.: Unsupervised novelty detection using deep autoencoders with density based clustering. Applied Sciences 8(9), 1468 (2018)
  * [4] An, J.: Variational Autoencoder based Anomaly Detection using Reconstruction Probability. In: SNU Data Mining Center 2015-2 Special Lecture on IE (2015)
  * [5] Bergman, L., Cohen, N., Hoshen, Y.: Deep nearest neighbor anomaly detection. arXiv preprint arXiv:2002.10445 (2020)
  * [6] Beyer, K., Goldstein, J., Ramakrishnan, R., Shaft, U.: When is “nearest neighbor” meaningful? In: International conference on database theory. pp. 217–235. Springer (1999)
  * [7] Bo, Z., Song, Q., Chen, H.: Deep autoencoding gaussian mixture model for unsupervised anomaly detection (2018)
  * [8] Breunig, M.M., Kriegel, H.P., Ng, R.T., Sander, J.: LOF: identifying density-based local outliers. In: SIGMOD. vol. 29, pp. 93–104. ACM (2000)
  * [9] Chen, J., Sathe, S., Aggarwal, C., Turaga, D.: Outlier detection with autoencoder ensembles. In: SDM. pp. 90–98. SIAM (2017)
  * [10] Chen, Y., Zhou, X.S., Huang, T.S.: One-class SVM for learning in image retrieval. In: ICIP. pp. 34–37. Citeseer (2001)
  * [11] Deng, A., Goodge, A., Lang, Y.A., Hooi, B.: CADET: Calibrated Anomaly Detection for Mitigating Hardness Bias. IJCAI (2022)
  * [12] Dinh, L., Sohl-Dickstein, J., Bengio, S.: Density estimation using real nvp. arXiv preprint arXiv:1605.08803 (2016)
  * [13] Feng, W., Han, C.: A Novel Approach for Trajectory Feature Representation and Anomalous Trajectory Detection. ISIF pp. 1093–1099 (2015)
  * [14] Goodge, A., Hooi, B., Ng, S.K., Ng, W.S.: Robustness of Autoencoders for Anomaly Detection Under Adversarial Impact. IJCAI (2020)
  * [15] Goodge, A., Hooi, B., Ng, S.K., Ng, W.S.: Lunar: Unifying local outlier detection methods via graph neural networks. In: Proceedings of the AAAI Conference on Artificial Intelligence (2022)
  * [16] inIT: Tool wear detection in cnc mill (2018), https://www.kaggle.com/init-owl/high-storage-system-data-for-energy-optimization
  * [17] Kim, K.H., Shim, S., Lim, Y., Jeon, J., Choi, J., Kim, B., Yoon, A.S.: RaPP: Novelty Detection with Reconstruction along Projection Pathway. In: ICLR (2019)
  * [18] Kirichenko, P., Izmailov, P., Wilson, A.G.: Why normalizing flows fail to detect out-of-distribution data. arXiv preprint arXiv:2006.08545 (2020)
  * [19] Lecun, Y.: Mnist (2012), http://yann.lecun.com/exdb/mnist/
  * [20] Liu, F.T., Ting, K.M., Zhou, Z.H.: Isolation-based anomaly detection. ACM Transactions on Knowledge Discovery from Data (TKDD) 6(1), 1–39 (2012)
  * [21] Maaten, L.v.d., Hinton, G.: Visualizing data using t-sne. Journal of machine learning research 9(Nov), 2579–2605 (2008)
  * [22] Otto, G.: Otto group product classification challenge (2015), https://www.kaggle.com/c/otto-group-product-classification-challenge
  * [23] Papamakarios, G., Pavlakou, T., Murray, I.: Masked autoregressive flow for density estimation. In: NeurIPS. pp. 2338–2347 (2017)
  * [24] Rezende, D.J., Mohamed, S.: Variational inference with normalizing flows. arXiv preprint arXiv:1505.05770 (2015)
  * [25] Ruff, L., Vandermeulen, R., Goernitz, N., Deecke, L., Siddiqui, S.A., Binder, A., Müller, E., Kloft, M.: Deep one-class classification. In: ICML. pp. 4393–4402 (2018)
  * [26] Sakurada, M., Yairi, T.: Anomaly detection using autoencoders with nonlinear dimensionality reduction. In: MLSDA. p. 4. ACM (2014)
  * [27] Shyu, M.L., Chen, S.C., Sarinnapakorn, K., Chang, L.: A novel anomaly detection scheme based on principal component classifier. Tech. rep., Miami Univ Coral Gables FL Dept of Electric and Computer Engineering (2003)
  * [28] Siffer, A., Fouque, P.A., Termier, A., Largouet, C.: Anomaly detection in streams with extreme value theory. In: SIGKDD. pp. 1067–1075 (2017)
  * [29] SMART: Tool wear detection in cnc mill (2018), https://www.kaggle.com/shasun/tool-wear-detection-in-cnc-mill
  * [30] Tax, D.M., Duin, R.P.: Support vector data description. Machine learning 54(1), 45–66 (2004)
  * [31] UCI: Sensorless drive diagnosis (2015)
  * [32] Xiao, H., Rasul, K., Vollgraf, R.: Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms (2017)
  * [33] Zenati, H., Foo, C.S., Lecouat, B., Manek, G., Chandrasekha, V.R.: Efficient gan-based anomaly detection (2019)
  * [34] Zhai, S., Cheng, Y., Lu, W., Zhang, Z.: Deep structured energy based models for anomaly detection. ICML 48, 1100–1109 (2016)
  * [35] Zimek, A., Gaudet, M., Campello, R.J.G.B., Sander, J.: Subsampling for efficient and effective unsupervised outlier detection ensembles. In: SIGKDD. pp. 428–436 (2013)


$\displaystyle\times\Big{[}e^{iy(p-l+k)}e^{-i\bar{y}(p+m+k)}e^{iz(l+m)}+e^{iy[p-l+(1-\xi)(k-q)]}e^{-i\bar{y}[p+m+(1-\xi)(k-q)]}e^{iz(l+m)}$
$\displaystyle-e^{iy(p-l+k)}e^{-i\bar{y}[p+m+(1\xi)(k-q)]}e^{iz[l+m-\xi
k-(1-\xi)q]}-e^{iy[p-l-(1-\xi)(k-q)]}e^{-i\bar{y}(p+m+q)}e^{iz[l+m+(1-\xi)k+\xi
q]}\Big{]}.$

After performing the integrals, this contribution reads

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q}_{\rm NLO,r}$
$\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{\xi_{0}}d\xi\,\frac{x_{p}}{(1-\xi)}\,f^{g}_{\mu_{0}^{2}}\left(\frac{x_{p}}{1-\xi}\right)\frac{1}{2}[\xi^{2}+(1-\xi)^{2}]\int_{k,q>\mu_{0}}s(k)s(q)$
(95)
$\displaystyle\times\,\bigg{\\{}\frac{1}{(p-q)^{2}}+\frac{1}{[p-(1-\xi)(q-k)]^{2}}-2\frac{(p-q)^{i}}{(p-q)^{2}}\frac{[p-(1-\xi)(q-k)]^{i}}{[p-(1-\xi)(q-k)]^{2}}\bigg{\\}}.$

At this point we introduce the fragmentation function and reorganize the
terms:

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,r}=\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{q}_{\mu_{0}^{2}}(\zeta)\int_{k,q>\mu_{0}}\int_{\xi_{0}}d\xi\,\frac{x_{F}}{\zeta(1-\xi)}\,f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta(1-\xi)}\right)\frac{1}{2}\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}s(k)s(q)$
$\displaystyle\times\bigg{\\{}\frac{1}{2}\frac{(k-q)^{2}}{(p/\zeta-k)^{2}(p/\zeta-q)^{2}}+\frac{1}{2}\int_{t}s(t)\frac{(1-\xi)^{2}[(q-k)-(q-t)]^{2}}{[p/\zeta-(1-\xi)(q-k)]^{2}[p/\zeta-(1-\xi)(q-t)]^{2}}\bigg{\\}}+({\rm
Gen.\,NLO})_{2}\ .$ (96)

The first term in eq. (D.1.2) is a contribution to the quark TMD PDF, the
second term is a contribution to gluon TMD FF and the last term is a genuine
NLO correction that do not include any large logarithms. The explicit
expression for the genuine NLO contribution is

$\displaystyle({\rm Gen.\,NLO})_{2}$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{q}_{\mu_{0}^{2}}(\zeta)\int_{k,q>\mu_{0}}\int_{\xi_{0}}d\xi\,\frac{x_{F}}{\zeta(1-\xi)}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta(1-\xi)}\right)\frac{1}{2}\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}s(k)s(q)$
(97)
$\displaystyle\times\bigg{\\{}\frac{(p/\zeta-q)^{i}}{(p/\zeta-q)^{2}}\bigg{[}\frac{(p/\zeta-k)^{i}}{(p/\zeta-k)^{2}}-\frac{[p/\zeta-(1-\xi)(q-k)]^{i}}{[p/\zeta-(1-\xi)(q-k)]^{2}}\bigg{]}$
$\displaystyle\hskip
17.07182pt+\int_{t}s(t)\frac{[p/\zeta-(1-\xi)(q-k)]^{i}}{[p/\zeta-(1-\xi)(q-k)]^{2}}\bigg{[}\frac{[p/\zeta-(1-\xi)(q-t)]^{i}}{[p/\zeta-(1-\xi)(q-t)]^{2}}-\frac{(p/\zeta-q)^{i}}{(p/\zeta-q)^{2}}\bigg{]}\bigg{\\}}.$

Using the same arguments introduced in eq. (91), the final expression for this
contribution can be written as

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,r}$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{q}_{\mu_{0}^{2}}(\zeta)\int_{k,q>\mu_{0}}\int_{\xi_{0}}d\xi\,\frac{x_{F}}{\zeta(1-\xi)}\,f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta(1-\xi)}\right)\frac{1}{2}\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\,\frac{1}{k^{2}}$
(98) $\displaystyle\times\,s(-k+p/\zeta)\,\bigg{[}1-\frac{k\cdot
q}{q^{2}}\bigg{]}\,s(-q+p/\zeta)$ $\displaystyle+$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,\frac{1}{(1-\xi)}D^{q}_{\mu_{0}^{2}}\left(\frac{\zeta}{1-\xi}\right)\int_{k,q,t>\mu_{0}}\int_{\xi_{0}}d\xi\,\frac{x_{F}}{\zeta}\,f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\frac{1}{2}\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\,\frac{1}{k^{2}}$
$\displaystyle\times\,s(t)s(p/\zeta+k-t)\,\bigg{[}1-\frac{k\cdot
q}{q^{2}}\bigg{]}\,s(p/\zeta+q-t)$ $\displaystyle+$ $\displaystyle({\rm
Gen.\,NLO})_{2}\ .$

####  D.1.3 Gluon initiated gluon production

The partonic level expression for this contribution is given by eq. (4.9) of
Ref. [37] which, at large $N_{c}$, reads

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm
NLO,r}=\frac{g^{2}}{(2\pi)^{3}}N_{c}\int_{\xi_{0}}d\xi\frac{x_{p}}{(1-\xi)}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{p}}{1-\xi}\right)\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}\int_{y,\bar{y},z}A^{i}_{\xi}(y-z)A^{i}_{\xi}(\bar{y}-z)$
$\displaystyle\times\Big{\\{}2s(y,\bar{y})s(\bar{y},y)+2s[(1-\xi)y,(1-\xi)\bar{y}]s[(1-\xi)\bar{y},(1-\xi)y]-s(y,z)s[(1-\xi)\bar{y}+\xi
z,y]s[z,(1-\xi)\bar{y}+\xi z]$ $\displaystyle\hskip
11.38092pt-s(z,y)s[y,(1-\xi)\bar{y}+\xi z]s[(1-\xi)\bar{y}+\xi
z,z]-s[\bar{y},(1-\xi)y+\xi z]s[(1-\xi)y+\xi z,z]s(z,\bar{y})$
$\displaystyle\hskip 11.38092pt-s[(1-\xi)y+\xi z,\bar{y}]s[z,(1-\xi)y+\xi
z]s(\bar{y},z)\Big{\\}}.$ (99)

Again using the translationally invariant dipoles and the definition of the
modified WW fields, we get

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm
NLO,r}=\frac{g^{2}}{(2\pi)^{3}}N_{c}\int_{\xi_{0}}d\xi\frac{x_{p}}{(1-\xi)}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{p}}{1-\xi}\right)\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}\frac{(-1)}{(2\pi)^{4}}\int_{y,\bar{y},z}\int_{l,m>\mu_{0}}\int_{k,q,t}\frac{l^{i}}{l^{2}}\frac{m^{i}}{m^{2}}$
$\displaystyle\times
s(k)s(q)s(t)\Big{\\{}2e^{iy(p-l+k-q)}e^{-i\bar{y}(p+m+k-q)}e^{iz(l+m)}+2e^{iy[p-l+(1-\xi)(k-q)]}e^{-i\bar{y}[p+m+(1-\xi)(k-q)]}e^{iz(l+m)}$
(100)
$\displaystyle-e^{iy(p-l+t-k)}e^{-i\bar{y}[p+m+(1-\xi)(q-k)]}e^{iz[l+m-t+\xi
k+(1-\xi)q]}-e^{iy(p-l-t+k)}e^{-i\bar{y}[p+m-(1-\xi)(q-k)]}e^{iz[l+m+t-\xi
k-(1-\xi)q]}$
$\displaystyle-e^{iy[p-l+(1-\xi)(q-k)]}e^{-i\bar{y}(p+m-k+t)}e^{iz[l+m+t-\xi
k-(1-\xi)q]}-e^{iy[p-l+(1-\xi)(k-q)]}e^{-i\bar{y}(p+m+k-t)}e^{iz[l+m-t+\xi
k+(1-\xi)q]}\Big{\\}}.$

After performing the integrals, this contribution can be written as

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm NLO,r}$
$\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{\xi_{0}}d\xi\frac{x_{p}}{(1-\xi)}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{p}}{1-\xi}\right)\,2N_{c}\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}\int_{k,q,t}s(k)\,s(q)\,s(t)$
(101)
$\displaystyle\times\bigg{\\{}\frac{1}{[p-(q-k)]^{2}}+\frac{1}{[p-(1-\xi)(q-k)]^{2}}-2\frac{(p-q+t)^{i}}{(p-q+t)^{2}}\frac{[p-(1-\xi)(q-k)]^{i}}{[p-(1-\xi)(q-k)]^{2}}\bigg{\\}}.$

Introducing the fragmentation function and using the symmetry between $k,q$
and $t$, it can be organized as

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g\to H}_{\rm
NLO,r}=\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{g}_{\mu_{0}^{2}}(\zeta)\int_{k>\mu_{0}}\int_{\xi_{0}}d\xi\frac{x_{F}}{\zeta(1-\xi)}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta(1-\xi)}\right)\,2N_{c}\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}$
$\displaystyle\times\int_{q,t}s(k)\,s(q)\,s(t)\bigg{\\{}\frac{1}{2}\frac{[(q-k)-(q-t)]^{2}}{[p/\zeta-(q-t)]^{2}[p/\zeta-(q-k)]^{2}}+\frac{1}{2}\frac{(1-\xi)^{2}[(q-k)-(q-t)]^{2}}{[p/\zeta-(1-\xi)(q-t)]^{2}[p/\zeta-(1-\xi)(q-k)]^{2}}\bigg{\\}}$
$\displaystyle+({\rm Gen.\,NLO})_{3}\ ,$ (102)

where the genuine NLO contribution reads

$\displaystyle({\rm
Gen.\,NLO})_{3}=\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{g}_{\mu_{0}^{2}}(\zeta)\int_{k>\mu_{0}}\int_{\xi_{0}}d\xi\frac{x_{F}}{\zeta(1-\xi)}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta(1-\xi)}\right)\,2N_{c}\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}$
$\displaystyle\times\int_{q,t}s(k)\,s(q)\,s(t)\bigg{[}\frac{[p/\zeta-(q-t)]^{i}}{[p/\zeta-(q-t)]^{2}}-\frac{[p/\zeta-(1-\xi)(q-t)]^{i}}{[p/\zeta-(1-\xi)(q-t)]^{2}}\bigg{]}\bigg{[}\frac{[p/\zeta-(q-k)]^{i}}{[p/\zeta-(q-k)]^{2}}-\frac{[p/\zeta-(1-\xi)(q-k)]^{i}}{[p/\zeta-(1-\xi)(q-k)]^{2}}\bigg{]}.$
(103)

One can again massage the first two terms in eq. (D.1.3) using similar
arguments to those introduced in the previous cases and the final result reads

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g\to H}_{\rm
NLO,r}$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{g}_{\mu_{0}^{2}}(\zeta)\int_{k>\mu_{0}}\int_{q,t}\int_{\xi_{0}}d\xi\frac{x_{F}}{\zeta(1-\xi)}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta(1-\xi)}\right)\,2N_{c}\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}$
(104)
$\displaystyle\times\,\frac{1}{k^{2}}\,s(t)\,s(p/\zeta+k-t)\bigg{[}1-\frac{k\cdot
q}{q^{2}}\bigg{]}\,s(p/\zeta+q-t)$ $\displaystyle+$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,\frac{1}{1-\xi}\,D^{g}_{\mu_{0}^{2}}(\zeta)\int_{k>\mu_{0}}\int_{q,t}\int_{\xi_{0}}d\xi\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,2N_{c}\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}$
$\displaystyle\times\,\frac{1}{k^{2}}\,s(t)\,s(p/\zeta+k-t)\bigg{[}1-\frac{k\cdot
q}{q^{2}}\bigg{]}\,s(p/\zeta+q-t)$ $\displaystyle+$ $\displaystyle({\rm
Gen.\,NLO})_{3}\ .$

####  D.1.4 Quark initiated gluon production

At large $N_{c}$, this contribution reads (eq. (4.8) in Ref. [37]):

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm q\to g}_{\rm NLO,r}$
$\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\int_{\xi_{0}}d\xi\,\frac{x_{p}}{\xi}\,f^{q}_{\mu_{0}^{2}}\left(\frac{x_{p}}{\xi}\right)\frac{N_{c}}{2}\bigg{[}\frac{1+(1-\xi)^{2}}{\xi}\bigg{]}\int_{z\bar{z}y}e^{ip(z-\bar{z})}A^{i}_{\xi}(z-y)A^{i}_{\xi}(\bar{z}-y)$
(105) $\displaystyle\times\Big{\\{}s(z,\bar{z})s(\bar{z},z)+s(\xi
z,\xi\bar{z})-s[z,(1-\xi)y+\xi\bar{z}]s(y,z)-s(\bar{z},y)s[(1-\xi)y+\xi
z,\bar{z}]\Big{\\}}.$

. Using the same arguments, we get

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm q\to g}_{\rm NLO,r}$
$\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\int_{\xi_{0}}d\xi\,\frac{x_{p}}{\xi}\,f^{q}_{\mu_{0}^{2}}\left(\frac{x_{p}}{\xi}\right)\frac{N_{c}}{2}\bigg{[}\frac{1+(1-\xi)^{2}}{\xi}\bigg{]}\frac{(-1)}{(2\pi)^{4}}\int_{z\bar{z}y}\int_{l,m>\mu_{0}}\int_{k,q}\frac{l^{i}}{l^{2}}\frac{m^{i}}{m^{2}}\,s(k)\,s(q)$
$\displaystyle\times\Big{\\{}e^{iz(p-l+k-q)}e^{-i\bar{z}(p+m+k-q)}e^{iy(l+m)}+e^{iz(p-l+\xi
q)}e^{-i\bar{z}(p+m+\xi q)}e^{iy(l+m)}$
$\displaystyle-e^{iz(p-l+q-k)}e^{-i\bar{z}(p+m+\xi
q)}e^{iy[l+m+k-(1-\xi)q]}-e^{iz(p-l+\xi
k)}e^{-i\bar{z}(p+m+k-q)}e^{iy[l+m-q+(1-\xi)k]}\Big{\\}}.$

After performing the integrals, the contribution to the partonic level cross
section from this channel reads

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm q\to g}_{\rm NLO,r}$
$\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{\xi_{0}}d\xi\,\frac{x_{p}}{\xi}\,f^{q}_{\mu_{0}^{2}}\left(\frac{x_{p}}{\xi}\right)\frac{N_{c}}{2}\bigg{[}\frac{1+(1-\xi)^{2}}{\xi}\bigg{]}\int_{k,q>\mu_{0}}s(k)\,s(q)$
(107) $\displaystyle\times\bigg{\\{}\frac{1}{[p-(q-k)]^{2}}+\frac{1}{(p-\xi
q)^{2}}-2\frac{[p-(q-k)]^{i}}{[p-(q-k)]^{2}}\frac{(p-\xi q)^{i}}{(p-\xi
q)^{2}}\bigg{\\}}.$

After introducing the fragmentation function and symmetrizing our result, we
get

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm q\to g\to H}_{\rm
NLO,r}$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{g}_{\mu_{0}^{2}}(\zeta)\int_{k,q>\mu_{0}}\int_{\xi_{0}}d\xi\,\frac{x_{F}}{\zeta\xi}\,f^{q}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta\xi}\right)\frac{N_{c}}{2}\bigg{[}\frac{1+(1-\xi)^{2}}{\xi}\bigg{]}s(k)s(q)$
$\displaystyle\times\bigg{\\{}\frac{1}{2}\int_{t}s(t)\frac{[(q-k)-(q-t)]^{2}}{[p/\zeta-(q-k)]^{2}[p/\zeta-(q-t)]^{2}}+\frac{1}{2}\frac{\xi^{2}(q-k)^{2}}{(p/\zeta-\xi
q)^{2}(p/\zeta-\xi k)^{2}}\Bigg{\\}}+({\rm Gen.\,NLO})_{4}\ ,$

with

$\displaystyle({\rm
Gen.\,NLO})_{4}=\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{g}_{\mu_{0}^{2}}(\zeta)\int_{k,q>\mu_{0}}\int_{\xi_{0}}d\xi\,\frac{x_{F}}{\zeta\xi}\,f^{q}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta\xi}\right)\frac{N_{c}}{2}\bigg{[}\frac{1+(1-\xi)^{2}}{\xi}\bigg{]}s(k)s(q)$
(109)
$\displaystyle\times\bigg{\\{}\int_{t}s(t)\frac{[p/\zeta-(q-k)]^{i}}{[p/\zeta-(q-k)]^{2}}\bigg{[}\frac{[p/\zeta-(q-t)]^{i}}{[p/\zeta-(q-t)]^{2}}-\frac{(p/\zeta-\xi
q)^{i}}{(p/\zeta-\xi q)^{2}}\bigg{]}+\frac{(p/\zeta-\xi q)^{i}}{(p/\zeta-\xi
q)^{2}}\bigg{[}\frac{(p/\zeta-\xi k)^{i}}{(p/\zeta-\xi
k)^{2}}-\frac{[p/\zeta-(q-k)]^{i}}{[p/\zeta-(q-k)]^{2}}\bigg{]}\bigg{\\}}.$

Performing the same manipulations for the first two terms in eq. (D.1.4), we
get the final expression as

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm q\to g\to H}_{\rm
NLO,r}$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{g}_{\mu_{0}^{2}}(\zeta)\int_{k>\mu_{0}}\int_{q,t}\int_{\xi_{0}}d\xi\,\frac{x_{F}}{\zeta\xi}\,f^{q}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta\xi}\right)\frac{N_{c}}{2}\bigg{[}\frac{1+(1-\xi)^{2}}{\xi}\bigg{]}\frac{1}{k^{2}}$
(110) $\displaystyle\times\,s(t)\,s(p/\zeta+k-t)\bigg{[}1-\frac{k\cdot
q}{q^{2}}\bigg{]}\,s(p/\zeta+q-t)$ $\displaystyle+$
$\displaystyle\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,\frac{1}{\xi}D^{g}_{\mu_{0}^{2}}\left(\frac{\zeta}{\xi}\right)\int_{k>\mu_{0}}\int_{q}\int_{\xi_{0}}d\xi\,\frac{x_{F}}{\zeta}\,f^{q}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\frac{N_{c}}{2}\bigg{[}\frac{1+(1-\xi)^{2}}{\xi}\bigg{]}\frac{1}{k^{2}}$
$\displaystyle\times\,s(-k+p/\zeta)\bigg{[}1-\frac{k\cdot
q}{q^{2}}\bigg{]}\,s(-q+p/\zeta)$ $\displaystyle+$ $\displaystyle({\rm
Gen.\,NLO})_{4}\ .$

####  D.1.5 Summing all real contributions

The definitions of full quark and gluon TMD PDFs and FFs are given in eqs.
(14), (16), (19) and (21). Here we ignore the contribution from antiquarks
that, as indicated in the main text, can be added straightforwardly. The
contributions to the quark channel arise from the first two terms in eq. (92),
the first term in eq. (92) and the second term in eq. (110). Adding these
terms and using the definitions for the quark TMD PDF and TMD FF, we get

$\displaystyle\frac{d\sigma}{d^{2}pd\eta}\bigg{|}_{\rm NLO,r}^{q}$
$\displaystyle=$ $\displaystyle
S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\int_{k>\mu_{0}}\bigg{[}D^{q}_{\mu_{0}^{2}}(\zeta)\frac{x_{F}}{\zeta}{\cal
T}_{q}\left(\frac{x_{f}}{\zeta},k^{2};k^{2},\xi_{0}\right)+\frac{x_{F}}{\zeta}f^{q}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right){\cal
F}_{q}(\zeta,k^{2};k^{2};\xi_{0})\bigg{]}$ (111)
$\displaystyle\times\,\int_{q}s(-k+p/\zeta)\,\bigg{[}1-\frac{k\cdot
q}{q^{2}}\bigg{]}\,s(-q+p/\zeta).$

The contributions to the gluon channel arise from the second term in eq. (98),
the first two terms in eq. (104) and the first term in eq. (110). Adding those
terms and using the definitions for the gluon TMD PDF and TMD FF, we get

$\displaystyle\frac{d\sigma}{d^{2}pd\eta}\bigg{|}_{\rm NLO,r}^{g}$
$\displaystyle=$ $\displaystyle
S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\int_{k>\mu_{0}}\bigg{[}D^{g}_{\mu_{0}^{2}}(\zeta)\frac{x_{F}}{\zeta}{\cal
T}_{g}\left(\frac{x_{f}}{\zeta},k^{2};k^{2},\xi_{0}\right)+\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right){\cal
F}_{g}(\zeta,k^{2};k^{2};\xi_{0})\bigg{]}$ (112)
$\displaystyle\times\,\int_{q,t}s(t)\,s(p/\zeta+k-t)\,\bigg{[}1-\frac{k\cdot
q}{q^{2}}\bigg{]}\,s(p/\zeta+q-t).$

Finally, the sum of the real corrections reads

$\displaystyle\frac{d\sigma}{d^{2}pd\eta}\bigg{|}_{\rm NLO}^{\rm
real}=\frac{d\sigma}{d^{2}pd\eta}\bigg{|}_{\rm
NLO,r}^{q}+\frac{d\sigma}{d^{2}pd\eta}\bigg{|}_{\rm
NLO,r}^{g}+\sum_{i=1}^{4}({\rm Gen.\,NLO})_{i}\ .$ (113)

###  D.2 Virtual corrections at NLO

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}_{\rm NLO}^{\rm
virtual}=\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}_{\rm NLO,v}^{\rm q\to q\to
H}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}_{\rm NLO,v}^{\rm g\to g\to
H}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}_{\rm NLO,v}^{\rm g\to q\to
H}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}_{\rm NLO,v}^{\rm g\to\bar{q}\to H}\
.$ (114)

####  D.2.1 Virtual corrections to the quark production

This is the first term in eq. (114) which was discussed earlier in the paper.

####  D.2.2 Virtual corrections to the gluon production

The virtual corrections in the gluon channel given by the last three terms in
eq. (114). The last term is contribution from the antiquark, we will not
calculate that contribution explicitly. The two terms we need to focus on are
the second and the third terms in eq. (114). Let us first consider the second
term. This contribution can be divided into two pieces at partonic level. In
the large $N_{c}$ limit, we can write it as (eq. (4.10) in Ref. [37])

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm
NLO,v}=\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm
NLO,v-1}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm NLO,v-2}\ ,$
(115)

where

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm
NLO,v-1}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}N_{c}x_{p}f^{g}_{\mu_{0}^{2}}(x_{p})\int_{0}^{1}d\xi\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}\int_{y,\bar{y},z}e^{ip(y-\bar{y})}A^{i}_{\xi}(y-z)A^{i}_{\xi}(y-z)$
(116) $\displaystyle\times$
$\displaystyle\frac{1}{2}\bigg{\\{}2s_{A}(y,\bar{y})-s[y+\xi(y-z),z+\xi(y-z)]s[z+\xi(y-z),\bar{y}]s[\bar{y},y+\xi(y-z)]$
$\displaystyle\hskip
14.22636pt-s[z+\xi(y-z),y+\xi(y-z)]s[\bar{y},z+\xi(y-z)]s[y+\xi(y-z),\bar{y}]\Big{\\}}$

and

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm
NLO,v-2}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}N_{c}x_{p}f^{g}_{\mu_{0}^{2}}(x_{p})\int_{0}^{1}d\xi\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}\int_{y,\bar{y},z}e^{ip(y-\bar{y})}A^{i}_{\xi,x_{p}}(\bar{y}-z)A^{i}_{\xi,x_{p}}(\bar{y}-z)$
(117) $\displaystyle\times$
$\displaystyle\frac{1}{2}\bigg{\\{}2s_{A}(y,\bar{y})-s[\bar{y}+\xi(\bar{y}-z),z+\xi(\bar{y}-z)]s[z+\xi(\bar{y}-z),y]s[y,\bar{y}+\xi(\bar{y}-z)]$
$\displaystyle\hskip
14.22636pt-s[z+\xi(\bar{y}-z),\bar{y}+\xi(\bar{y}-z)]s[y,z+\xi(\bar{y}-z)]s[\bar{y}+\xi(\bar{y}-z),y]\Big{\\}}.$

Following the same procedure indicated previously, these two contributions can
be written as

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm
NLO,v-1}=-\frac{g^{2}}{(2\pi)^{3}}\,N_{c}\,x_{p}\,f^{g}_{\mu_{0}^{2}}(x_{p})\int_{\xi_{0}}^{1-\xi_{0}}d\xi\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}\frac{(-1)}{(2\pi)^{4}}\int_{y,\bar{y},z}\int_{l,m>\mu_{0}}\int_{k,q,t}\frac{l^{i}}{l^{2}}\frac{m^{i}}{m^{2}}$
$\displaystyle\hskip 17.07182pt\times
s(k)s(q)s(t)\frac{1}{2}\Big{\\{}2e^{iy(p-l-m+k-q)}e^{-i\bar{y}(p+k-q)}e^{iz(l+m)}-e^{iy[p-l-m+k+\xi
q-(1+\xi)t]}e^{-i\bar{y}(p+q-t)}e^{iz[l+m-k+(1-\xi)q+\xi t]}$
$\displaystyle\hskip 91.04872pt-\,e^{iy[p-l-m-k-\xi
q+(1+\xi)t]}e^{-i\bar{y}(p-q+t)}e^{iz[l+m+k-(1-\xi)q-\xi t]}\Big{\\}}$ (118)

and

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm
NLO,v-2}=-\frac{g^{2}}{(2\pi)^{3}}\,N_{c}\,x_{p}\,f^{g}_{\mu_{0}^{2}}(x_{p})\int_{\xi_{0}}^{1-\xi_{0}}d\xi\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}\frac{(-1)}{(2\pi)^{4}}\int_{y,\bar{y},z}\int_{l,m>\mu_{0}}\int_{k,q,t}\frac{l^{i}}{l^{2}}\frac{m^{i}}{m^{2}}$
$\displaystyle\hskip 17.07182pt\times
s(k)s(q)s(t)\frac{1}{2}\Big{\\{}2e^{iy(p+k-q)}e^{-i\bar{y}(p+l+m+k-q)}e^{iz(l+m)}-e^{iy(p-q+t)}e^{-i\bar{y}[p+l+m-k-\xi
q+(1+\xi)t]}e^{iz[l+m-k+(1-\xi)q+\xi t]}$ $\displaystyle\hskip
91.04872pt-\,e^{iy(p+q-t)}e^{-i\bar{y}[p+l+m+k+\xi
q-(1+\xi)t]}e^{iz[l+m+k-(1-\xi)q-\xi t]}\Big{\\}}.$ (119)

After performing the integrals, these two contributions can be combined with
some change of variables. The final result at partonic level reads

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g}_{\rm NLO,v}$
$\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}x_{p}f^{g}_{\mu_{0}^{2}}(x_{p})\int_{\xi_{0}}^{1-\xi_{0}}d\xi\,2N_{c}\,\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}$
(120)
$\displaystyle\times\,\int_{l>\mu_{0}}\int_{q,t}s(q)s(t-p)s(t)\bigg{\\{}\frac{1}{l^{2}}+\frac{l^{i}}{l^{2}}\frac{[(1-\xi)p-l-
q-t]^{i}}{[(1-\xi)p-l-q-t]^{2}}\bigg{\\}}.$

Upon introducing the fragmentation function, this contribution becomes

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g\to H}_{\rm
NLO,v}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{\xi_{0}}^{1-\xi_{0}}d\xi\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,2N_{c}\,\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}$
$\displaystyle\times\,\int_{l>\mu_{0}}\int_{q,t}s(q)s(t-p/\zeta)s(t)\bigg{\\{}\frac{1}{l^{2}}+\frac{l^{i}}{l^{2}}\frac{[(1-\xi)p/\zeta-
l-q-t]^{i}}{[(1-\xi)p/\zeta-l-q-t]^{2}}\bigg{\\}}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{\xi_{0}}^{1-\xi_{0}}d\xi\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,2N_{c}\,\int_{l>\mu_{0}}\int_{q,t}s(q)s(t-p/\zeta)s(t)$
$\displaystyle\hskip
5.69046pt\times\,\Bigg{\\{}\Bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}\bigg{\\{}\left[\frac{1}{l^{2}}+\frac{l^{i}}{l^{2}}\frac{[p/\zeta-
l-q-t]^{i}}{[p/\zeta-l-q-t]^{2}}\right]$ $\displaystyle\hskip
22.76228pt+\left[\frac{l^{i}}{l^{2}}\frac{[(1-\xi)p/\zeta-l-
q-t]^{i}}{[(1-\xi)p/\zeta-l-q-t]^{2}}-\frac{l^{i}}{l^{2}}\frac{[p/\zeta-l-
q-t]^{i}}{[p/\zeta-l-q-t]^{2}}\right]\bigg{\\}}\Bigg{\\}}.$

Let us now consider the third term in eq. (114) (the last term in eq. (114)
corresponds to the antiquark and we will not consider it explicitly). This
term can also be grouped into two contributions

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q}_{\rm
NLO,v}=\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q}_{\rm
NLO,v-1}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q}_{\rm NLO,v-2}\ ,$
(122)

and each of them are given by (eq. (4.11) in Ref. [37])

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q}_{\rm
NLO,v-1}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}x_{p}f^{g}_{\mu_{0}^{2}}(x_{p})\int_{\xi_{0}}d\xi\frac{1}{2}\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\int_{y,\bar{y},z}e^{ip(y-\bar{y})}A^{i}_{\xi}(y-z)A^{i}_{\xi}(y-z)$
(123)
$\displaystyle\times\Big{\\{}s_{A}(y,\bar{y})-s[y+\xi(y-z),\bar{y}]s[\bar{y},z+\xi(y-z)]\Big{\\}}$

and

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q}_{\rm
NLO,v-2}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}x_{p}f^{g}_{\mu_{0}^{2}}(x_{p})\int_{\xi_{0}}d\xi\frac{1}{2}\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\int_{y,\bar{y},z}e^{ip(y-\bar{y})}A^{i}_{\xi}(\bar{y}-z)A^{i}_{\xi}(\bar{y}-z)$
(124)
$\displaystyle\times\Big{\\{}s_{A}(y,\bar{y})-s[\bar{y}+\xi(\bar{y}-z),y]s[y,z+\xi(\bar{y}-z)]\Big{\\}}.$

Within our approximations, these two contributions can be written

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q}_{\rm
NLO,v-1}=-\frac{g^{2}}{(2\pi)^{3}}x_{p}f^{g}_{\mu_{0}^{2}}(x_{p})\int_{\xi_{0}}d\xi\,\frac{1}{2}\,\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\frac{(-1)}{(2\pi)^{4}}\int_{y,\bar{y},z}\int_{l,m>\mu_{0}}\int_{k,q}\frac{l^{i}}{l^{2}}\frac{m^{i}}{m^{2}}\,s(k)s(q)$
$\displaystyle\times\Big{\\{}e^{iy(p-l-m+k+q)}e^{-i\bar{y}(p+k-q)}e^{iz(l+m)}-e^{iy[p-l-m+(1+\xi)k-\xi
q]}e^{-i\bar{y}(p+k-q)}e^{iz[l+m-\xi k-(1-\xi)q]}\Big{\\}}$ (125)

and

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q}_{\rm
NLO,v-1}=-\frac{g^{2}}{(2\pi)^{3}}x_{p}f^{g}_{\mu_{0}^{2}}(x_{p})\int_{\xi_{0}}d\xi\,\frac{1}{2}\,\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\frac{(-1)}{(2\pi)^{4}}\int_{y,\bar{y},z}\int_{l,m>\mu_{0}}\int_{k,q}\frac{l^{i}}{l^{2}}\frac{m^{i}}{m^{2}}\,s(k)s(q)$
$\displaystyle\times\Big{\\{}e^{iy(p+k-q)}e^{-i\bar{y}(p+l+m+k-q)}e^{iz(l+m)}-e^{iy(p-k+q)}e^{-i\bar{y}[p+l+m-(1+\xi)k+\xi
q]}e^{iz[l+m-\xi k-(1-\xi)q]}\Big{\\}}.$ (126)

These two contributions can be combined after performing the integrals and the
final result at partonic level reads

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q}_{\rm NLO,v}$
$\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}x_{p}f^{g}_{\mu_{0}^{2}}(x_{p})\int_{\xi_{0}}d\xi\,\frac{1}{2}\,\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\int_{l>\mu_{0}}\int_{q}s(p-q)s(q)$
(127)
$\displaystyle\times\,\bigg{\\{}\frac{2}{l^{2}}+\frac{l^{i}}{l^{2}}\frac{[\xi
p-l-q]^{i}}{[\xi
p-l-q]^{2}}+\frac{l^{i}}{l^{2}}\frac{[(1-\xi)p-l-q]^{i}}{[(1-\xi)p-l-q]^{2}}\bigg{\\}}.$

Introducing the fragmentation function and reorganizing the terms, we get

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,v}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{\xi_{0}}^{1}\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,\frac{1}{2}\,\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\int_{l>\mu_{0}}\int_{q,t}s(t)s(t-p/\zeta)s(q)$
(128)
$\displaystyle\times\,\bigg{\\{}\frac{2}{l^{2}}+\frac{l^{i}}{l^{2}}\frac{[\xi
p/\zeta-l-t]^{i}}{[\xi p/\zeta-
l-t]^{2}}+\frac{l^{i}}{l^{2}}\frac{[(1-\xi)p/\zeta-l-t]^{i}}{[(1-\xi)p/\zeta-
l-t]^{2}}\bigg{\\}}.$

Note that eqs. (D.2.2) and (128) are the final expressions for the virtual
contributions in the gluon channel. In the rest of the discussion we will
follow the same arguments adopted in the analysis of the virtual contributions
in the quark channel. Namely, we will separate each of these two contributions
into a piece containing a large logarithm which will be included in the
evolution of the TMDs and a piece without any large logarithms which therefore
will be referred to as new genuine NLO corrections.

Let us start with eq. (D.2.2). The last term in this equation is finite - it
does not contain either transverse or longitudinal logs, and we set it aside.
Just like in eq. (III.3.2), we can perform the angular integration over the
angle of vector $l$ in the first tern in eq. (D.2.2),

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g\to H(1)}_{\rm
NLO,v}$ $\displaystyle=$
$\displaystyle-2\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{q,t}\int_{\mu_{0}^{2}}^{[q+t-p/\zeta]^{2}}\frac{d^{2}l}{l^{2}}\int_{\xi_{0}}^{1}d\xi\,\frac{x_{F}}{\zeta}\,f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,$
(129)
$\displaystyle\times\,N_{c}\,\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}s(q)s(t-p/\zeta)s(t).$

Again, following the same arguments in eq. (III.3.2) we divide the integration
region over $l$ into two pieces as

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g\to H(1)}_{\rm
NLO,v}$ $\displaystyle=$
$\displaystyle-2\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{q,t}\bigg{[}\int_{\mu_{0}^{2}}^{\mu^{2}}+\int_{\mu^{2}}^{[q+t-p/\zeta]^{2}}\bigg{]}\frac{d^{2}l}{l^{2}}\int_{\xi_{0}}^{1}d\xi\,\frac{x_{F}}{\zeta}\,f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,$
(130)
$\displaystyle\times\,N_{c}\,\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}s(q)s(t-p/\zeta)s(t).$

The first integral over $l$ in eq. (130) will give a large log which will
contribute to the evolution of the gluon TMD PDFs and TMD FFs, the second
integral does not give a large log and therefore it is a genuine NLO
correction. We write it as

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g\to H}_{\rm
NLO,v}$ $\displaystyle=$
$\displaystyle-2\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{t}\,\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d^{2}l}{l^{2}}\int_{\xi_{0}}^{1}d\xi\,\frac{x_{F}}{\zeta}\,f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,$
(131)
$\displaystyle\times\,N_{c}\,\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}s(t-p/\zeta)s(t)+({\rm
Gen.NLO})_{5}\ ,$

with

$\displaystyle({\rm Gen.NLO})_{5}$ $\displaystyle=$
$\displaystyle-2\frac{g^{2}N_{c}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{q,t}\,\int_{\mu^{2}}^{[q+t-(1-p/\zeta]^{2}}\bigg{]}\frac{d^{2}l}{l^{2}}\int_{\xi_{0}}^{1}d\xi\,\frac{x_{F}}{\zeta}\,f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,$
$\displaystyle\times\,\,\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}s(q)s(t-p/\zeta)s(t)$
$\displaystyle-$ $\displaystyle
2\frac{g^{2}N_{c}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{\xi_{0}}^{1}d\xi\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,\,\bigg{[}\frac{1-\xi}{\xi}+\frac{\xi}{1-\xi}+\xi(1-\xi)\bigg{]}$
$\displaystyle\times\,\int_{l>\mu_{0}}\int_{q,t}s(q)s(t-p/\zeta)s(t)\bigg{\\{}\left[\frac{l^{i}}{l^{2}}\frac{[(1-\xi)p/\zeta-
l-q-t]^{i}}{[(1-\xi)p/\zeta-l-q-t]^{2}}-\frac{l^{i}}{l^{2}}\frac{[p/\zeta-l-
q-t]^{i}}{[p/\zeta-l-q-t]^{2}}\right]\bigg{\\}}.$

Let us now consider eq. (128). We write it as

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,v}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{\xi_{0}}^{1}\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,\frac{1}{2}\,\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\int_{l>\mu_{0}}\int_{q,t}s(t)s(t-p/\zeta)s(q)$
(133)
$\displaystyle\times\,\bigg{\\{}\bigg{[}\frac{1}{l^{2}}+\frac{l^{i}}{l^{2}}\frac{[\xi
p/\zeta-l-t]^{i}}{[\xi p/\zeta-
l-t]^{2}}\bigg{]}+\bigg{[}\frac{1}{l^{2}}+\frac{l^{i}}{l^{2}}\frac{[(1-\xi)p/\zeta-
l-t]^{i}}{[(1-\xi)p/\zeta-l-t]^{2}}\bigg{]}\bigg{\\}}.$

We can now perform the angular integration over the angle of vector $l$,
separately for the first two terms and the last two terms. The result reads

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,v}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{t}\bigg{[}\int_{\mu_{0}^{2}}^{[t-(1-\xi)p/\zeta]^{2}}+\int_{\mu_{0}^{2}}^{[t-\xi
p/\zeta]^{2}}\bigg{]}\frac{d^{2}l}{l^{2}}\int_{\xi_{0}}^{1}\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,$
(134)
$\displaystyle\times\,\frac{1}{2}\,\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\,s(t)s(t-p/\zeta).$

Following the same arguments, we divide the integration region over $l$ into
two pieces:

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,v}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{t}\bigg{[}\int_{\mu_{0}^{2}}^{\mu^{2}}+\int_{\mu^{2}}^{[t-(1-\xi)p/\zeta]^{2}}+\int_{\mu_{0}^{2}}^{\mu^{2}}+\int_{\mu^{2}}^{[t-\xi
p/\zeta]^{2}}\bigg{]}\frac{d^{2}l}{l^{2}}\int_{\xi_{0}}^{1}\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,$
(135)
$\displaystyle\times\,\frac{1}{2}\,\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\,s(t)s(t-p/\zeta).$

The first and the third terms in eq. (135) contain large logarithms and they
contribute to the evolution of the TMDs while the second and fourth terms are
genuine NLO corrections. Thus the final expression can be written as

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,v}$ $\displaystyle=$
$\displaystyle-2\,\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{t}\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d^{2}l}{l^{2}}\int_{\xi_{0}}^{1}\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,\frac{1}{2}\,\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\,s(t)s(t-p/\zeta)$
(136) $\displaystyle+$ $\displaystyle({\rm Gen.NLO})_{6}\ ,$

where the genuine NLO correction reads

$\displaystyle({\rm Gen.NLO})_{6}$ $\displaystyle=$
$\displaystyle-\,\frac{g^{2}}{(2\pi)^{3}}\,S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}D^{g}_{\mu_{0}^{2}}(\zeta)\int_{t}\bigg{[}\int_{\mu^{2}}^{[t-(1-\xi)p/\zeta]^{2}}+\int_{\mu^{2}}^{[t-\xi
p/\zeta]^{2}}\bigg{]}\frac{d^{2}l}{l^{2}}\int_{\xi_{0}}^{1}\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\,$
(137)
$\displaystyle\times\,\frac{1}{2}\,\Big{[}\xi^{2}+(1-\xi)^{2}\Big{]}\,s(t)s(t-p/\zeta).$

###  D.3 Combining all contributions: gluon channel

In the gluon channel, the leading order expression reads

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to H}_{\rm LO}$
$\displaystyle=$ $\displaystyle
S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\,D^{g}_{\mu_{0}^{2}}\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\int_{t}\,s(t)s(t-p/\zeta).$
(138)

Combining the LO expression (eq. (138)) with the virtual corrections (eqs.
(131) and (136)) in the gluon channel, we get

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to H}_{\rm
LO}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g\to H}_{\rm
NLO,v}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,v}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to\bar{q}\to H}_{\rm
NLO,v}=S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\int_{t}\int_{0}^{\mu_{0}^{2}}d^{2}k$
$\displaystyle\hskip
51.21504pt\times\,\bigg{[}D^{g}_{\mu_{0}^{2}}(\zeta)\frac{x_{F}}{\zeta}{\cal
T}_{g}\left(\frac{x_{F}}{\zeta},k^{2};\mu^{2};\xi_{0}\right)+{\cal
F}_{g}(\zeta,k^{2};\mu^{2};\xi_{0})\frac{x_{F}}{\zeta}f^{g}_{\mu_{0}^{2}}\left(\frac{x_{F}}{\zeta}\right)\bigg{]}\,s(t)s(t-p/\zeta)$
$\displaystyle\hskip 51.21504pt+\;({\rm Gen.NLO})_{5}+({\rm Gen.NLO})_{6}\ .$
(139)

Here, we have used the evolution equations for the TMDs over the resolution
scale given in eqs. (17) and (22). Note that we have not calculated the last
term on the left hand side of eq. (D.3). We can rewrite eq. (D.3) in a more
convenient way (by following the discussions in the quark channel):

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to H}_{\rm
LO}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g\to H}_{\rm
NLO,v}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,v}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to\bar{q}\to H}_{\rm
NLO,v}\simeq
S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\int_{t,q}\int_{0}^{\mu_{0}^{2}}d^{2}k\int_{0}^{\mu_{0}^{2}}d^{2}l$
$\displaystyle\hskip 22.76228pt\times\,{\cal
F}_{g}(\zeta,l^{2};\mu^{2};\xi_{0})\frac{x_{F}}{\zeta}{\cal
T}_{g}\left(\frac{x_{F}}{\zeta},k^{2};\mu^{2};\xi_{0}\right)s(t)\,s\big{(}p/\zeta-t+(k+l)\big{)}\bigg{[}1-\frac{(k+l)\cdot
q}{q^{2}}\bigg{]}\,s(p/\zeta-t+q)$ $\displaystyle\hskip 22.76228pt+\;({\rm
Gen.NLO})_{5}+({\rm Gen.NLO})_{6}\ .$ (140)

Finally, combining this result with the real NLO correction given in eq.
(112), we get

$\displaystyle\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to H}_{\rm
LO}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g}_{\rm
NLO,r}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to g\to H}_{\rm
NLO,v}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to q\to H}_{\rm
NLO,v}+\frac{d\sigma}{d^{2}p\,d\eta}\bigg{|}^{\rm g\to\bar{q}\to H}_{\rm
NLO,v}\simeq S_{\perp}\int_{x_{F}}^{1}\frac{d\zeta}{\zeta^{2}}\int_{t,q}\int
d^{2}k\int d^{2}l$ $\displaystyle\hskip 22.76228pt\times\,{\cal
F}_{g}(\zeta,l^{2};\mu^{2};\xi_{0})\frac{x_{F}}{\zeta}{\cal
T}_{g}\left(\frac{x_{F}}{\zeta},k^{2};\mu^{2};\xi_{0}\right)s(t)\,s\big{(}p/\zeta-t+(k+l)\big{)}\bigg{[}1-\frac{(k+l)\cdot
q}{q^{2}}\bigg{]}\,s(p/\zeta-t+q)$ $\displaystyle\hskip 22.76228pt+\;({\rm
Gen.NLO})_{5}+({\rm Gen.NLO})_{6}\ .$ (141)

## References

  * [1] Y. V. Kovchegov and E. Levin, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 33 (2012), 1-350 Cambridge University Press, 2022, ISBN 978-1-00-929144-6, 978-1-00-929141-5, 978-1-00-929142-2, 978-0-521-11257-4, 978-1-139-55768-9 doi:10.1017/9781009291446
  * [2] M. Hentschinski, C. Royon, C. Royon, M. A. Peredo, C. Baldenegro, A. Bellora, R. Boussarie, F. G. Celiberto, S. Cerci and G. Chachamis, et al. Acta Phys. Polon. B 54, no.3, 2 (2023) doi:10.5506/APhysPolB.54.3-A2 [arXiv:2203.08129 [hep-ph]].
  * [3] F. Gelis, E. Iancu, J. Jalilian-Marian and R. Venugopalan, Ann. Rev. Nucl. Part. Sci. 60 (2010), 463-489 doi:10.1146/annurev.nucl.010909.083629 [arXiv:1002.0333 [hep-ph]].
  * [4] I. Arsene et al. [BRAHMS], Phys. Rev. Lett. 93 (2004), 242303 doi:10.1103/PhysRevLett.93.242303 [arXiv:nucl-ex/0403005 [nucl-ex]].
  * [5] J. Adams et al. [STAR], Phys. Rev. Lett. 97 (2006), 152302 doi:10.1103/PhysRevLett.97.152302 [arXiv:nucl-ex/0602011 [nucl-ex]].
  * [6] B. Abelev et al. [ALICE], Phys. Rev. Lett. 110 (2013) no.8, 082302 doi:10.1103/PhysRevLett.110.082302 [arXiv:1210.4520 [nucl-ex]].
  * [7] G. Aad et al. [ATLAS], Phys. Lett. B 763 (2016), 313-336 doi:10.1016/j.physletb.2016.10.053 [arXiv:1605.06436 [hep-ex]].
  * [8] R. Aaij et al. [LHCb], JHEP 01 (2022), 166 doi:10.1007/JHEP01(2022)166 [arXiv:2107.10090 [hep-ex]].
  * [9] R. Aaij et al. [LHCb], Phys. Rev. Lett. 128 (2022) no.14, 142004 doi:10.1103/PhysRevLett.128.142004 [arXiv:2108.13115 [hep-ex]].
  * [10] A. Dumitru, A. Hayashigaki and J. Jalilian-Marian, Nucl. Phys. A 765 (2006), 464-482 doi:10.1016/j.nuclphysa.2005.11.014 [arXiv:hep-ph/0506308 [hep-ph]].
  * [11] T. Altinoluk and A. Kovner, Phys. Rev. D 83, 105004 (2011) doi:10.1103/PhysRevD.83.105004 [arXiv:1102.5327 [hep-ph]].
  * [12] G. A. Chirilli, B. W. Xiao and F. Yuan, Phys. Rev. Lett. 108 (2012), 122301 doi:10.1103/PhysRevLett.108.122301 [arXiv:1112.1061 [hep-ph]].
  * [13] G. A. Chirilli, B. W. Xiao and F. Yuan, Phys. Rev. D 86 (2012), 054005 doi:10.1103/PhysRevD.86.054005 [arXiv:1203.6139 [hep-ph]].
  * [14] V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. 15 (1972), 438-450 IPTI-381-71.
  * [15] G. Altarelli and G. Parisi, Nucl. Phys. B 126 (1977), 298-318 doi:10.1016/0550-3213(77)90384-4
  * [16] Y. L. Dokshitzer, Sov. Phys. JETP 46 (1977), 641-653
  * [17] I. Balitsky, Nucl. Phys. B 463 (1996), 99-160 doi:10.1016/0550-3213(95)00638-9 [arXiv:hep-ph/9509348 [hep-ph]].
  * [18] J. Jalilian-Marian, A. Kovner and H. Weigert, Phys. Rev. D 59 (1998), 014015 doi:10.1103/PhysRevD.59.014015 [arXiv:hep-ph/9709432 [hep-ph]].
  * [19] A. Kovner, J. G. Milhano and H. Weigert, Phys. Rev. D 62 (2000), 114005 doi:10.1103/PhysRevD.62.114005 [arXiv:hep-ph/0004014 [hep-ph]].
  * [20] A. Kovner and J. G. Milhano, Phys. Rev. D 61 (2000), 014012 doi:10.1103/PhysRevD.61.014012 [arXiv:hep-ph/9904420 [hep-ph]].
  * [21] H. Weigert, Nucl. Phys. A 703 (2002), 823-860 doi:10.1016/S0375-9474(01)01668-2 [arXiv:hep-ph/0004044 [hep-ph]].
  * [22] E. Iancu, A. Leonidov and L. D. McLerran, Nucl. Phys. A 692 (2001), 583-645 doi:10.1016/S0375-9474(01)00642-X [arXiv:hep-ph/0011241 [hep-ph]].
  * [23] E. Ferreiro, E. Iancu, A. Leonidov and L. McLerran, Nucl. Phys. A 703 (2002), 489-538 doi:10.1016/S0375-9474(01)01329-X [arXiv:hep-ph/0109115 [hep-ph]].
  * [24] I. Balitsky and G. A. Chirilli, Phys. Rev. D 77 (2008), 014019 doi:10.1103/PhysRevD.77.014019 [arXiv:0710.4330 [hep-ph]].
  * [25] A. Kovner, M. Lublinsky and Y. Mulian, Phys. Rev. D 89 (2014) no.6, 061704 doi:10.1103/PhysRevD.89.061704 [arXiv:1310.0378 [hep-ph]].
  * [26] A. Kovner, M. Lublinsky and Y. Mulian, JHEP 08 (2014), 114 doi:10.1007/JHEP08(2014)114 [arXiv:1405.0418 [hep-ph]].
  * [27] M. Lublinsky and Y. Mulian, JHEP 05 (2017), 097 doi:10.1007/JHEP05(2017)097 [arXiv:1610.03453 [hep-ph]].
  * [28] I. Balitsky, Phys. Rev. Lett. 81 (1998), 2024-2027 doi:10.1103/PhysRevLett.81.2024 [arXiv:hep-ph/9807434 [hep-ph]].
  * [29] I. Balitsky, Phys. Rev. D 60 (1999), 014020 doi:10.1103/PhysRevD.60.014020 [arXiv:hep-ph/9812311 [hep-ph]].
  * [30] Y. V. Kovchegov, Phys. Rev. D 61 (2000), 074018 doi:10.1103/PhysRevD.61.074018 [arXiv:hep-ph/9905214 [hep-ph]].
  * [31] A. M. Stasto, B. W. Xiao and D. Zaslavsky, Phys. Rev. Lett. 112 (2014) no.1, 012302 doi:10.1103/PhysRevLett.112.012302 [arXiv:1307.4057 [hep-ph]].
  * [32] A. M. Staśto, B. W. Xiao, F. Yuan and D. Zaslavsky, Phys. Rev. D 90 (2014) no.1, 014047 doi:10.1103/PhysRevD.90.014047 [arXiv:1405.6311 [hep-ph]].
  * [33] D. Zaslavsky, [arXiv:1409.8259 [hep-ph]].
  * [34] Z. B. Kang, I. Vitev and H. Xing, Phys. Rev. Lett. 113 (2014), 062002 doi:10.1103/PhysRevLett.113.062002 [arXiv:1403.5221 [hep-ph]].
  * [35] B. W. Xiao and F. Yuan, [arXiv:1407.6314 [hep-ph]].
  * [36] E. Iancu, A. H. Mueller and D. N. Triantafyllopoulos, JHEP 12 (2016), 041 doi:10.1007/JHEP12(2016)041 [arXiv:1608.05293 [hep-ph]].
  * [37] T. Altinoluk, N. Armesto, G. Beuf, A. Kovner and M. Lublinsky, Phys. Rev. D 91, no.9, 094016 (2015) doi:10.1103/PhysRevD.91.094016 [arXiv:1411.2869 [hep-ph]].
  * [38] K. Watanabe, B. W. Xiao, F. Yuan and D. Zaslavsky, Phys. Rev. D 92, no.3, 034026 (2015) doi:10.1103/PhysRevD.92.034026 [arXiv:1505.05183 [hep-ph]].
  * [39] H. Y. Liu, Y. Q. Ma and K. T. Chao, Phys. Rev. D 100, no.7, 071503 (2019) doi:10.1103/PhysRevD.100.071503 [arXiv:1909.02370 [nucl-th]].
  * [40] H. Y. Liu, Z. B. Kang and X. Liu, Phys. Rev. D 102, no.5, 051502 (2020) doi:10.1103/PhysRevD.102.051502 [arXiv:2004.11990 [hep-ph]].
  * [41] B. W. Xiao and F. Yuan, Phys. Lett. B 788, 261-269 (2019) doi:10.1016/j.physletb.2018.10.067 [arXiv:1806.03522 [hep-ph]].
  * [42] Y. Shi, L. Wang, S. Y. Wei and B. W. Xiao, Phys. Rev. Lett. 128 (2022) no.20, 202302 doi:10.1103/PhysRevLett.128.202302 [arXiv:2112.06975 [hep-ph]].
  * [43] H. y. Liu, K. Xie, Z. Kang and X. Liu, JHEP 07 (2022), 041 doi:10.1007/JHEP07(2022)041 [arXiv:2204.03026 [hep-ph]].
  * [44] L. Wang, L. Chen, Z. Gao, Y. Shi, S. Y. Wei and B. W. Xiao, Phys. Rev. D 107, no.1, 016016 (2023) doi:10.1103/PhysRevD.107.016016 [arXiv:2211.08322 [hep-ph]].
  * [45] E. Iancu and Y. Mulian, Nucl. Phys. A 985 (2019), 66-127 doi:10.1016/j.nuclphysa.2019.02.003 [arXiv:1809.05526 [hep-ph]].
  * [46] E. Iancu and Y. Mulian, JHEP 03 (2021), 005 doi:10.1007/JHEP03(2021)005 [arXiv:2009.11930 [hep-ph]].
  * [47] P. Caucal, F. Salazar and R. Venugopalan, JHEP 11 (2021), 222 doi:10.1007/JHEP11(2021)222 [arXiv:2108.06347 [hep-ph]].
  * [48] P. Taels, T. Altinoluk, G. Beuf and C. Marquet, JHEP 10 (2022), 184 doi:10.1007/JHEP10(2022)184 [arXiv:2204.11650 [hep-ph]].
  * [49] P. Caucal, F. Salazar, B. Schenke and R. Venugopalan, JHEP 11 (2022), 169 doi:10.1007/JHEP11(2022)169 [arXiv:2208.13872 [hep-ph]].
  * [50] P. Caucal, F. Salazar, B. Schenke, T. Stebel and R. Venugopalan, [arXiv:2304.03304 [hep-ph]].
  * [51] F. Bergabo and J. Jalilian-Marian, JHEP 01 (2023), 095 doi:10.1007/JHEP01(2023)095 [arXiv:2210.03208 [hep-ph]].
  * [52] F. Bergabo and J. Jalilian-Marian, Phys. Rev. D 106 (2022) no.5, 054035 doi:10.1103/PhysRevD.106.054035 [arXiv:2207.03606 [hep-ph]].
  * [53] E. Iancu and Y. Mulian, [arXiv:2211.04837 [hep-ph]].
  * [54] F. Bergabo and J. Jalilian-Marian, Phys. Rev. D 107 (2023) no.5, 054036 doi:10.1103/PhysRevD.107.054036 [arXiv:2301.03117 [hep-ph]].
  * [55] K. J. Golec-Biernat and M. Wusthoff, Phys. Rev. D 59 (1998), 014017 doi:10.1103/PhysRevD.59.014017 [arXiv:hep-ph/9807513 [hep-ph]].
  * [56] J. Collins, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 32 (2011), 1-624 Cambridge University Press, 2013, ISBN 978-1-107-64525-7, 978-1-107-64525-7, 978-0-521-85533-4, 978-1-139-09782-6
  * [57] R. Boussarie, M. Burkardt, M. Constantinou, W. Detmold, M. Ebert, M. Engelhardt, S. Fleming, L. Gamberg, X. Ji and Z. B. Kang, et al. [arXiv:2304.03302 [hep-ph]].
  * [58] F. Hautmann, H. Jung, A. Lelek, V. Radescu and R. Zlebcik, JHEP 01, 070 (2018) doi:10.1007/JHEP01(2018)070 [arXiv:1708.03279 [hep-ph]].
  * [59] A. B. Martinez, [arXiv:2307.06704 [hep-ph]].
  * [60] J. C. Collins and D. E. Soper, Nucl. Phys. B 194, 445-492 (1982) doi:10.1016/0550-3213(82)90021-9
  * [61] J. C. Collins, D. E. Soper and G. F. Sterman, Nucl. Phys. B 250, 199-224 (1985) doi:10.1016/0550-3213(85)90479-1
  * [62] J. w. Qiu and X. f. Zhang, Phys. Rev. D 63, 114011 (2001) doi:10.1103/PhysRevD.63.114011 [arXiv:hep-ph/0012348 [hep-ph]].
  * [63] E. L. Berger and J. w. Qiu, Phys. Rev. D 67, 034026 (2003) doi:10.1103/PhysRevD.67.034026 [arXiv:hep-ph/0210135 [hep-ph]].
  * [64] A. Bacchetta and A. Prokudin, Nucl. Phys. B 875, 536-551 (2013) doi:10.1016/j.nuclphysb.2013.07.013 [arXiv:1303.2129 [hep-ph]].
  * [65] M. A. Ebert, J. K. L. Michel, I. W. Stewart and Z. Sun, JHEP 07, 129 (2022) doi:10.1007/JHEP07(2022)129 [arXiv:2201.07237 [hep-ph]].
  * [66] J. C. Collins and D. E. Soper, Nucl. Phys. B 193, 381 (1981) [erratum: Nucl. Phys. B 213, 545 (1983)] doi:10.1016/0550-3213(81)90339-4
  * [67] S. M. Aybat and T. C. Rogers, Phys. Rev. D 83, 114042 (2011) doi:10.1103/PhysRevD.83.114042 [arXiv:1101.5057 [hep-ph]].
  * [68] T. Becher and M. Neubert, Eur. Phys. J. C 71, 1665 (2011) doi:10.1140/epjc/s10052-011-1665-7 [arXiv:1007.4005 [hep-ph]].
  * [69] T. Becher, M. Neubert and D. Wilhelm, JHEP 02, 124 (2012) doi:10.1007/JHEP02(2012)124 [arXiv:1109.6027 [hep-ph]].
  * [70] M. G. Echevarria, A. Idilbi and I. Scimemi, JHEP 07, 002 (2012) doi:10.1007/JHEP07(2012)002 [arXiv:1111.4996 [hep-ph]].
  * [71] J. Y. Chiu, A. Jain, D. Neill and I. Z. Rothstein, JHEP 05, 084 (2012) doi:10.1007/JHEP05(2012)084 [arXiv:1202.0814 [hep-ph]].
  * [72] M. G. Echevarría, A. Idilbi and I. Scimemi, Phys. Lett. B 726, 795-801 (2013) doi:10.1016/j.physletb.2013.09.003 [arXiv:1211.1947 [hep-ph]].
  * [73] T. Becher, M. Neubert and D. Wilhelm, JHEP 05, 110 (2013) doi:10.1007/JHEP05(2013)110 [arXiv:1212.2621 [hep-ph]].
  * [74] M. G. Echevarria, A. Idilbi and I. Scimemi, Phys. Rev. D 90, no.1, 014003 (2014) doi:10.1103/PhysRevD.90.014003 [arXiv:1402.0869 [hep-ph]].
  * [75] M. A. Ebert, I. W. Stewart and Y. Zhao, JHEP 03, 099 (2020) doi:10.1007/JHEP03(2020)099 [arXiv:1910.08569 [hep-ph]].
  * [76] J. Jalilian-Marian and A. H. Rezaeian, Phys. Rev. D 85, 014017 (2012) doi:10.1103/PhysRevD.85.014017 [arXiv:1110.2810 [hep-ph]].
  * [77] E. A. Kuraev, L. N. Lipatov and V. S. Fadin, Sov. Phys. JETP 45 (1977), 199-204
  * [78] I. I. Balitsky and L. N. Lipatov, Sov. J. Nucl. Phys. 28 (1978), 822-829