# An Empirical Study on Model-agnostic Debiasing Strategies
for Robust Natural Language Inference

Tianyu Liu1 Xin Zheng3 11footnotemark: 1 Xiaoan Ding4 Baobao Chang1 2 Zhifang
Sui1 2
1 Peking University, Beijing, China 2 Peng Cheng Laboratory, Shenzhen, China
3 Beijing University of Posts and Telecommunications, Beijing, China
4 University of Chicago, IL, USA
{tianyu0421, chbb<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>Equal contribution.

###### Abstract

The prior work on natural language inference (NLI) debiasing mainly targets at
one or few known biases while not necessarily making the models more robust.
In this paper, we focus on the model-agnostic debiasing strategies and explore
how to (or is it possible to) make the NLI models robust to multiple distinct
adversarial attacks while keeping or even strengthening the models’
generalization power. We firstly benchmark prevailing neural NLI models
including pretrained ones on various adversarial datasets. We then try to
combat distinct known biases by modifying a mixture of experts (MoE) ensemble
method Clark et al. (2019) and show that it’s nontrivial to mitigate multiple
NLI biases at the same time, and that model-level ensemble method outperforms
MoE ensemble method. We also perform data augmentation including text swap,
word substitution and paraphrase and prove its efficiency in combating various
(though not all) adversarial attacks at the same time. Finally, we investigate
several methods to merge heterogeneous training data (1.35M) and perform model
ensembling, which are straightforward but effective to strengthen NLI models.

##  1 Introduction

Natural language inference (NLI) (also known as recognizing textual
entailment) is a widely studied task which aims to infer the relationship
(e.g., _entailment_ , _contradiction_ , _neutral_) between two fragments of
text, known as _premise_ and _hypothesis_ Dagan et al. (2005, 2013). Recent
works have found that NLI models are sensitive to the compositional features
Nie et al. (2019), syntactic heuristics McCoy et al. (2019), stress test
Geiger et al. (2018); Naik et al. (2018) and human artifacts in the data
collection phase Gururangan et al. (2018); Poliak et al. (2018b); Tsuchiya
(2018). Accordingly, several adversarial datasets are proposed for these known
biases111In this paper, we use the term ‘bias’ to refer to these known dataset
biases in NLI following Clark et al. (2019). In other context, ‘bias’ may
refer to systematic mishandling of gender or evidences of racial stereotypes
Rudinger et al. (2017) in NLI datasets or models..

Through our preliminary trials on specific adversarial datasets, we find that
although the model specific or dataset specific debiasing methods could
increase the model performance on the paired adversarial dataset, they might
hinder the model performance on other adversarial datasets, as well as hurt
the model generalization power, i.e. deficient scores on cross-datasets or
cross-domain settings. These phenomena motivate us to investigate if it exists
a unified model-agnostic debiasing strategy which can mitigate distinct (or
even all) known biases while keeping or strengthening the model generalization
power.

We begin with NLI debiasing models. To make our trials more generic, we adopt
a mixture of experts (MoE) strategy Clark et al. (2019), which is known for
being model-agnostic and is adaptable to various kinds of known biases, as
backbone. Specifically we treat three known biases, namely word overlap,
length mismatch and partial input heuristics as independent experts and train
corresponding debiasing models. Our results show that the debiasing methods
tied to one particular known bias may not be sufficient to build a
generalized, robust model. This motivates us to investigate a better solution
to integrate the advantages of distinct debiasing models. We find model-level
ensemble is more effective than other MoE ensemble methods. Although our
findings are based on the MoE backbone due to the prohibitive exhaustive
studies on the all existing debiasing strategies, we provide actionable
insights on combining distinct NLI debiasing methods to the practitioners.

Then we explore model agnostic and generic data augmentation methods in NLI,
including text swap, word substitution and paraphrase. We find these methods
could help NLI models combat multiple (though not all) adversarial attacks,
e.g. augmenting training data by swapping hypothesis and premise could boost
the model performance on stress tests and lexical inference test, and data
augmentation by paraphrasing the hypothesis sentences could help the models
resist the superficial patterns from syntactic and partial input heuristics.
We also observe that increasing training size by incorporating heterogeneous
training resources is a simple but effective method to build robust and
generalized models. Specifically we investigate how to incorporate different
training data with different sizes and annotation processes, as well as the
best way to perform model ensembling.

Datasets | Paper | Categories | Labels | Size
---|---|---|---|---
PI-CD | (a) | $1\|3\|7$ | (E,N,C) | 3.2k
PI-SP | (b) | $1\|3\|7$ | (E,N,C) | .37k
IS-SD | (c) | $2\|5\|8$ | ($\neg$E, E) | 30k
IS-CS | (d) | $2\|3\|7$ | (E,N,C) | .65k
LI-LI | (e)(f) | $2\|4\|9$ | (E,C) | 9.9K
LI-TS | (g)(h) | $2\|6\|10$ | ($\neg$C, C) | 9.8K
ST-WO | (e) | $2\|4\|11$ | (E,N,C) | 9.8K
ST-NE | (e) | $2\|4\|11$ | (E,N,C) | 9.8K
ST-LM | (e) | $2\|4\|11$ | (E,N,C) | 9.8K
ST-SE | (e) | $2\|4\|12$ | (E,N,C) | 31K

(a) Gururangan et al. (2018) | (b) Liu et al. (2020)
---|---
(c) McCoy et al. (2019) | (d) Nie et al. (2019)
(e) Naik et al. (2018) | (f) Glockner et al. (2018)
(g) Wang et al. (2019c) | (h) Minervini and Riedel

Category | First-level | Second-level
---|---|---
1 | $(\mathbf{I})$ | Partial input heuristics
2 | $(\mathbf{I})$ | Inter-sentence heuristics
3 | $(\mathbf{II})$ | Instance selection
4 | $(\mathbf{II})$ | Single Sentence Modification
5 | $(\mathbf{II})$ | Sentence Pair Modification
6 | $(\mathbf{II})$ | Sentence Pair Swapping
7 | $(\mathbf{III})$ | Lexical Statistical Irregularity
8 | $(\mathbf{III})$ | Syntactic Statistical Irregularity
9 | $(\mathbf{III})$ | Lexical Inference
10 | $(\mathbf{III})$ | First Order Logic
11 | $(\mathbf{III})$ | Stress Test - Distraction Test
12 | $(\mathbf{III})$ | Stress Test - Noise Test

$(\mathbf{I})$ Where are the heuristics?
$(\mathbf{II})$ How did the dataset constructed?
$(\mathbf{III})$ Which aspect did the dataset detect?

Table 1: The information of adversarial datasets (Sec 2) we use in this paper.
We categorize and rename these datasets as discussed in Sec 2.1.

##  2 Benchmark Datasets

Our benchmark datasets include the adversarial datasets222Some datasets listed
in Table 1 were originally proposed to probe for systematicity. Here we call
them ‘adversarial’ datasets in the sense that the NLI models can not reach the
same performance on these datasets as the in-domain test sets. and some widely
used general-purpose NLI datasets which test the generalization power of NLI
models. 333The datasets used in this paper can be found in the following
github repository https://github.com/tyliupku/nli-debiasing-datasets

###  2.1 Adversarial Datasets

Categorization: to provide more insights on how the adversarial datasets
attack the models, we roughly categorize them in Table 1 according to their
characteristics and elaborate the categorization in this section. To
facilitate the narrative of following sections, we rename the adversarial
datasets according to their prominent features.

Comparability: all the following datasets are collected based on the public
available resources proposed by their authors, thus the experimental results
in this paper are comparable to the numbers reported in the original papers
and the other papers that use these datasets444The ownership of these datasets
belong to their authors. We encourage the readers to acknowledge and cite the
original papers listed in Table 1 when using them..

####  2.1.1 Partial-input (PI) Heuristics

Partial-input heuristics refer to the hypothesis-only bias Poliak et al.
(2018b) in NLI.

Classifier Detected Datasets (PI-CD): Gururangan et al. (2018) trained a
neural classifier (fastText555https://fasttext.cc/) on the hypothesis
sentences and then treated those instances in the SNLI test sets which can not
be correctly classified as ‘hard’ instances.

Surface Pattern Datasets (PI-SP): Liu et al. (2020) recognized surface
patterns which are highly correlated to the specific labels and
correspondingly proposed adversarial test sets which are against surface
patterns’ indications. We use their ‘hard’ instances for MultiNLI mismatched
dev set as adversarial datasets.

####  2.1.2 Inter-sentences (IS) Heuristics

Syntactic Diagnostic Datasets (IS-SD): The HANS dataset McCoy et al. (2019)
includes lexical overlap, subsequence and constituent heuristics between the
hypothesis and premises sentences, e.g. the model might incorrectly predict ‘
_entailment_ ’ for instance like ‘The actor was paid by the judge’ and ‘The
actor paid the judge’.

Compositionality-sensitivity Datasets (IS-CS): Nie et al. (2019) trained a
softmax regression model using unigram pattern pair features across two
sentences as well as unigram features in hypothesis and premise sentences to
obtain the ’lexically misleading scores (LMS)’ for each instance in the test
sets. We use $\mathbf{CS}_{0.7}$ in their paper which denotes the subsets
whose LMS are larger that 0.7.

####  2.1.3 Logical Inference Ability (LI)

Lexical Inference Test (LI-LI): A proper NLI system should recognize hypernyms
and hyponyms; synonym and antonyms. We merge the “antonym” category in Naik et
al. (2018) and Glockner et al. (2018) to assess the models’ capability to
model lexical inference.

Text-fragment Swap Test (LI-TS): NLI system should also follow the first-order
logic constraints Wang et al. (2019c); Minervini and Riedel (2018). For
example, if the premise sentence $s_{p}$ entails the hypothesis sentence
$s_{h}$, then $s_{h}$ must not be contradicted by $s_{p}$. We then swap the
two sentences in the original MultiNLI mismatched dev sets. If the gold label
is ‘ _contradiction_ ’, the corresponding label in the swapped instance
remains unchanged, otherwise it becomes ‘ _non-contradicted_ ’.

####  2.1.4 Stress Test (ST)

We also include the “word overlap” (ST-WO), “negation” (ST-NE), “length
mismatch” (ST-LM) and “spelling errors” (ST-SE) in Naik et al. (2018), in
which ST-WO aims at detecting lexical overlap heuristics described in McCoy et
al. (2019) (IS-SD in Sec 2.1.2); ST-NE aims at detecting strong negative
lexical cues in partial-input sentences like PI-SP in Sec 2.1.2.

Figure 1: The surrogate correlations between different adversarial datasets.
We show the Pearson’s correlation coefficients of model performance on
different adversarial datasets in different runs (Sec 2.1.5).

####  2.1.5 Insights within Adversarial Tests

To provide actionable insights to NLP practitioners, we list how these
adversarial instances constructed and why they might fail NLI models in Table
1. Those adversarial datasets are potentially correlated with each other due
to similar constructing process or constructing goals. For example, ‘PI-CD’,
‘PI-SP’ and ‘IS-CS’ are all created with instance selection from original test
sets in order to attack the models which improperly rely on the superficial
lexical patterns, thus they might be potentially correlated. Although we could
analytically assess the correlation between adversarial datasets, it is hard
to demonstrate their underlying relationships from a quantitative perspective.
We instead try to utilize the model performances on these adversarial datasets
as surrogates to visualize their correlations. Concretely, we first collect
the model accuracy scores on each adversarial dataset according to 30 runs of
10 baseline models (3 runs each) listed in Table 3. Then we show the pearson
correlation coefficients of the model scores on any two distinct adversarial
datasets in Fig 1. According to Fig 1, ‘IS-SD’ (HANS) has higher correlation
with ‘IS-CS’ and ‘LI-TS’ compared with other adversarial datasets, we assume
this is because they are constructed based on cross sentence heuristics in the
natural occurring settings, as opposed to stress test datasets which add
tautology like ‘and true is true’ to the end of hypothesis sentences Naik et
al. (2018). ‘LI-LI’ instances are created by few lexical changes on premise
sentence which would easily fall into ‘word overlap’ heuristics as elaborated
in the ‘IS-SD’ dataset, thus ‘LI-LI’ has low correlation with ‘IS-SD’.

| SNLI | MNLI | DNLI | ANLI
---|---|---|---|---
Train | 549362 | 392702 | 249947 | 162765
Valid | 9842 | 9832 | 31696 | 2200
Test | 9824 | 9815 | 31232 | 2200
Table 2: Statistics for datasets used in Sec 5. For MNLI, we utlize the matched dev and mismatched dev sets as valid and test sets respectively. | Adversarial Test | Generalization Power Test |
---|---|---|---
| PI-CD | PI-SP | IS-SD | IS-CS | LI-LI | LI-TS | ST | Avg. | RTE | DIAG | SICK | SciTail | Avg. | MNLI
InferSent | 52.1 | 55.3 | 53.9 | 33.5 | 43.6 | 70.5 | 53.3 | 51.7 | 61.8 | 10.6 | 25.4 | 24.7 | 30.6 | 70.5
+ELMO | 48.6 | 59.8 | 55.2 | 42.1 | 38.5 | 72.4 | 52.7 | 52.8 | 62.5 | 9.8 | 24.6 | 18.5 | 28.9 | 72.5
DAM | 55.0 | 54.4 | 50.2 | 35.7 | 62.7 | 74.3 | 53.0 | 55.0 | 62.7 | 10.3 | 27.0 | 30.0 | 32.5 | 70.3
ESIM | 55.1 | 66.3 | 49.8 | 52.7 | 63.2 | 79.6 | 53.8 | 60.1 | 66.2 | 11.3 | 25.1 | 27.5 | 32.5 | 77.3
$\mathrm{BERT}_{B}$ | 72.2 | 73.9 | 63.8 | 65.4 | 85.6 | 82.6 | 63.5 | 72.4 | 75.4 | 36.2 | 54.2 | 66.1 | 58.0 | 83.5
$\mathrm{BERT}_{L}$ | 74.7 | 75.5 | 70.4 | 70.6 | 87.9 | 83.8 | 67.3 | 75.7 | 77.6 | 39.4 | 55.5 | 68.3 | 60.2 | 85.7
$\mathrm{XLNet}_{B}$ | 73.1 | 77.9 | 71.2 | 70.4 | 85.5 | 84.8 | 68.5 | 75.9 | 78.0 | 39.2 | 55.8 | 66.7 | 59.9 | 86.6
$\mathrm{XLNet}_{L}$ | 78.8 | 81.7 | 76.7 | 77.3 | 93.4 | 88.5 | 72.4 | 81.3 | 83.4 | 45.9 | 57.6 | 73.0 | 65.0 | 89.3
$\mathrm{RoBERTa}_{B}$ | 76.6 | 80.9 | 72.0 | 74.1 | 89.6 | 85.3 | 66.4 | 77.8 | 80.9 | 42.1 | 55.9 | 69.0 | 62.0 | 87.4
$\mathrm{RoBERTa}_{L}$ | 80.0 | 79.2 | 80.0 | 77.0 | 92.4 | 88.6 | 73.4 | 81.5 | 84.4 | 50.5 | 57.3 | 72.2 | 66.1 | 89.9
Table 3: The performance of models on adversarial and generalization power
tests (Sec 2) trained on MultiNLI. B and L in the subscript denote base and
large versions of pretrained models. We use bold and underlined numbers to
represent the highest scores in each column/block. Same marks are also used in
Table 4, 5 and 6.

###  2.2 Other Data Resources

Generalization Power Test: we test the models on several general purpose
datasets, including NLI diagnostic dataset (Diag) Wang et al. (2019b), for
which we use ‘Matthews correlation coefficient’ Matthews (1975) as the
evaluation metric. We also incorporate RTE Dagan et al. (2005), SICK Marelli
et al. (2014) and SciTail Khot et al. (2018) in our testing.

Training Resources: apart from SNLI Bowman et al. (2015), and MultiNLI
Williams et al. (2018), we also incorporate Diverse NLI (DNLI) Poliak et al.
(2018a) and Adversarial NLI (ANLI) Nie et al. (2020) datasets for training.
For DNLI, we merge the subsets to form unified train/valid/test sets. Dataset
Statistics are shown in Table 2.

###  2.3 Model Performance on the Benchmark

We show the performance of different models trained on MultiNLI in Table 3.
The general trend is that more powerful model which has higher performance on
the original (in-domain) test sets (RoBERTa (large)) outperforms most models
in both adversarial and general purpose settings.

In the following sections, we investigate several model agnostic methods for
debiasing NLI models. Specifically, we are interested in: 1) how to (or is it
possible to) make the NLI models robust to multiple distinct adversarial
attacks using a unified debiasing method and 2) how the debiasing methods
influence model generalization power of NLI.

| Baseline | Word Overlap | Partial Input | Sentence Length | Debiasing Combination
---|---|---|---|---|---
| ($\mathrm{BERT}_{base}$) | ReW | BiasProd | ReW | BiasProd | ReW | BiasProd | MixW | AddProd | BestEn
PI-CD | 72.2 | 70.9 | 71.4 | 72.6 | 71.8 | 72.6 | 72.3 | 71.9 | 71.3 | 72.6
PI-SP | 73.9 | 70.6 | 70.1 | 74.7 | 73.0 | 75.2 | 73.3 | 71.7 | 70.4 | 73.9
IS-SD | 63.8 | 69.2 | 71.0 | 65.7 | 63.8 | 56.9 | 59.5 | 54.6 | 61.5 | 72.5
IS-CS | 65.4 | 64.8 | 64.2 | 67.1 | 68.9 | 64.9 | 66.9 | 65.4 | 68.9 | 64.9
LI-LI | 85.6 | 87.0 | 87.8 | 86.0 | 85.0 | 85.7 | 85.5 | 86.8 | 88.4 | 87.7
LI-TS | 82.6 | 81.8 | 81.7 | 82.0 | 82.3 | 81.3 | 83.7 | 82.3 | 81.9 | 84.5
ST-LM | 82.2 | 82.3 | 81.7 | 81.6 | 81.1 | 82.6 | 82.7 | 82.6 | 79.9 | 83.1
Gen. Avg. | 58.0 | 56.8 | 56.6 | 57.5 | 56.7 | 57.9 | 57.5 | 57.1 | 55.9 | 58.1
MNLI | 83.5 | 84.2 | 82.8 | 84.3 | 83.3 | 80.3 | 80.9 | 84.0 | 81.2 | 84.5
Table 4: The performance of debiasing methods (Sec 3) based on BERT base model
(baseline) trained on MultiNLI. ReW, BiasProd refer to instance reweighting
and bias product ensemble methods in Sec 3.1. Word overlap, partial input and
sentence length are the known biases in NLI (Sec 3.2). MixW, AddProd, BestEn
are our trials to combine distinct debiasing methods (Sec 3.3). ‘Gen. Avg’ is
the average score of test sets in generalization power test. Bold numbers mark
the highest score (besting debiasing model) _in each row_.

##  3 Mixture of Experts (MoE) Debiasing

We utilize the MoE ensemble model Clark et al. (2019) as the backbone to
mitigate three known biases in NLI. Concretely, we implement the ‘instance
reweighting’ and ‘bias product’ methods in Clark et al. (2019). Based on these
methods, we perform several trials on combating several distinct NLI biases at
the same time.

###  3.1 Debiasing Methods

Notations: for a known NLI bias, they firstly train a bias-only model $B$ and
then use its output $\mathbf{b}$ as a guidance to train the prime model. In
the context of three-way NLI training, $\mathbf{b_{i}}$ is a normalized
3-element vector which represents the predicted possibility of each NLI label
for i-th training example. Suppose $\mathbf{p_{i}}$ is output of the prime
model which has the same meaning as $\mathbf{b_{i}}$.

Instance Reweighting: suppose $b_{i}^{y_{i}}$ is the possibility that the
bias-only model assigns to the correct label $y_{i}$ for $i$-th training
example. They trained the models in a weighted version of the data, where the
weight $\alpha_{i}$ for the $i$-th training example is (1-$b_{i}^{y_{i}}$).
The loss function for a training batch with $k$ examples is a weighted sum of
instance-level loss $l_{i}$:
$\mathbf{L}_{batch}=\alpha_{i}*l_{i}/(\sum_{i=1}^{k}\alpha_{i})$.

Bias Product Ensemble: an ensemble method that is a product of experts
$\mathbf{\hat{p_{i}}}=softmax(log(\mathbf{p_{i}})+log(\mathbf{b_{i}}))$.

By doing so, the prime model would be encouraged to learn all the information
except the specific bias. An intuitive justification from the probabilistic
view can be found in Clark et al. (2019). Note that while training, only the
prime model is updated while the bias-only model remains unchanged.

###  3.2 Known Biases in NLI

Word overlap heuristics: To combat the word overlap heuristics (HANS McCoy et
al. (2019), renamed as IS-SD in Sec 2.1.2), Clark et al. (2019) used the
following features to train a bias-only model: (1) whether the hypothesis is a
sub-sequence of the premise, (2) whether all words in the hypothesis appear in
the premise, (3) the percent of words from the hypothesis that appear in the
premise, (4) the average and the max of the minimum distance between each
premise word with each hypothesis word. We use their trained bias-only model
output for experiments.

Partial input heuristics: To combat the hypothesis-only bias in NLI (PI-CD and
PI-SP in Sec 2.1.1), we use RoBERTa (base) model to train a bias-only model by
taking only hypothesis sentences as inputs. Our hypothesis-only model gets
60.4% accuracy on the mismatched dev set of MultiNLI, which is higher than the
reported numbers in Gururangan et al. (2018) (52.3%) and Poliak et al. (2018b)
(55.18%).

Sentence length heuristics: Gururangan et al. (2018) shows that the length of
hypothesis and premise over different labels is not evenly distributed (ST-LM
in Sec 2.1.4). So we trained a bias-only classifier based on the following
sentence length related features: 1) the sentence lengths of hypothesis and
premise sentences, 2) the mean and difference of these lengths. Our classifier
achieves 41.3% accuracy on the mismatched dev set of MultiNLI, which
outperforms the majority class baseline by 6.1%.

| Adversarial Test | Generalization Power Test |
---|---|---|---
| PI-CD | PI-SP | IS-SD | IS-CS | LI-LI | LI-TS | ST | Avg. | RTE | DIAG | SICK | SciTail | Avg. | MNLI
Baseline | 72.2 | 73.9 | 63.8 | 65.4 | 85.6 | 82.6 | 63.5 | 72.4 | 75.4 | 36.2 | 54.2 | 66.1 | 58.0 | 83.5
Text Swap | 71.7 | 72.8 | 63.5 | 67.4 | 86.3 | 86.8 | 66.5 | 73.6 | 73.3 | 35.3 | 54.7 | 66.8 | 57.6 | 83.7
Sub (synonym) | 69.8 | 72.0 | 62.4 | 65.8 | 85.2 | 82.8 | 64.3 | 71.8 | 74.4 | 34.2 | 55.1 | 65.8 | 57.4 | 83.5
Sub (MLM) | 71.0 | 72.8 | 64.4 | 65.9 | 85.6 | 83.3 | 64.9 | 72.6 | 74.8 | 34.7 | 55.4 | 65.7 | 57.7 | 83.6
Paraphrase | 72.1 | 74.6 | 66.5 | 66.4 | 85.7 | 83.1 | 64.8 | 73.3 | 75.8 | 35.1 | 55.0 | 65.0 | 57.7 | 83.7
Table 5: The performance of BERT base model under different data augmentation
strategies (Sec 4).

###  3.3 Combating Distinct Biases

Suppose we already have $m$ bias-only models $\\{B^{1},B^{2},\cdots,B^{m}\\}$
and the corresponding output
$\\{\mathbf{b^{1}},\mathbf{b^{2}},\cdots,\mathbf{b^{m}}\\}$ at hand, we test
three different approaches to integrate these models.

MixWeight: Using the product of weights from different debiasing models while
performing instance reweighting. We replace the weight for the $i$-th training
example ($\alpha_{i}$ in Sec 3.1) with $\prod_{j=1}^{m}(1-b_{i}^{y_{i}})$ and
utilize the same loss function as ‘instance reweighting’ in Sec 3.1).

AddProduct: We view different bias-only models as multiple independent experts
and then apply the bias product ensemble as ‘bias product ensemble’ in Sec
3.2:
$\mathbf{\hat{p_{i}}}=softmax(log(\mathbf{p_{i}})+\sum_{j=1}^{m}log(\mathbf{b_{i}^{j}}))$.

BestEnsemble: We also try to ensemble the best single debiasing models. In our
experiments (Table 4), we ensemble the three reweighting models (‘ReW’ models
in column 2,4 and 6) for each bias to form the BestEnsemble model.

###  3.4 Discussions for MoE Methods

For mixture of experts model, we summarize our findings from Table 4 below:

1) For all three known biases in Sec 3.2, we find that the debiasing methods
targeting at specific known biases increase the model performance on the
corresponding adversarial datasets, e.g. for the word overlap heuristics,
BiasProd model gets 71.0% accuracy on IS-SD (HANS) test set, 7.2% higher than
baseline.

2) The bias-specific methods might not make the NLI models more robust and
generalized. For example, the methods designed for word overlap heuristics get
lower scores on PI-CD, PI-SP, IC-CS, LI-TS test sets than the baseline model.

3) The proposed debiasing merging methods BestEn (Sec 3.3) inherits the
advantages of the 4 bias-specific methods on PI-CD, IS-SD, LI-TS and ST-LM
compared with other MoE debiasing models.

##  4 Data Augmentation

In this section, we explore 3 automatic augmentation ways without collecting
new data. For fair comparison, in all the following settings, we double the
training size by automatically generating the same number of augmented
instances as the original training sets as shown in Table 5.

###  4.1 Methods

Text Swap: It is an easy-to-implement method which swaps the premise $p$ and
hypothesis $h$ sentences in the original datasets. It might be an potential
solution to combat the partial-input heuristics (Sec 2.1.1) as the superficial
patterns are not observed in the premise sentences. According to the first-
order logic rules (LI-TS in Sec 2.1.3), we can only determine the gold labels
for the swapped sentence pairs whose original labels are _contradiction_. For
the _entailment_ and _neutral_ instances, we using the ensembled RoBERTa large
model trained on ‘all4’ training set (Table 6) to label the swapped sentence
pairs.

Word Substitution: We also tried to create new training instances by flipping
the words in the hypothesis sentences. We try two ways to perform
substitution: 1) synonym: We use NLTK Bird and Loper (2004) to firstly find
the synonym candidates of the content words (including nouns, verbs and
adjectives) in the hypothesis sentences, and then we replace the content words
with their synonyms if the cosine similarity ([-1,1]) between the original
window and the window after replacement is larger than 0. The window contains
at most 3 words including the replaced word and its neighbours. We represent
that window by max-pooling over the 300d Glove Pennington et al. (2014)
embedding of the words in that window. 2) Masked LM: we randomly select 30%
content words and then load the pretrained BERT large model to perform masked
LM task. We uniformly sample from top-100 ranking candidate words (excluding
the original word) and then replace the original content word with the sampled
one.

| Adversarial Test | Generalization Power Test | Original Test Sets
---|---|---|---
| PI-CD | PI-SP | IS-SD | IS-CS | LI-LI | LI-TS | ST | Avg. | RTE | DIAG | SICK | SciTail | Avg. | DNLI | ANLI | SNLI | MNLI
RoBERTa (base) Model
D(only) | 38.5 | 48.2 | 55.6 | 40.9 | 12.6 | 72.9 | 40.9 | 44.2 | 54.9 | 9.1 | 40.9 | 39.4 | 36.1 | 92.9 | 32.6 | 42.1 | 47.0
A(only) | 64.6 | 60.6 | 57.9 | 66.9 | 92.6 | 80.8 | 68.1 | 70.2 | 80.6 | 33.8 | 51.2 | 63.7 | 57.3 | 58.9 | 49.1 | 73.6 | 78.5
S(only) | 82.2 | 64.4 | 67.4 | 62.2 | 93.2 | 80.7 | 64.6 | 73.5 | 72.5 | 36.0 | 57.8 | 49.6 | 54.0 | 58.8 | 31.3 | 91.3 | 79.9
M(only) | 76.6 | 80.9 | 72.0 | 74.1 | 89.6 | 85.3 | 66.4 | 77.8 | 80.9 | 42.1 | 55.9 | 69.0 | 62.0 | 59.3 | 29.4 | 84.2 | 87.4
M+S | 82.8 | 80.1 | 73.3 | 74.4 | 91.8 | 85.6 | 67.8 | 79.4 | 81.2 | 40.7 | 57.5 | 67.4 | 61.7 | 60.5 | 28.3 | 91.7 | 87.4
M+S+D | 82.7 | 79.8 | 75.1 | 72.9 | 92.1 | 84.7 | 68.1 | 79.3 | 80.4 | 40.9 | 57.1 | 68.3 | 61.8 | 92.8 | 30.3 | 91.7 | 87.7
All4 | 82.6 | 81.7 | 77.0 | 74.7 | 94.7 | 85.3 | 69.1 | 80.7 | 83.7 | 41.9 | 57.3 | 70.5 | 63.4 | 93.0 | 49.2 | 91.9 | 87.7
All4+SR | 82.6 | 82.5 | 74.7 | 73.8 | 95.2 | 86.0 | 69.0 | 80.5 | 83.9 | 41.3 | 57.3 | 69.6 | 63.0 | 92.8 | 49.1 | 91.7 | 87.8
All4+PR | 83.4 | 79.5 | 75.5 | 73.8 | 94.6 | 85.5 | 69.1 | 80.2 | 83.8 | 44.0 | 57.5 | 70.5 | 64.0 | 92.9 | 51.2 | 91.9 | 87.6
RoBERTa (large) Model
All4 | 84.6 | 83.8 | 79.6 | 79.3 | 94.9 | 88.6 | 71.6 | 83.2 | 87.6 | 50.2 | 57.9 | 73.1 | 67.2 | 93.2 | 55.5 | 92.7 | 90.4
All4+ME | 85.0 | 81.4 | 80.1 | 77.7 | 95.7 | 88.7 | 72.2 | 83.0 | 87.2 | 47.4 | 58.0 | 73.7 | 66.6 | 93.3 | 54.8 | 93.0 | 90.2
All4+SE | 85.0 | 81.9 | 77.5 | 77.9 | 95.4 | 89.2 | 72.5 | 82.8 | 88.5 | 49.3 | 57.9 | 73.9 | 67.4 | 93.3 | 55.7 | 93.0 | 90.6
Table 6: Performance of RoBERTa model trained on different datasets using
multiple reweighting and ensemble strategries (Sec 5). ‘D’, ‘A’, ‘S’, ‘M’,
‘All4’ denotes DNLI, ANLI, SNLI, MNLI and the merge of all 4 datasets
respectively. ‘M+S’ is created by merging MNLI and SNLI datasets, same
principle in other settings. ‘ME’ and ‘SE’ denote the ensemble strategies in
Sec 5.2: the ensemble of 3 distinct models: BERT(large), XLNet(large) and
RoBERTa(large) and the ensemble of 3 RoBERTa(large) models. ‘SR’ and ‘PR’
refer to the size-based and performance-based reweighting in Sec 5.1. Here for
‘PR’ we use the average score of all the listed tests in ‘D(only)’, ‘A(only)’,
‘S(only)’ and ‘M(only)’ rows as their corresponding performance.

Paraphrase: We create the paraphrases for the original hypothesis sentences by
back translation Wieting and Gimpel (2018); Hu et al. (2019) using the
pretrained English-German and German-English machine translation models Ng et
al. (2019). To increase the diversity, we use beam search (size=5) for German-
English translation and get the paraphrase by sampling from the candidate
sentences.

###  4.2 Quality Analysis

To assess the quality of augmented data, we conduct both automatic and human
evaluation. For automatic evaluation, we use the best NLI model
(RoBERTa(large) model with ‘All4+SinEN’ in Table 6) in this paper to judge if
the labels of augmented data are consistent with the predictions of our best
NLI model. For human evaluation, we firstly sample 50 instances from each
augmented training data and then hire 3 human annotators to decide the
relation for the sentences pairs. We shuffle the 200 instances without showing
the annotators the augmentation method for certain instances. We also ask the
annotators to be objective and not to guess the augmentation methods and then
use the majority vote for final annotation. The accuracy of text swap, word
substitution (synonym), word substitution (MLM) and paraphrase are 84.0%,
82.0%, 88.1% and 92.9% respectively based on human-annotated gold labels.
Correspondingly, word substitution (synonym), word substitution (MLM) and
paraphrase get 76.9 %, 83.5% and 94.5% accuracy on the automatic evaluation.
Paraphrase augmentation is shown to have the highest quality among the four
methods.

###  4.3 Discussions for Data Augmentations

For Data Augmentation, we show the performance of a BERT base model using
different data augmentation methods in Table 5.

Text swap method increases the model performance on IS-CS, LI-LI, LI-TS and ST
test sets, as it can make the data distribution in the premises and hypotheses
more balanced. It is also an easy-to-implement method which could serve as a
baseline to evaluate other automatic data augmentation methods. For the other
two methods, the fragility of NLI models to partial input and inter-sentence
heuristics is partially due to the rigid word-label concurrence (PI-SP in Sec
2.1.1) or word-to-word mapping (IS-SD, IS-CS in Sec 2.1.2). More diverse
lexical choices via word substitution or paraphrase might help to relieve the
biases caused by these heuristics. We see that ‘word sub’ in Table 5
outperforms baseline on IS-CS, LI-TS and ST; ‘paraphrase’ outperforms the
baseline on IS-SD, LI-TS. However, these two methods get lower scores on other
adversarial and general purpose datasets as these debiasing techniques bias
the model towards being robust to a specific bias, so it compensates by
trading off performance.

##  5 Dataset Merging and Model Ensemble

In this section we explore 1) to what extend larger dataset and ensemble would
make the NLI models more robust to distinct adverserial datasets. 2) what is
the best way to combine the large-scale NLI training sets in very different
domains.

###  5.1 Merging Heterogeneous Datasets

To set up more diverse and stronger baselines for the proposed benchmark
datasets, we use 4 large-scale training datasets: SNLI, MNLI, DNLI and ANLI
for the following experiments. Those training sets are created using different
strategies. Specifically, SNLI and MNLI are created in a human elicited way
Poliak et al. (2018b): the human annotators are asked to write a hypothesis
sentence according to the given premise and label. DNLI recasts other NLP
tasks to fit in the form of NLI. ANLI is created as hard datasets that may
fail the models. Since those datasets vary in sizes, domains and collection
processes, they might have different contribution to the final predictions.
Here we investigate two instance reweighting methods accordingly.

Notations: suppose we have $k$ training sets $\\{T_{i}\\}_{i=1}^{k}$ whose
sizes are $\\{n_{i}\\}_{i=1}^{k}$. The accuracies of a baseline model trained
on $\\{T_{i}\\}_{i=1}^{k}$ are $\\{p_{i}\\}_{i=1}^{k}$ respectively. $p_{i}$
can be the average scores of multiple test sets or the score on an single in-
domain/ out-of-domain/ adversarial test set.

Size-based reweighting (SR): Smaller training sets might have less influence
on the models than larger ones. In this setting, we try to increase the weight
of smaller datasets so that each dataset contributes more equally to the final
predictions. We implement this reweighting method by replacing the
$\alpha_{i}$ in Sec 3.1 with $(\sum_{k}n_{k})/n_{i}(i\in T_{i})$.

Performance-based reweighting (PR): Different training sets may vary in
annotation quality and collection process thus have distinct model
performance. In this setting, we reweight the training instances with the
performance of a baseline model on the specific training sets. We still use
the instance weights in Sec 3.1 with $\alpha_{i}=p_{i}/(\sum_{k}p_{k})(i\in
T_{i})$.

###  5.2 Model Ensemble

We try two modes for model ensemble: mixed and single mode. In the mixed mode,
we ensemble three different models (BERT, XLNet, RoBERTa) while in the single
mode, we ensemble three same models (RoBERTa*3).

###  5.3 Discussions

For Dataset merging and model ensemble, according to Table 6, We find that:

1) Incorporating heterogeneous training data is a straightforward method to
enhance the robustness of NLI models. Empirically we see incorporating
datasets with adversarial human-in-the-loop annotating (e.g. ANLI) is more
efficient that incorporating automatically constructed dataset without human
curation (e.g. DNLI).

2) In RoBERTa base model, the ‘All4+PR’ model get higher scores on diagnostic
and ANLI test sets than ‘All4’ baseline, which shows that increasing the
weight of higher quality dataset may help to increase accuracy on certain test
sets. Notably, performance based reweighting helps the model gain 2 points
(49.2 vs 51.2) on ANLI compared with baseline model while keeping the
inference ability on DNLI, SNLI and MNLI test sets.

3) In RoBERTa large model, we see that on some datasets, like IS-SD, the mixed
ensemble model may even outperform the single ensemble model even if its two
components (XLNet and BERT) are less powerful than those (RoBERTa) in single
ensemble mode.

Figure 2: Per-layer analysis for RoBERTa(base) model trained on MultiNLI. Darker blue denotes higher score. ‘max’ represents the maxpooled vector across all layers. Nearly all test sets except ANLI get higher scores by using higher layers. On ANLI, the performance of the first 4 layers are close to random guess while that of higher layers is about 4 point lower than random guess. Labels | Transformation | Datasets
---|---|---
($\neg$E, E) | C$\Rightarrow\neg$E, N$\Rightarrow\neg$E | IS-SD, RTE, DNLI
($\neg$C, C) | E$\Rightarrow\neg$C, N$\Rightarrow\neg$C | LI-TS
(E, C) | - | LI-LI
(N, E) | - | SciTail
Table 7: How we evaluate the test sets with only two labels in 3-way NLI
classification. E,C,N,$\neg$ means _entailment_ , _contradiction_ , _neutral_
and _not_ respectively. $\Rightarrow$ means changing the left-hand side model
prediction with the right-hand side label while evaluation.

##  6 Experimental Settings

###  6.1 Implementation Details

We set up both pretrained and non-pretrained model baselines for the proposed
evaluation bechmarks. We rerun their public available codebases Wolf et al.
(2019), including InferSent Conneau et al. (2017)
666https://github.com/facebookresearch/InferSent (w/ and w/o Elmo Peters et
al. (2018)), DAM Parikh et al. (2016) 777https://github.com/harvardnlp/decomp-
attn, ESIM Chen et al. (2017)888https://github.com/coetaur0/ESIM, BERT
(uncased) Devlin et al. (2019), XLNet (cased) Yang et al. (2019) and RoBERTa
Liu et al. (2019), 999https://github.com/huggingface/transformers. we map the
vector at the position of the ‘[CLS]’ token in the pretrained models to three-
way NLI classification via linear transformation. We show the per-layer
analyses for RoBERTa model in Table 2. We try to reduce the randomness of our
experiments by 3 runs using different random seeds. We report the median of
the 3 runs for all the tables except the ensemble-related (Sec 5.2)
experiments in Table 6. Table 7 shows how we evaluate the test sets with only
two labels in 3-way NLI classification.

| RTE | SICK | SciTail | DNLI | ANLI | SNLI | MNLI
---|---|---|---|---|---|---|---
Origin | 75.4 | 54.2 | 66.1 | 54.2 | 27.7 | 80.0 | 83.5
Mixed | 75.5 | 54.3 | 67.3 | 54.8 | 27.4 | 79.9 | 83.4
Oracle | 75.5 | 55.2 | 67.3 | 56.7 | 28.0 | 80.3 | 83.5
Table 8: The performance of BERT base model under different model selection
strategies.

###  6.2 Model Selection Strategy

Since we test the NLI models on multiple general-purpose dataset. it is an
important question how we choose the dev set. We explore 3 different model
selection settings:

1) Origin: using the original in-domain dev set.

2) Mixed: using the merged dev sets which include all the instances in the in-
domain and extra dev sets in generalization power tests.

3) Oracle: tuning the model for each generalization power test using its own
dev set.

We show the performance of a BERT base model trained on MultiNLI utilizing the
above mentioned model selection strategies in Table 8. In this paper we use
the ‘origin’ mode, as it is too expensive to use the ‘oracle’ strategy in all
experiments, besides we did not see much difference between the ‘mixed’ and
‘origin’ modes. Notably when we merge different training sets, we also merge
their dev sets correspondingly to form a unified in-domain dev set in Table 6.

##  7 Related Work

Bias in NLI: The bias in the data annotation exists in many tasks, e.g.
lexical inference Levy et al. (2015), visual question answering Goyal et al.
(2017), ROC story cloze Cai et al. (2017); Schwartz et al. (2017) etc. The NLI
models are shown to be sensitive to the compositional features in premises and
hypothesesNie et al. (2019); Dasgupta et al. (2018), data permutations
Schluter and Varab (2018); Wang et al. (2019c) and vulnerable to adversarial
examples Iyyer et al. (2018); Minervini and Riedel (2018); Glockner et al.
(2018) and crafted stress test Geiger et al. (2018); Naik et al. (2018). Other
evidences of artifacts include sentence occurrence Zhang et al. (2019),
syntactic heuristics between hypotheses and premises McCoy et al. (2019) and
black-box clues derived from neural models Gururangan et al. (2018); Poliak et
al. (2018b); He et al. (2019). Rudinger et al. (2017) showed hypotheses in
SNLI has the evidence of gender, racial stereotypes, etc. Sanchez et al.
(2018) analysed the behaviour of NLI models and the factors to be more robust.
Feng et al. (2019) discussed how to use partial-input baseline in future
dataset creation. Belinkov et al. (2019); Clark et al. (2019); He et al.
(2019); Yaghoobzadeh et al. (2019); Ding et al. (2020); Utama et al. (2020a)
proposed efficient methods to mitigate a particular known bias in NLI. Wu et
al. (2020); Utama et al. (2020b) investigates the method to improve model
generalization by debiasing unknown dataset bias in natural language
understanding and question answering tasks.

Benchmark collection in NLI: GLUE Wang et al. (2019b, a) benchmark contains
several NLI-related benchmark datasets. However it does not include
adversarial test sets, domain specific test Romanov and Shivade (2018);
Ravichander et al. (2019). Researchers create NLI datasets using different
collection criteria, such as recasting other NLP tasks to NLI Poliak et al.
(2018a), iteratively filtering adversarial training data by model decisions
Bras et al. (2020) (model-in-the-loop), counterfactually augmenting training
data by human editing examples to break the model Kaushik et al. (2020)
(human-in-the-loop) and multi-round annotating depending on both human and
model decisions Nie et al. (2020).

##  8 Conclusions

We try to investigate how to build robust and generalized NLI models by model-
agnostic debiasing strategies, including mixture of experts ensemble (MoE),
data augmentation (DA), dataset merging and model ensemble, and benchmark
these methods on various adversarial and general purpose datasets. Our
findings suggest model-level MoE ensemble, text swap DA and performance based
dataset merging would effectively combat multiple (though not all) distinct
biases.

Although we haven’t found a debiasing strategy that can guarantee the NLI
models to be more robust on every adversarial dataset used in this paper, we
leave the question of whether such a debiasing method exists for future
research.

## Acknowledgments

We would like to thank Sam Wiseman and Kevin Gimpel for very thoughtful
discussions, and the anonymous reviewers for their helpful feedback. This
project is supported by NSFC (No. 61876004, No. U19A2065) and Beijing Academy
of Artificial Intelligence (BAAI).

## References

  * Belinkov et al. (2019) Yonatan Belinkov, Adam Poliak, Stuart Shieber, Benjamin Van Durme, and Alexander Rush. 2019.  On adversarial removal of hypothesis-only bias in natural language inference.  pages 256–262, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Bird and Loper (2004) Steven Bird and Edward Loper. 2004.  NLTK: The natural language toolkit.  In _Proceedings of the ACL Interactive Poster and Demonstration Sessions_ , pages 214–217, Barcelona, Spain. Association for Computational Linguistics.
  * Bowman et al. (2015) Samuel R. Bowman, Gabor Angeli, Christopher Potts, and Christopher D. Manning. 2015\.  A large annotated corpus for learning natural language inference.  In _Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing_ , pages 632–642, Lisbon, Portugal. Association for Computational Linguistics.
  * Bras et al. (2020) Ronan Le Bras, Swabha Swayamdipta, Chandra Bhagavatula, Rowan Zellers, Matthew E. Peters, Ashish Sabharwal, and Yejin Choi. 2020.  Adversarial filters of dataset biases.  _CoRR_ , abs/2002.04108.
  * Cai et al. (2017) Zheng Cai, Lifu Tu, and Kevin Gimpel. 2017.  Pay attention to the ending:strong neural baselines for the ROC story cloze task.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers)_ , pages 616–622, Vancouver, Canada. Association for Computational Linguistics.
  * Chen et al. (2017) Qian Chen, Xiaodan Zhu, Zhen-Hua Ling, Si Wei, Hui Jiang, and Diana Inkpen. 2017\.  Enhanced LSTM for natural language inference.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1657–1668, Vancouver, Canada. Association for Computational Linguistics.
  * Clark et al. (2019) Christopher Clark, Mark Yatskar, and Luke Zettlemoyer. 2019.  Don’t take the easy way out: Ensemble based methods for avoiding known dataset biases.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 4069–4082, Hong Kong, China. Association for Computational Linguistics.
  * Conneau et al. (2017) Alexis Conneau, Douwe Kiela, Holger Schwenk, Loïc Barrault, and Antoine Bordes. 2017.  Supervised learning of universal sentence representations from natural language inference data.  In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_ , pages 670–680, Copenhagen, Denmark. Association for Computational Linguistics.
  * Dagan et al. (2005) Ido Dagan, Oren Glickman, and Bernardo Magnini. 2005.  The PASCAL recognising textual entailment challenge.  In _Machine Learning Challenges, Evaluating Predictive Uncertainty, Visual Object Classification and Recognizing Textual Entailment, First PASCAL Machine Learning Challenges Workshop, MLCW 2005, Southampton, UK, April 11-13, 2005, Revised Selected Papers_ , volume 3944 of _Lecture Notes in Computer Science_ , pages 177–190. Springer.
  * Dagan et al. (2013) Ido Dagan, Dan Roth, Mark Sammons, and Fabio Massimo Zanzotto. 2013.  _Recognizing Textual Entailment: Models and Applications_.  Synthesis Lectures on Human Language Technologies. Morgan & Claypool Publishers.
  * Dasgupta et al. (2018) Ishita Dasgupta, Demi Guo, Andreas Stuhlmüller, Samuel Gershman, and Noah D. Goodman. 2018.  Evaluating compositionality in sentence embeddings.  In _Proceedings of the 40th Annual Meeting of the Cognitive Science Society, CogSci 2018, Madison, WI, USA, July 25-28, 2018_. cognitivesciencesociety.org.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Ding et al. (2020) Xiaoan Ding, Tianyu Liu, Baobao Chang, Zhifang Sui, and Kevin Gimpel. 2020.  Discriminatively-tuned generative classifiers for robust natural language inference.  _CoRR_ , abs/2010.03760.
  * Feng et al. (2019) Shi Feng, Eric Wallace, and Jordan Boyd-Graber. 2019.  Misleading failures of partial-input baselines.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 5533–5538, Florence, Italy. Association for Computational Linguistics.
  * Geiger et al. (2018) Atticus Geiger, Ignacio Cases, Lauri Karttunen, and Christopher Potts. 2018.  Stress-testing neural models of natural language inference with multiply-quantified sentences.  _CoRR_ , abs/1810.13033.
  * Glockner et al. (2018) Max Glockner, Vered Shwartz, and Yoav Goldberg. 2018.  Breaking NLI systems with sentences that require simple lexical inferences.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers)_ , pages 650–655, Melbourne, Australia. Association for Computational Linguistics.
  * Goyal et al. (2017) Yash Goyal, Tejas Khot, Douglas Summers-Stay, Dhruv Batra, and Devi Parikh. 2017\.  Making the V in VQA matter: Elevating the role of image understanding in visual question answering.  In _2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017_ , pages 6325–6334. IEEE Computer Society.
  * Gururangan et al. (2018) Suchin Gururangan, Swabha Swayamdipta, Omer Levy, Roy Schwartz, Samuel Bowman, and Noah A. Smith. 2018.  Annotation artifacts in natural language inference data.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 2 (Short Papers)_ , pages 107–112, New Orleans, Louisiana. Association for Computational Linguistics.
  * He et al. (2019) He He, Sheng Zha, and Haohan Wang. 2019.  Unlearn dataset bias in natural language inference by fitting the residual.  In _Proceedings of the 2nd Workshop on Deep Learning Approaches for Low-Resource NLP (DeepLo 2019)_ , pages 132–142, Hong Kong, China. Association for Computational Linguistics.
  * Hu et al. (2019) J. Edward Hu, Huda Khayrallah, Ryan Culkin, Patrick Xia, Tongfei Chen, Matt Post, and Benjamin Van Durme. 2019.  Improved lexically constrained decoding for translation and monolingual rewriting.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 839–850, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Iyyer et al. (2018) Mohit Iyyer, John Wieting, Kevin Gimpel, and Luke Zettlemoyer. 2018.  Adversarial example generation with syntactically controlled paraphrase networks.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 1875–1885, New Orleans, Louisiana. Association for Computational Linguistics.
  * Kaushik et al. (2020) Divyansh Kaushik, Eduard H. Hovy, and Zachary Chase Lipton. 2020.  Learning the difference that makes A difference with counterfactually-augmented data.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net.
  * Khot et al. (2018) Tushar Khot, Ashish Sabharwal, and Peter Clark. 2018.  Scitail: A textual entailment dataset from science question answering.  In _Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, (AAAI-18), the 30th innovative Applications of Artificial Intelligence (IAAI-18), and the 8th AAAI Symposium on Educational Advances in Artificial Intelligence (EAAI-18), New Orleans, Louisiana, USA, February 2-7, 2018_ , pages 5189–5197. AAAI Press.
  * Levy et al. (2015) Omer Levy, Steffen Remus, Chris Biemann, and Ido Dagan. 2015.  Do supervised distributional methods really learn lexical inference relations? In _Proceedings of the 2015 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 970–976, Denver, Colorado. Association for Computational Linguistics.
  * Liu et al. (2020) Tianyu Liu, Zheng Xin, Baobao Chang, and Zhifang Sui. 2020.  HypoNLI: Exploring the artificial patterns of hypothesis-only bias in natural language inference.  In _Proceedings of the 12th Language Resources and Evaluation Conference_ , pages 6852–6860, Marseille, France. European Language Resources Association.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  RoBERTa: A robustly optimized BERT pretraining approach.  _CoRR_ , abs/1907.11692.
  * Marelli et al. (2014) Marco Marelli, Luisa Bentivogli, Marco Baroni, Raffaella Bernardi, Stefano Menini, and Roberto Zamparelli. 2014.  SemEval-2014 task 1: Evaluation of compositional distributional semantic models on full sentences through semantic relatedness and textual entailment.  In _Proceedings of the 8th International Workshop on Semantic Evaluation (SemEval 2014)_ , pages 1–8, Dublin, Ireland. Association for Computational Linguistics.
  * Matthews (1975) Brian W Matthews. 1975.  Comparison of the predicted and observed secondary structure of t4 phage lysozyme.  _Biochimica et Biophysica Acta (BBA)-Protein Structure_ , 405(2):442–451.
  * McCoy et al. (2019) Tom McCoy, Ellie Pavlick, and Tal Linzen. 2019.  Right for the wrong reasons: Diagnosing syntactic heuristics in natural language inference.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 3428–3448, Florence, Italy. Association for Computational Linguistics.
  * Minervini and Riedel (2018) Pasquale Minervini and Sebastian Riedel. 2018.  Adversarially regularising neural NLI models to integrate logical background knowledge.  In _Proceedings of the 22nd Conference on Computational Natural Language Learning_ , pages 65–74, Brussels, Belgium. Association for Computational Linguistics.
  * Naik et al. (2018) Aakanksha Naik, Abhilasha Ravichander, Norman Sadeh, Carolyn Rose, and Graham Neubig. 2018.  Stress test evaluation for natural language inference.  In _Proceedings of the 27th International Conference on Computational Linguistics_ , pages 2340–2353, Santa Fe, New Mexico, USA. Association for Computational Linguistics.
  * Ng et al. (2019) Nathan Ng, Kyra Yee, Alexei Baevski, Myle Ott, Michael Auli, and Sergey Edunov. 2019\.  Facebook FAIR’s WMT19 news translation task submission.  In _Proceedings of the Fourth Conference on Machine Translation (Volume 2: Shared Task Papers, Day 1)_ , pages 314–319, Florence, Italy. Association for Computational Linguistics.
  * Nie et al. (2019) Yixin Nie, Yicheng Wang, and Mohit Bansal. 2019.  Analyzing compositionality-sensitivity of NLI models.  In _The Thirty-Third AAAI Conference on Artificial Intelligence, AAAI 2019, The Thirty-First Innovative Applications of Artificial Intelligence Conference, IAAI 2019, The Ninth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2019, Honolulu, Hawaii, USA, January 27 - February 1, 2019_ , pages 6867–6874. AAAI Press.
  * Nie et al. (2020) Yixin Nie, Adina Williams, Emily Dinan, Mohit Bansal, Jason Weston, and Douwe Kiela. 2020.  Adversarial NLI: A new benchmark for natural language understanding.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 4885–4901, Online. Association for Computational Linguistics.
  * Parikh et al. (2016) Ankur Parikh, Oscar Täckström, Dipanjan Das, and Jakob Uszkoreit. 2016.  A decomposable attention model for natural language inference.  In _Proceedings of the 2016 Conference on Empirical Methods in Natural Language Processing_ , pages 2249–2255, Austin, Texas. Association for Computational Linguistics.
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher Manning. 2014.  GloVe: Global vectors for word representation.  In _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1532–1543, Doha, Qatar. Association for Computational Linguistics.
  * Peters et al. (2018) Matthew Peters, Mark Neumann, Mohit Iyyer, Matt Gardner, Christopher Clark, Kenton Lee, and Luke Zettlemoyer. 2018.  Deep contextualized word representations.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 2227–2237, New Orleans, Louisiana. Association for Computational Linguistics.
  * Poliak et al. (2018a) Adam Poliak, Aparajita Haldar, Rachel Rudinger, J. Edward Hu, Ellie Pavlick, Aaron Steven White, and Benjamin Van Durme. 2018a.  Collecting diverse natural language inference problems for sentence representation evaluation.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 67–81, Brussels, Belgium. Association for Computational Linguistics.
  * Poliak et al. (2018b) Adam Poliak, Jason Naradowsky, Aparajita Haldar, Rachel Rudinger, and Benjamin Van Durme. 2018b.  Hypothesis only baselines in natural language inference.  In _Proceedings of the Seventh Joint Conference on Lexical and Computational Semantics_ , pages 180–191, New Orleans, Louisiana. Association for Computational Linguistics.
  * Ravichander et al. (2019) Abhilasha Ravichander, Aakanksha Naik, Carolyn Rose, and Eduard Hovy. 2019.  EQUATE: A benchmark evaluation framework for quantitative reasoning in natural language inference.  In _Proceedings of the 23rd Conference on Computational Natural Language Learning (CoNLL)_ , pages 349–361, Hong Kong, China. Association for Computational Linguistics.
  * Romanov and Shivade (2018) Alexey Romanov and Chaitanya Shivade. 2018.  Lessons from natural language inference in the clinical domain.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 1586–1596, Brussels, Belgium. Association for Computational Linguistics.
  * Rudinger et al. (2017) Rachel Rudinger, Chandler May, and Benjamin Van Durme. 2017.  Social bias in elicited natural language inferences.  In _Proceedings of the First ACL Workshop on Ethics in Natural Language Processing_ , pages 74–79, Valencia, Spain. Association for Computational Linguistics.
  * Sanchez et al. (2018) Ivan Sanchez, Jeff Mitchell, and Sebastian Riedel. 2018.  Behavior analysis of NLI models: Uncovering the influence of three factors on robustness.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 1975–1985, New Orleans, Louisiana. Association for Computational Linguistics.
  * Schluter and Varab (2018) Natalie Schluter and Daniel Varab. 2018.  When data permutations are pathological: the case of neural natural language inference.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing_ , pages 4935–4939, Brussels, Belgium. Association for Computational Linguistics.
  * Schwartz et al. (2017) Roy Schwartz, Maarten Sap, Ioannis Konstas, Leila Zilles, Yejin Choi, and Noah A. Smith. 2017.  The effect of different writing tasks on linguistic style: A case study of the ROC story cloze task.  In _Proceedings of the 21st Conference on Computational Natural Language Learning (CoNLL 2017)_ , pages 15–25, Vancouver, Canada. Association for Computational Linguistics.
  * Tsuchiya (2018) Masatoshi Tsuchiya. 2018.  Performance impact caused by hidden bias of training data for recognizing textual entailment.  In _Proceedings of the Eleventh International Conference on Language Resources and Evaluation (LREC 2018)_ , Miyazaki, Japan. European Language Resources Association (ELRA).
  * Utama et al. (2020a) Prasetya Ajie Utama, Nafise Sadat Moosavi, and Iryna Gurevych. 2020a.  Mind the trade-off: Debiasing NLU models without degrading the in-distribution performance.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 8717–8729, Online. Association for Computational Linguistics.
  * Utama et al. (2020b) Prasetya Ajie Utama, Nafise Sadat Moosavi, and Iryna Gurevych. 2020b.  Towards debiasing NLU models from unknown biases.  _CoRR_ , abs/2009.12303.
  * Wang et al. (2019a) Alex Wang, Yada Pruksachatkun, Nikita Nangia, Amanpreet Singh, Julian Michael, Felix Hill, Omer Levy, and Samuel R. Bowman. 2019a.  Superglue: A stickier benchmark for general-purpose language understanding systems.  In _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, 8-14 December 2019, Vancouver, BC, Canada_ , pages 3261–3275.
  * Wang et al. (2019b) Alex Wang, Amanpreet Singh, Julian Michael, Felix Hill, Omer Levy, and Samuel R. Bowman. 2019b.  GLUE: A multi-task benchmark and analysis platform for natural language understanding.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_. OpenReview.net.
  * Wang et al. (2019c) Haohan Wang, Da Sun, and Eric P. Xing. 2019c.  What if we simply swap the two text fragments? A straightforward yet effective way to test the robustness of methods to confounding signals in nature language inference tasks.  In _The Thirty-Third AAAI Conference on Artificial Intelligence, AAAI 2019, The Thirty-First Innovative Applications of Artificial Intelligence Conference, IAAI 2019, The Ninth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2019, Honolulu, Hawaii, USA, January 27 - February 1, 2019_ , pages 7136–7143. AAAI Press.
  * Wieting and Gimpel (2018) John Wieting and Kevin Gimpel. 2018.  ParaNMT-50M: Pushing the limits of paraphrastic sentence embeddings with millions of machine translations.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 451–462, Melbourne, Australia. Association for Computational Linguistics.
  * Williams et al. (2018) Adina Williams, Nikita Nangia, and Samuel Bowman. 2018.  A broad-coverage challenge corpus for sentence understanding through inference.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 1112–1122, New Orleans, Louisiana. Association for Computational Linguistics.
  * Wolf et al. (2019) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Rémi Louf, Morgan Funtowicz, and Jamie Brew. 2019.  Huggingface’s transformers: State-of-the-art natural language processing.  _CoRR_ , abs/1910.03771.
  * Wu et al. (2020) Mingzhu Wu, Nafise Sadat Moosavi, Andreas Rücklé, and Iryna Gurevych. 2020.  Improving QA generalization by concurrent modeling of multiple biases.  _CoRR_ , abs/2010.03338.
  * Yaghoobzadeh et al. (2019) Yadollah Yaghoobzadeh, Remi Tachet des Combes, Timothy J. Hazen, and Alessandro Sordoni. 2019.  Robust natural language inference models with example forgetting.  _CoRR_ , abs/1911.03861.
  * Yang et al. (2019) Zhilin Yang, Zihang Dai, Yiming Yang, Jaime G. Carbonell, Ruslan Salakhutdinov, and Quoc V. Le. 2019.  XLNet: Generalized autoregressive pretraining for language understanding.  In _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, 8-14 December 2019, Vancouver, BC, Canada_ , pages 5754–5764.
  * Zhang et al. (2019) Guanhua Zhang, Bing Bai, Jian Liang, Kun Bai, Shiyu Chang, Mo Yu, Conghui Zhu, and Tiejun Zhao. 2019.  Selection bias explorations and debias methods for natural language sentence matching datasets.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 4418–4429, Florence, Italy. Association for Computational Linguistics.


∎

e1e-mail<EMAIL_ADDRESS>
11institutetext: Institute of Theoretical Physics, Jagiellonian University,
ul. Łojasiewicza 11, 30-348 Kraków, Poland

# Finite volume effects in the McLerran-Venugopalan initial condition for the
JIMWLK equation

Piotr Korcyle1,addr1

(Received: date / Accepted: date)

###### Abstract

We revisit the numerical construction of the initial condition for the dipole
amplitude from the McLerran-Venugopalan model in the context of the JIMWLK
evolution equation. We observe large finite volume effects induced by the
Poisson equation formulated on a torus. We show that the situation can be
partially cured by introducing an infrared regularization. We propose a
procedure which has negligible finite volume corrections. The control of the
finite volume and finite lattice spacings effects is crucial when considering
the numerical solutions of the JIMWLK evolution equation with the collinear
improvement.

###### Keywords:

Heavy Ion Phenomenology QCD Phenomenology gluon saturation Langevin equation
Numerical simulations

###### pacs:

12.38.Lg 12.38.Mh 12.39.St

††journal: Eur. Phys. J. C

##  1 Introduction

Evolution equations are necessary elements of the theoretical frameworks
developed to describe data from Deep Inelastic Scattering (DIS) experiments
gathered from various facilities worldwide. Experimental data was collected
for a wide range of kinematical variables describing this process, in
particular the momentum transfer $Q^{2}$ and longitudinal momentum fraction
$x$. Our understanding of Quantum Chromodynamics, the theory behind these
experiments, provides us with predictions on the dependence of the DIS cross-
section on $Q^{2}$ and $x$ and is summarized by the mentioned evolution
equations (see Ref. Kovchegov:2012mbw  for a review). Due to the complexity of
QCD the known evolution equations are based on perturbative kernels and are
valid only in a limited region of phase space. In many situations these
evolution equations can be cast into the form of parabolic differential
equation and introduce a Hamiltonian governing the dynamics in the evolution
variable. In that case, solving the problem requires the knowledge of the
initial condition at some initial value of $x_{0}$ or $Q_{0}$. That object is
inherently non-perturbative and therefore usually it is not available and has
to be taken from some approximate model.

In the case of the JIMWLK equation Balitsky:1995ub ; JalilianMarian:1997jx ;
JalilianMarian:1997gr ; JalilianMarian:1997dw ; Kovner:2000pt ; Kovner:1999bj
; Weigert:2000gi ; Iancu:2000hn ; Ferreiro:2001qy , which motivated our
present work, the initial condition is provided by the parton probability
distribution in the transverse momentum. The description of the DIS data
requires only the simplest, two-point gluon probability distribution also
known as the dipole distribution (which also enters other processes, see for
example Ref. 2019 ). Typically, one resorts to the Color Glass Condensate
effective theory of QCD (see e.g. Gelis:2010nm ) which offers a computational
basis where the dipole distribution can be modelled for instance in the
McLerran-Venugopalan model McLerran:1993ni ; McLerran:1993ka .

The JIMWLK evolution equation being a non-linear equation does not admit
general analytically known solutions. One can employ numerical methods
Rummukainen:2003ns  to study its solutions for a particular initial condition.
The first step of such calculation is the numerical construction of the
initial dipole amplitude in a discretized setup. In what follows we show the
details of that step and expose the problems associated with finite volume
effects which were so far overlooked. We then propose a simple solution and
present numerical evidence for its effectiveness.

The construction of the initial condition in the MV model was already
discussed in the Literature (see initial1995 ; initial1995b ; initial1999 ;
initial2000 ; initial2001 ; initial2007 ; PhysRevC.79.024909 ), also from the
numerical perspective (for example in Ref. Fukushima:2007yk ;
fukushima2008randomness ), however we found that some aspects concerning
finite volume effects were missing and we aim to close this gap with this
work. In Ref. Lappi:2007ku ; Lappi:2016qhk  the numerical values of the
parameters in the MV model were set in physical units through numerical
calculations. The construction of initial condition based on phenomenological
considerations were discussed for example in Refs. albacete1 ; albacete2  and
many others.

Our findings and observations are particularly important in the context of
collinearly improved JIMWLK evolution equation cl2009 ; cl2014 ; cl2015 ;
cl2015b . It is well-known since many years that the JIMWLK evolution equation
alone is not compatible with experimental data and the remedy to that is
expected to come from resummation of single and double logarithms in $Q^{2}$.
As was discussed in Ref. Hatta:2016ujq , such resummation can be introduced
into the Langevin formulation of the JIMWKL equation with the help of time-
ordering of the consecutive gluon interaction during the evolution and some
modifications of the JIMWLK kernel. In order to be able to study the effects
of these improvements one has to be able to disentangle them from the various
effects associated with the initial condition. This would not be possible with
the initial dipole amplitude motivated by the McLerran-Venugopalan model
derived numerically with the numerical setup commonly used due to large finite
volume effects. In contrast, our proposal for the initial condition does not
suffer from any spurious finite volume and finite lattice spacing effects.
Hence, any signs of such effects can be readily associated to the physics of
the JIMWLK equation and the collinear improvement. The discussion of the
solutions of the JIMWLK equation with the collinear improvement are beyond the
scope of this work and will be presented elsewhere.

General information about the open-access numerical code used in this work can
be found in Ref. Korcyl:2020orf  and in the associated `git` repository.
Complementary study of the systematic effects of the solutions of the JIMWLK
equations was presented in Ref. Cali:2021tsh . All calculations in that
reference were conducted at a single value of the volume, hence the
conclusions presented there do not touch the issues contained in the current
work.

##  2 Dipol gluon probability distribution

In the McLerran-Venugopalan (MV) model McLerran:1993ka ; McLerran:1993ni ; MV3
the color charge density $\rho^{a}(\mathbf{x})$ is defined on the infinite
transverse plane and follows an independent Gaussian distribution at each
point,

$\langle\rho^{a}(\mathbf{x})\rho^{b}(\mathbf{y})\rangle=\delta_{a,b}\
\delta^{(2)}(\mathbf{x}-\mathbf{y})\ \mu^{2},$ (1)

where we dropped the dependence on $x^{+}$ and $x^{-}$ light-cone coordinates
because it will not play any role in the following discussion. The $\mu$
parameter is the main parameter of the MV model and has the interpretation of
the density of color charges per area element. Following the philosophy of the
CGC effective theory, the static color charges represent the large-$x$ gluons
inside the probed hadron whereas the small-$x$ gluons result from these
charges via the Yang-Mills equations. The latter reduce to the Poisson
equation,

$\nabla^{2}A^{a}(\mathbf{x})=\rho^{a}(\mathbf{x})\,,$ (2)

which can be solved for the color potential $A^{a}(\mathbf{x})$ using the
known Green function of the two-dimensional Poisson equation for any given
$\rho^{a}(\mathbf{x})$. Since the interplay of boundary conditions with the
existence of solutions will be crucial in what follows let us be more
explicit.

On the infinite transverse plane (following the original MV construction)
there exists a Green function for the two-dimensional Poisson equation,

$G(\mathbf{x},\mathbf{z})=\frac{1}{2\pi}\ln\big{|}\mathbf{x}-\mathbf{z}\big{|}.$
(3)

Therefore, for any color charges configuration we can obtain the color
potential,

$A^{a}(\mathbf{x})=\int_{-\infty}^{\infty}d\mathbf{z}\ \rho^{a}(\mathbf{z})\
G(\mathbf{x},\mathbf{z}).$ (4)

Proceeding, we can calculate the expectation value of the product of two color
potentials,

$\langle A^{a}(\mathbf{x})A^{b}(\mathbf{y})\rangle=\\\
=\int_{-\infty}^{\infty}d\mathbf{z}\ d\mathbf{z^{\prime}}\
G(\mathbf{x},\mathbf{z})\ G(\mathbf{y},\mathbf{z^{\prime}})\
\langle\rho^{a}(\mathbf{z})\rho^{b}(\mathbf{z^{\prime}})\rangle=\\\ =\mu^{2}\
\delta_{a,b}\ \Gamma(\mathbf{x},\mathbf{y}),$ (5)

where

$\Gamma(\mathbf{x},\mathbf{y})=\int_{-\infty}^{\infty}d\mathbf{z}\
G(\mathbf{x},\mathbf{z})\ G(\mathbf{y},\mathbf{z}).$ (6)

Because,

$\frac{1}{2\pi}\ln\big{|}\mathbf{x}\big{|}=\int\frac{d^{2}\mathbf{k}}{4\pi^{2}}\frac{e^{i\mathbf{k}\cdot\mathbf{x}}}{k^{2}},$
(7)

which is the free gluon propagator, one can also see Eq.(5) as a two-gluon
interaction of the hadron with the static color potential. Introducing Wilson
lines $U(\mathbf{x})$,

$U(\mathbf{x})=e^{iA^{a}(\mathbf{x})\lambda^{a}},$ (8)

one can calculate

$C_{\textrm{MV}}(|\mathbf{x}-\mathbf{y}|)=\langle\textrm{Tr}U^{\dagger}(\mathbf{x})U(\mathbf{y})\rangle=\\\
=\exp\Big{(}-\frac{(\mathbf{x}-\mathbf{y})^{2}}{2R^{2}}\ln\big{(}\sqrt{2}\frac{\kappa
R}{|x|}+e\big{)}\Big{)},$ (9)

where the dimensionful combination $\kappa R$ corresponds to the infrared
regulator and can be thought of being of the order of
$\Lambda_{\textrm{QCD}}$. We express it using the distance $R$ which is in
fact a function of the MV parameters, $R=R(g^{2}\mu)$. $R$ is the most
convenient quantity to be used to set the scale. Preferably, $R$ and $\kappa$
should be set by fitting the gluon distribution to some experimental data.

One should note that there exist also a simpler model for the two-point gluon
initial distribution, the Golec-Biernat-Wusthoff (GBW) model
GolecBiernat:1998js  where the following shape of the initial correlation
function is postulated,

$C_{\textrm{GBW}}(|\mathbf{x}-\mathbf{y}|)=\langle\textrm{Tr}U^{\dagger}(\mathbf{x})U(\mathbf{y})\rangle=\\\
=\exp\Big{(}-\frac{(\mathbf{x}-\mathbf{y})^{2}}{2R^{2}}\Big{)}.$ (10)

Figure 1: Comparison of the dipole amplitude as a function of the distance for
the GBW and MV models.

In both cases the dimensionful parameter $R$ has the interpretation of the
saturation radius. It is typical to define it for a distribution of any shape
as the distance at which the two-point gluon correlator has the value
$e^{-\frac{1}{2}}$,

$C(|\mathbf{x}-\mathbf{y}|=R_{s})=e^{-\frac{1}{2}}.$ (11)

We illustrate the differences between $C_{\textrm{GBW}}$ and $C_{\textrm{MV}}$
in Fig. 1. These two definitions, Eqs. (9) and (10), are the usual starting
point of analytical calculations of various gluon and quark distributions
needed for estimating more cross-sections of process like DIS or with more
complex kinematics. If one intends to employ numerical formulation and solve
evolution equations numerically, the initial gluon distribution has to be
constructed numerically.

##  3 Discretization of the McLerran-Venugopalan model

In order to treat the problem numerically one has to discretize the transverse
plane, i.e. introduce a finite spacing between neighbouring sites. The
numerical framework imposes also another fundamental constraint: the
calculations are performed in a finite volume. As is typically done, in order
to preserve translational symmetry and reduce boundary effects, one applies
periodic boundary conditions which transform the infinite transverse plane
into a torus. It should be stressed that the finite lattice spacing and the
finite volume are _technical_ tools which allow to tackle the problem on a
computer and hence should not have any effects on the final results, because
they do not appear in the original formulation of the models discussed in the
previous section. In particular, final results should be obtained after the
continuum extrapolation of the lattice spacing to zero and after the infinite
volume extrapolation with the radii of the torus going to infinity. In the
current problem the two limits are independent and can be taken in any order.
Let us finally remark that the results can only be correct if they are
properly defined at finite values of the lattice spacing and volume. In order
to better illustrate the last statement let us describe the way the MV model
is typically discretized.

One starts with the discretized torus with the color charges
$\rho^{a}(\mathbf{x})$ located at each site. Due to periodicity we have,

$\rho^{a}(\mathbf{x})\equiv\rho^{a}(x,y)=\rho^{a}(x+L_{x},y)=\rho^{a}(x,y+L_{y}),$
(12)

where $L_{x}$ and $L_{y}$ are the radii of the torus in physical units, i.e.
fm. In order to simplify the discussion we set them equal, $L_{x}=L_{y}\equiv
L$. The basic object in the calculation are again the Wilson lines at fixed
position $\mathbf{x}$, $U(\mathbf{x})$. This time $\mathbf{x}$ belongs to a
discrete set,

$\mathbf{x}=(x,y)=a(n_{x},n_{y}),\qquad n_{x},n_{y}\in\mathbb{Z}$ (13)

Wilson line variables $U(\mathbf{x})$ can be obtained from the gauge
potentials $A^{a}$

$U(\mathbf{x})=\exp\left(-igA^{a}(\mathbf{x})\lambda^{a}\right)\,,$ (14)

which in turn are related to the color charges through the Poisson equation,
as in Eq.(2)

$U(\mathbf{x})=\exp\left(-igA^{a}(\mathbf{x})\lambda^{a}\right)=\exp\left(-i\frac{g\rho^{a}(\mathbf{x})\lambda^{a}}{\nabla^{2}}\right)\,.$
(15)

At this point the finite volume and the imposed boundary conditions play a
critical role. The finite volume and in particular the fact that the torus is
a compact manifold imply that only a total zero charge is allowed by the Gauss
law. As a consequence, the Poisson equation does not admit solutions on a
torus for a non-zero probe charge, i.e. the Green function Eq. (3) does not
exist. This remains true for any value of the lattice spacing, also in the
continuum limit. There are two ways of circumventing this problem.

One can impose the global condition on the random color charges that the total
charge is zero. That is physically well motivated as all hadrons are color
singlets, and hence it seems reasonable to model them with a charge
distribution which has zero net charge. However, this introduces global
correlations among different positions in the transverse plane which are not
taken into account in the original MV model and are not present in
$S_{\textrm{MV}}$ given by Eq. (9). On the technical level this is usually
achieved by dropping the zero mode of the Fourier-transformed charge
distribution fukushima2008randomness ; Marquet:2016cgx . Note that this
procedure does not suppress long-range gluon interactions and hence can be
expected to exhibit large finite volume effects.

The second way is to regularize the Poisson equation in such a way that would
allow solutions for a charge distribution with a non-zero total charge on a
compact manifold. Note that our final results correspond to the infinite
volume limit where the boundary conditions become irrelevant and the
compactness in lost. Hence, the relevant physical correlations should be
recovered when the regulator is removed and the infinite volume limit is
taken. Obviously, these two limits do not commute, for a given value of the
regulator one has to take the infinite volume limit first and only afterwards
remove the regulator.

Most of the Literature follows the combination of the two approaches. One can
circumvent the problem with finite volume effects by limiting the charge
distribution to a small region of the transverse plane. This was investigated
in Ref. PhysRevD.98.034013  where the size of the proton $B$ was introduced
which exponentially suppresses the long range interactions. In that case, when
$B$ is much smaller then $L$, finite volume effects are under control.

In this work we keep the Fourier transform unaltered and we introduce the
infrared regulator into the Poisson equation. In that case Eq. (15) becomes

$U(\mathbf{x})=\exp\left(-igA^{a}(\mathbf{x})\lambda^{a}\right)=\exp\left(-i\frac{g\rho^{a}(\mathbf{x})\lambda^{a}}{\nabla^{2}-m^{2}}\right)\,.$
(16)

One can provide a physical interpretation of the infrared regulator parameter
$m$ in terms of the gluon propagator. We substitute $1/k^{2}$ by
$1/(k^{2}+m^{2})$ with the aim of introducing by hand the effects of color
confinement of gluons at the scale of hadron size. We recover the original MV
model when $L\rightarrow\infty$ and $m\rightarrow 0$. We demonstrate how this
works in practice in the next Section.

Let us highlight a point that turns out to be important in the following
discussion. Namely, that the Green function for the regularized Poisson
equation on an infinite plane is known and is proportional to the $K_{0}$
Bessel function,

$G(\mathbf{x},\mathbf{z},m^{2})=\frac{1}{2\pi}K_{0}\big{(}m^{2}|\mathbf{x}-\mathbf{z}|\big{)},$
(17)

which for small $m$ has the asymptotics,
$G(\mathbf{x},\mathbf{z},m^{2})\sim\ln\frac{1}{2}m^{2}|\mathbf{x}-\mathbf{z}|$.
Hence, for sufficiently small values of the infrared regulator we expect to
recover the original MV model on the infinite plane.

Let us close this Section by providing last details needed for the numerical
implementation of the MV model. Technically, on each site $\mathbf{x}$ we
construct a random matrix from the SU(3) algebra by

$\rho(\mathbf{x})=\rho(\mathbf{x})^{a}\lambda^{a}\,,$ (18)

where the color sources $\rho(\mathbf{x})^{a}$ are normally distributed and
generated by the Box-Muller method. For their standard deviation, we have

$\langle
g\rho(\mathbf{x})^{a}g\rho(\mathbf{y})^{b}\rangle=g^{4}\mu^{2}\delta^{ab}\delta(\mathbf{x}-\mathbf{y})\,.$
(19)

Above, $g^{2}$ and $\mu$ are input parameters of the MV model and it turns out
that the only relevant dimensionful combination is $g^{2}\mu$. All
dimensionful quantities can be expressed in units of $g^{2}\mu$, in particular
the saturation radius $R_{s}g^{2}\mu$ and the circumference of the torus
$Lg^{2}\mu$ become dimensionless numbers. $\mu$ being the value of the charge
density parametrizes the line of constant physics, so larger $Lg^{2}\mu$
corresponds to a larger circumference. Hence, the infinite volume limit is
obtained when $Lg^{2}\mu\rightarrow\infty$.

##  4 Numerical results for the MV model on the torus

Figure 2: Dipole amplitude generated from the discretized MV model on the
torus with an infrared regulator in a small volume $g^{2}\mu L=30.72$.
Calculation done using a linear extend of $L/a=512$. Figure 3: Dipole
amplitude generated from the discretized MV model on the torus with an
infrared regulator in a large volume $g^{2}\mu L=983.04$. Calculation done
using a linear extend of $L/a=16384$.

With the MV model implemented as described above one can calculate the dipole
amplitude,

$C(|\mathbf{x}-\mathbf{y}|)=\langle\textrm{Tr}U^{\dagger}(\mathbf{x})U(\mathbf{y})\rangle$
(20)

and check if the resulting shape corresponds to the predictions of Eq. (9). We
present the results in Figs. 2 and 3 for a wide range of the infrared
regulator $m$ values. The former shows the situation in a small volume whereas
the latter in a large volume, where small and large corresponds to the value
of the $g^{2}\mu L$ parameter. The same physics should be obtained with the
same value of $g^{2}\mu$ in physical units, hence if $g^{2}\mu L_{1}<g^{2}\mu
L_{2}$ corresponds to $L_{1}<L_{2}$ in physical units. The typical values used
in the Literature vary between $5\leq g^{2}\mu L\leq 200$ initial2000 ;
initial2001 . For Fig.2 we set $g^{2}\mu L=30.72$ as was done in recent
publications Ref. Marquet:2016cgx ; Lappi:2014wya ; Mantysaari:2018zdd ,
whereas in Fig. 3 we set $g^{2}\mu L=983.04$ which is 32 times larger than in
Fig. 2.

Apart of data sets we also show the results of our attempts to identify the
shapes predicted by Eqs. (9) and (10). We have chosen the lattice spacing
small enough that the difference between the continuum extrapolated data is
smaller than the statistical uncertainties which in turn are smaller than the
symbol size. Hence, if the numerical construction is valid and is performed at
the volume large enough one should recover unintegrated gluon distribution of
Eq.(9) for some value of $R$ and $\kappa$. After performing the fit, we see
that in the small volume neither of the theoretical predictions describes the
data for any of the tested values of the infrared regulator. We note that this
parameter was changed by four orders of magnitude and only the largest value
exhibits some deviations; other three data sets lie on top of each other.

Fig. 2 allows us to conclude the for the three smallest values of the
regulator we are in the logarithmic regime of the Bessel function and hence
the gluons are not constrained from traveling around the torus. We expect that
such long range effects introduce large finite volume effects, especially when
the volume is small. Only the data set with the the largest value of the $m$
parameter shows some effect of that regulator. However, for neither of the
values of the regulator do we recover the expected shape of the gluon
distribution.

Let us now discuss Fig. 3 where, for the same lattice spacing and the same
values of the $m$ parameter as in Fig. 2, we show the results from a much
larger volume. We notice that the differences between the different values of
$m$ are much larger. The two smallest values fall on top of each other,
whereas for the larger two the gluon distributions exhibit different shapes.
At small values of $m$ the data can be described by the $C_{\textrm{MV}}$
prediction and _cannot_ be described by the $C_{\textrm{GBW}}$ formula. We
take this as the evidence, that from the point of view of correlations between
the neighbouring Wilson lines the volume is large enough so that it can be
approximated by the infinite plane, as in the original MV model. For the
largest value of $m$ shown in the plot, the resulting gluon distribution is
very well described by the Gaussian shape postulated by the $C_{\textrm{GBW}}$
formula. We also take this as a hint that the volume is large enough so that
finite volume effects are not visible and also that the large value of $m$ has
suppressed the correlations described by the MV model.

We finish this Section by describing the results for the MV model at various
intermediate volumes and with values of the infrared regulator and
correlations small enough to fall into the logarithmic regime of the Green
function Eq. (17). Using that data we can infer the speed of the convergence
of the gluon distribution to the infinite volume limit. In Fig. 4 we show the
gluon distribution for volumes differing by a factor 2, ranging from $g^{2}\mu
L=30.72$ up to $g^{2}\mu L=1966.08$. We see slow convergence for increasing
volumes. In order to quantify that convergence we use condition given by Eq.
(11) to define the saturation radius for each volume separately and plot it
against the volume. We do it in Fig. 5, where on the horizontal axis we show
the square root of the linear extend of the torus. We also add a tentative
linear fit, which seems to be consistent with the data at hand. Thus, we can
conclude that the finite volume effects vanish very slowly, as
$\sqrt{L}=V^{\frac{1}{4}}$, and therefore one needs very large lattices in
order to correctly reproduce the MV model from the discretized setting. In the
next Section we describe a different approach which is needed to numerically
simulate the GBW model which turns out to be free of finite volume effects and
in the last Section we will demonstrate that the same method can be applied to
the MV model, thus allowing calculations using moderate computer resources.

Figure 4: Volume dependence of the dipole amplitude in the MV model on the
torus. Increasing torus size pushes the distribution to the left. Consecutive
volumes differ by a factor 4 (factor 2 in linear extend). Convergence to some
limiting distribution can be seen for very large volumes. Figure 5: Dependence
of the saturation radius extracted from the data from Fig. 4 as a function of
the fourth root of the volume. Approximate linear scaling can be seen pointing
to a very mild suppression of finite volume effects.

##  5 Discretization of the Golec-Biernat-Wusthoff model

In the GBW model the shape of the correlation function Eq. (20) is postulated
to be Gaussian, i.e. given by Eq. (10) and parametrized solely by $R$. A
numerical recipe to generate any correlation in position space between two
Wilson lines was proposed in Ref. Rummukainen:2003ns . Here, we describe a
slightly modified prescription. Note that there is no physical interpretation
of the following steps and quantities introduced below. We are interested in a
_technical_ algorithm to generate elements of the $SU(3)$ group in each site
of the torus such that their correlation function in position space takes a
desired form.

We start with the desired correlation function $C(\mathbf{x})$ defined on the
torus. Hence, $\mathbf{x}$ is understood as the shortest distance between two
sites on the torus. In order to proceed we calculate numerically the Fourier
transform of $C$ getting $C(\mathbf{k})$. Note that in the case of the GBW
model we need the Fourier transform of a Gaussian which is a Gaussian, and
hence not only $\textrm{Im}\ C(\mathbf{k})=0$ but also $C(\mathbf{k})>0$. Now,
we couple the square root of $C(\mathbf{k})$ to a white noise in the $SU(3)$
algebra. For each momentum $\mathbf{k}$ we generate random variables
$\sigma^{a}(\mathbf{k})$ that are normally distributed with a zero mean and a
variance scaled by $c\sqrt{L}$,

$\langle\sigma^{a}(\mathbf{k})\sigma^{b}(\mathbf{k^{\prime}})\rangle=c^{2}L\
\delta^{ab}\delta(\mathbf{k}-\mathbf{k^{\prime}})\,.$ (21)

where $c$ is a technical parameter which controls the ergodicity of the random
walk on the $SU(3)$ algebra. After some numerical experiments we set it to a
constant value $c=4$. The rescaling of the color sources by the desired
correlation function in momentum space

$\sigma^{a}(\mathbf{k})\rightarrow\chi^{a}(\mathbf{k})=\sqrt{C(\mathbf{k})}\
\sigma^{a}(\mathbf{k})$ (22)

guarantees that we obtain in position space the correlation given by
$C(\mathbf{x})$. The last missing step is the construction of Wilson lines
which are in the $SU(3)$ group. Hence, after a Fourier transform back to
position space we exponentiate the so obtained elements of the algebra,

$U(\mathbf{x})=\exp\left(-i\chi(\mathbf{x})\right)=\exp\left(-i\chi^{a}(\mathbf{x})\lambda^{a}\right)\,,$
(23)

In order to increase the ergodicity we found that the final Wilson line has to
be obtained from a product of a large number of elementary Wilson lines given
by Eq. (23). These should not be interpreted as multiple gluon rescatterings,
as the above procedure is purely algorithmical and the particular steps do not
correspond to any physical process. In our calculations we use 500 elementary
Wilson lines.

We provide results of this procedure in Fig. 6 where we plot the unintegrated
gluon distribution estimated with various lattice spacings and different
volumes. We see clear evidence that at the statistical precision at which we
work the procedure is free of finite lattice spacing and finite volume
effects.

Figure 6: Volume and lattice spacing dependence of the GBW dipole amplitude.
Results from different volumes and lattice spacings fall on top of each other
providing evidence of negligible finite volume and finite lattice spacing
effects.

##  6 Discretized McLerran-Venugopalan model revisited

Figure 7: Volume and lattice spacing dependence of the MV dipole amplitude.
Results from different volumes and lattice spacings fall on top of each other
providing evidence of negligible finite volume and finite lattice spacing
effects. $\kappa$ was set to 1.0 in this example.

It was found that from several phenomenological perspectives the MV model
provides a better model for the unintegrated gluon distribution than the GBW
model. In particular the behaviour at very small distances, reachable in
perturbation theory, has large impact on quantities which can be compared with
experimental data. In this region the differences between the MV and GBW
models are significant. It would be therefore desirable to dispose of the
dipole amplitude given by Eq. (9) with a similar quality as shown in Fig. 6
for the GBW model. We propose to follow the algorithm presented in the
previous Section also for the correlation function Eq. (9). The Fourier
transform is real and positive Giraud:2005vx ; Giraud:2016lgg  and therefore
it is possible to directly apply Eq. (21). All steps of the recipe can be
safely performed leading to satisfying results shown in Fig. 7.

##  7 Conclusions

In summary, in this work we have revisited the numerical construction of
unintegrated gluon distribution from the McLerran-Venugopalan model. We
reminded the results from the model defined on an infinite plane. We then
described the differences induced by a compact manifold, for instance two-
dimensional torus, which is used typically in numerical calculations. We
observed large finite volume effects present in the saturation radius
extracted from dipole amplitudes. We identified the origin of these finite
volume effects to be the Poisson equation, which is ill-defined for a non-zero
net charge on a compact manifold. We then proposed different algorithm to
generated the desired dipole amplitude which does not suffer of finite volume
and finite lattice spacing corrections. We provided numerical evidence for the
correctness of our proposal.

The contribution presented in this work provides a stable numerical setup for
studying the effects of the evolution equations. In particular, in view of the
extension of the JIMWLK equation by the collinear improvement proposed in Ref.
Hatta:2016ujq  it is crucial to have full control over the initial condition.
The presented setup allows to study the effects of the collinear improvement
alone, without any remnants of the finite volume induced by the initial
condition. We will present such results in a separate publication.

## Acknowledgements

The Author thanks K. Cichy, K. Kutak and L. Motyka for many useful
discussions. Computer time allocations ’nspt’, ’pionda’, ’tmdlangevin’,
’plgtmdlangevin2’ and ’plgquasi’ on the Prometheus supercomputer and trial
access to the Ares supercomputer, both hosted by AGH Cyfronet in Kraków,
Poland were used through the polish PLGRID consortium. The Author acknowledges
partial support by the Polish National Science Center, grant No. 2017/27/B/
ST2/02755.

## References

  * (1) Y.V. Kovchegov, E. Levin, _Quantum chromodynamics at high energy_ , vol. 33 (Cambridge University Press, 2012).  URL http://www.cambridge.org/de/knowledge/isbn/item6803159
  * (2) I. Balitsky, Nucl. Phys. B463, 99 (1996).  DOI 10.1016/0550-3213(95)00638-9
  * (3) J. Jalilian-Marian, A. Kovner, A. Leonidov, H. Weigert, Nucl. Phys. B504, 415 (1997).  DOI 10.1016/S0550-3213(97)00440-9
  * (4) J. Jalilian-Marian, A. Kovner, A. Leonidov, H. Weigert, Phys. Rev. D59, 014014 (1998).  DOI 10.1103/PhysRevD.59.014014
  * (5) J. Jalilian-Marian, A. Kovner, H. Weigert, Phys. Rev. D59, 014015 (1998).  DOI 10.1103/PhysRevD.59.014015
  * (6) A. Kovner, J.G. Milhano, H. Weigert, Phys. Rev. D62, 114005 (2000).  DOI 10.1103/PhysRevD.62.114005
  * (7) A. Kovner, J.G. Milhano, Phys. Rev. D61, 014012 (2000).  DOI 10.1103/PhysRevD.61.014012
  * (8) H. Weigert, Nucl. Phys. A703, 823 (2002).  DOI 10.1016/S0375-9474(01)01668-2
  * (9) E. Iancu, A. Leonidov, L.D. McLerran, Nucl. Phys. A692, 583 (2001).  DOI 10.1016/S0375-9474(01)00642-X
  * (10) E. Ferreiro, E. Iancu, A. Leonidov, L. McLerran, Nucl. Phys. A703, 489 (2002).  DOI 10.1016/S0375-9474(01)01329-X
  * (11) A. van Hameren, P. Kotko, K. Kutak, S. Sapeta, Physics Letters B 795, 511–515 (2019).  DOI 10.1016/j.physletb.2019.06.055.  URL http://dx.doi.org/10.1016/j.physletb.2019.06.055
  * (12) F. Gelis, E. Iancu, J. Jalilian-Marian, R. Venugopalan, Ann. Rev. Nucl. Part. Sci. 60, 463 (2010).  DOI 10.1146/annurev.nucl.010909.083629
  * (13) L.D. McLerran, R. Venugopalan, Phys. Rev. D49, 2233 (1994).  DOI 10.1103/PhysRevD.49.2233
  * (14) L.D. McLerran, R. Venugopalan, Phys. Rev. D49, 3352 (1994).  DOI 10.1103/PhysRevD.49.3352
  * (15) K. Rummukainen, H. Weigert, Nucl. Phys. A739, 183 (2004).  DOI 10.1016/j.nuclphysa.2004.03.219
  * (16) A. Kovner, L. McLerran, H. Weigert, Physical Review D 52(7), 3809–3814 (1995).  DOI 10.1103/physrevd.52.3809.  URL http://dx.doi.org/10.1103/PhysRevD.52.3809
  * (17) A. Kovner, L. McLerran, H. Weigert, Physical Review D 52(11), 6231–6237 (1995).  DOI 10.1103/physrevd.52.6231.  URL http://dx.doi.org/10.1103/PhysRevD.52.6231
  * (18) A. Krasnitz, R. Venugopalan, Nuclear Physics B 557(1-2), 237–270 (1999).  DOI 10.1016/s0550-3213(99)00366-1.  URL http://dx.doi.org/10.1016/S0550-3213(99)00366-1
  * (19) A. Krasnitz, R. Venugopalan, Physical Review Letters 84(19), 4309–4312 (2000).  DOI 10.1103/physrevlett.84.4309.  URL http://dx.doi.org/10.1103/PhysRevLett.84.4309
  * (20) A. Krasnitz, R. Venugopalan, Physical Review Letters 86(9), 1717–1720 (2001).  DOI 10.1103/physrevlett.86.1717.  URL http://dx.doi.org/10.1103/PhysRevLett.86.1717
  * (21) K. Fukushima, F. Gelis, L. McLerran, Nuclear Physics A 786(1-4), 107–130 (2007).  DOI 10.1016/j.nuclphysa.2007.01.086.  URL http://dx.doi.org/10.1016/j.nuclphysa.2007.01.086
  * (22) H. Fujii, K. Fukushima, Y. Hidaka, Phys. Rev. C 79, 024909 (2009).  DOI 10.1103/PhysRevC.79.024909.  URL https://link.aps.org/doi/10.1103/PhysRevC.79.024909
  * (23) K. Fukushima, Phys. Rev. C 76, 021902 (2007).  DOI 10.1103/PhysRevC.76.021902.  [Erratum: Phys.Rev.C 77, 029901 (2007)]
  * (24) K. Fukushima, Physical Review D 77(7), 074005 (2008)
  * (25) T. Lappi, Eur. Phys. J. C55, 285 (2008).  DOI 10.1140/epjc/s10052-008-0588-4
  * (26) T. Lappi, EPJ Web Conf. 137, 07013 (2017).  DOI 10.1051/epjconf/201713707013
  * (27) J.L. Albacete, N. Armesto, J.G. Milhano, C.A. Salgado, Physical Review D 80(3) (2009).  DOI 10.1103/physrevd.80.034031.  URL http://dx.doi.org/10.1103/PhysRevD.80.034031
  * (28) J.L. Albacete, N. Armesto, J.G. Milhano, P. Quiroga Arias, C.A. Salgado, The European Physical Journal C 71(7) (2011).  DOI 10.1140/epjc/s10052-011-1705-3.  URL http://dx.doi.org/10.1140/epjc/s10052-011-1705-3
  * (29) L. Motyka, A.M. Staśto, Physical Review D 79(8) (2009).  DOI 10.1103/physrevd.79.085016.  URL http://dx.doi.org/10.1103/PhysRevD.79.085016
  * (30) G. Beuf, Physical Review D 89(7) (2014).  DOI 10.1103/physrevd.89.074039.  URL http://dx.doi.org/10.1103/PhysRevD.89.074039
  * (31) E. Iancu, J. Madrigal, A. Mueller, G. Soyez, D. Triantafyllopoulos, Physics Letters B 744, 293–302 (2015).  DOI 10.1016/j.physletb.2015.03.068.  URL http://dx.doi.org/10.1016/j.physletb.2015.03.068
  * (32) E. Iancu, J. Madrigal, A. Mueller, G. Soyez, D. Triantafyllopoulos, Physics Letters B 750, 643–652 (2015).  DOI 10.1016/j.physletb.2015.09.071.  URL http://dx.doi.org/10.1016/j.physletb.2015.09.071
  * (33) Y. Hatta, E. Iancu, JHEP 08, 083 (2016).  DOI 10.1007/JHEP08(2016)083
  * (34) P. Korcyl, arXiv p. 2009.02045 (2020)
  * (35) S. Cali, K. Cichy, P. Korcyl, P. Kotko, K. Kutak, C. Marquet, (2021).  DOI 10.1140/epjc/s10052-021-09380-6
  * (36) L. McLerran, R. Venugopalan, Physical Review D 50(3), 2225–2233 (1994).  DOI 10.1103/physrevd.50.2225.  URL http://dx.doi.org/10.1103/PhysRevD.50.2225
  * (37) K.J. Golec-Biernat, M. Wusthoff, Phys. Rev. D59, 014017 (1998).  DOI 10.1103/PhysRevD.59.014017
  * (38) C. Marquet, E. Petreska, C. Roiesnel, JHEP 10, 065 (2016).  DOI 10.1007/JHEP10(2016)065
  * (39) H. Mäntysaari, B. Schenke, Phys. Rev. D 98, 034013 (2018).  DOI 10.1103/PhysRevD.98.034013.  URL https://link.aps.org/doi/10.1103/PhysRevD.98.034013
  * (40) T. Lappi, H. Mäntysaari, Nucl. Phys. A932, 69 (2014).  DOI 10.1016/j.nuclphysa.2014.07.009
  * (41) H. Mäntysaari, B. Schenke, Phys. Rev. D98(3), 034013 (2018).  DOI 10.1103/PhysRevD.98.034013
  * (42) B.G. Giraud, R.B. Peschanski, Acta Phys. Polon. B 37, 331 (2006)
  * (43) B.G. Giraud, R. Peschanski, Phys. Lett. B 760, 26 (2016).  DOI 10.1016/j.physletb.2016.06.033