# Unsupervised, Bottom-up Category Discovery for Symbol Grounding with a
Curious Robot

Catherine Henry  Department of Computer Science
Boise State University
Boise Idaho, 83725
<EMAIL_ADDRESS>Casey Kennington  Department of Computer
Science
Boise State University
Boise Idaho, 83725
<EMAIL_ADDRESS>
###### Abstract

Towards addressing the Symbol Grounding Problem and motivated by early
childhood language development, we leverage a robot which has been equipped
with an approximate model of curiosity with particular focus on bottom-up
building of unsupervised categories grounded in the physical world. That is,
rather than starting with a top-down symbol (e.g., a word referring to an
object) and providing meaning through the application of predetermined
samples, the robot autonomously and gradually breaks up its exploration space
into a series of increasingly specific unlabeled categories at which point an
external expert may optionally provide a symbol association. We extend prior
work by using a robot that can observe the visual world, introducing a higher
dimensional sensory space, and using a more generalizable method of category
building. Our experiments show that the robot learns categories based on
actions and what it visually observes, and that those categories can be
symbolically grounded into.

##  I Introduction

The symbol grounding problem posits that the meaning behind all words cannot
be entirely composed of other words (e.g., definitions); at some point
something else has to give meaning to a sufficient subset of words that are
grounded into the physical world [12, 35]. The meaning of such foundational
words is often derived from physical experience with concrete objects within
the physical world in a bottom-up conceptual learning progression [25]. In
other words, the learning of conceptual categories comes first and the symbol
(i.e., word label) comes later, whereas in most research on symbol grounding,
the symbol comes first. This unsupervised, bottom-up approach, where a
category develops independent of the symbol which is then later attached to
the category, follows the early language learning progression of human
children. According to Smith and Gasser [33], children live in a physical
world full of rich regularities that organize perception, action, and thought
and they interact in a social world to learn a shared linguistic, symbolic
communicative system. This means that physical embodiment is very likely a
prerequisite for symbol grounding (i.e., learning categories bottom-up) and,
therefore, language acquisition (see also [13, 7, 22, 10, 8]). We infer,
therefore, that the most likely platform for language acquisition that matches
closely to humans will be on robotic platforms because they are embodied and
operate in the same world that humans do.111This is in contrast to how
language models learn using a distributional representation of word meaning
largely from text or text-image pairs; the learning progression is the
opposite that of humans [32].

But what motivates children to explore and categorize their world? In this
paper, we focus on curiosity. As Ginsburg and Sylvia Opper [11] explained,
“curiosity is a function of the relation between the new object and the
individual’s previous experience.” As a child explores and engages with the
space around them, they are likely to encounter objects which may be entirely
novel or which they have not yet interacted with while a situational expert
(i.e., adult) is present to assign a symbol (i.e., a word) to the object of
interest–following the symbol grounding progression. For example, a child
picks up an un-before-seen object to see and touch it, but does not learn
until later when an adult utters spoon that there is a label for that
category.

In this paper, we leverage an existing model of curiosity from Oudeyer et al
[30, 29, 28] and a Cozmo robot which can physically move through a space and
explore the visual properties of objects contained within that space in order
to build unsupervised bottom-up categories. The physical movement of the
robot, per action/perception turn, consists of a single rotation and linear
movement forward or backward. These movements are an approximation of an
infant’s motor actions of looking and crawling. Exemplars (realized as
vectors) collected throughout the exploration process result in a
representation of the sensorimotor space $\mathbf{SM(t)}$ are composed of the
motor action and the sensory experience of the robot associated with each
movement. The sensory experience of the robot is a visual vector and can also
include supplementary positional information. As the robot explores, there is
a process in which categories are carved out of the sensorimotor space. The
robot’s autonomous, unsupervised discovery and categorization of these regions
of sensorimotor space is the primary focus of our work as it matches more
directly with the bottom-up learning progression of children and, we argue,
the symbol grounding problem of identifying categories before they are
labeled.

With this curious system, it is our intent to flip traditional methods of
grounding in that we begin with an existing, unlabeled, category which the
robot later associates with a symbol rather than starting with a symbol and
associating a relevant category (which constitutes most work on symbol
grounding). It is our hope that by enabling the robot to autonomously identify
unsupervised categories throughout the space without proactively or
immediately assigning labels, we are able to explore meaning with a bottom-up
perspective. Furthermore, this independent exploration of a space allows for
the reduced manual effort—perhaps even removal—of defining and building tasks
of varying complexity for the purpose of aggregating data to be used in
training systems. Our experiments build directly on and systematically extend
prior work until we test our hypothesis that the categories are developed
enough for symbol grounding by learning word-level classifiers for the robot’s
discovered categories.

##  II Related Work

Our work is partially related to a growing literature on few-shot learning in
robotic scenarios, including datasets [20] and methods. Early methods included
[36, 34] which used different Bayesian approaches for word learning. Seminal
work by [19] proposed a method to learn novel categories, focusing on
character recognition using few examples. More recently, [2] proposed a
Gaussian Mixture Model method for few-shot learning to mitigate catastrophic
forgetting, and [26] proposed online active control for lifelong learning of
object categories based on visual information. Finally, [21] used incremental
learning and active class selection to improve few-shot learning. Though
related to this literature, our research is fundamentally different in that we
are directly focusing on resolving an aspect of the symbol grounding problem
by bottom-up learning of unlabeled categories discovered as the robot
curiously explores the world. This brings to bear an enactive element of
learning in that the robot is actively exploring its physical space instead of
only being exposed to pre-collected data [14, 5]. In contrast, this related
work is focused on categories with top-down labels.

More directly related to our work is Intelligent Adaptive Curiosity (IAC),
which has been extended in a variety of ways. Aly et al [1] leverage one such
extension, RL-IAC [9], in that they model segmenting objects in an
unsupervised manner. Our work shares additional parallels with theirs as they
explore grounding spatial concepts and object categories through visual
perception.

Our work is also related to developmental robotics (see [8]) in that we use
similar methods to enable a robot to explore and build a degree of
understanding about its environment, but we deviate from the developmental
robotics agenda in that we focus on using these methods for ultimately
building concrete categories that form the basis of first language
acquisition. Thus our work is related to the long and still growing literature
on addressing Harnad’s Symbol Grounding Problem. Recent work showed that
abstract categories can be built on more basic, concrete categories, but the
categories were pre-determined [15], whereas in our work categories are
discovered autonomously by the robot as it explores.

##  III Method

In this section, we give a background of the work on curiosity that we build
on, then explain our system that uses curiosity to explore and learn about its
space.

###  III-A Background

The learning pipeline in our system includes technologies which are expanded
on in the following subsections. Figure 1 is a high level system overview
referencing each of these components in our model.

Figure 1: High level overview of the Explauto + Retico pipeline.
Initialization point is denoted by the star. In the Perception phase, the
Sensorimotor model is updated with the Motor Goal vs Actual Sensory Effect and
the Interest Model is updated with the Motor Goal + Predicted Sensory Effect
vs Motor Goal + Actual Sensory Effect. Production phase uses the InterestModel
to for the Motor Goal and the Sensorimotor Model to predict Sensory Effect.

#### Cozmo

The robot we use in our experiments is Cozmo.222Originally produced by Anki,
Cozmo is now Licensed by Digital Dream Labs. Cozmo is a small tabletop robot
which moves using a tread system and has a built-in camera capable of
streaming at 320x240 (QVGA) resolution at a rate of approximately 15 frames
per second. We leverage the Anki provided Python SDK for interacting directly
with Cozmo for sensory input and robot control.

#### Retico

We use the incremental framework Retico and the robot-ready spoken dialogue
system extensions [23, 17] to support the integration of Cozmo with our
system. While Retico is typically oriented to Spoken Dialogue Systems, low
granularity transfer of information (e.g., camera frames) allows for a
straightforward execution of and data pass-through to and from Explauto
functions. Furthermore, we are able to leverage pre-existing robot-ready
spoken dialogue system modules including object detection, and object
representation. Importantly, leveraging Retico’s framework for this work will
enable integration of capabilities that will allow the robot to learn object
category names (i.e., words) from a human interlocutor, which we leave for
future work.

#### Object Detection & Representation (v1): YOLOv4 & CLIP

YOLO (You only look once) is an object detection model designed for real-time
object detection tasks [6]. YOLO deviates other from other detection systems
by re-framing object detection as a single regression problem going straight
from image pixels to bounding box coordinates and class label probabilities.
YOLO allows us to quickly identify objects within an image and select observed
objects with the highest confidence score, in our case the robot only observes
a single object at a time in a given frame.333We are using the YOLOv4
implementation from https://github.com/yuto3o/yolox with pretrained weights
from https://github.com/AlexeyAB/darknet/releases We only use YOLO to identify
that an object exists and its bounding box; we ignore its label. YOLO is meant
to act as a proxy for a child’s capability to isolate individual objects, even
if their categories are unknown. Figure 2 is an example of how we use YOLO to
extract a sub-image with the object from a camera frame.

Figure 2: Using the YOLO bounding box to extract a sub-image of the object in
view. This sub-image is passed to the CLIP module. The object name assigned by
YOLO is dropped entirely.

CLIP is a multi-modal vision and language model that can be used for image-
text similarity and for zero-shot image classification [31]. We leverage the
image encoder to encode the sub-image from the YOLO bounding box with a
resulting feature vector of size 512. This vector, supplemented with
additional positional information including width, height, and object location
relative to the overall image, serves as the backbone representation of our
unlabeled categories. We refer to this vector as CLIP+.

YOLO is trained using supervised data with few object categories and CLIP is a
model that is highly trained on symbolic language data, meaning they are
somewhat at odds with our long-term research goals that adhere to child
language acquisition settings. YOLO and CLIP act as placeholders in our first
experiments to allow us us to systematically ensure that our pipeline has the
highest likelihood of success as we examine systematic changes from prior
work. For our final experiment, we use models for object detection and
representation that are modeled and trained using unsupervised methods,
explained presently.

#### Object Detection & Representation (v2): Segment Anything & DINOv2

The segment anything model (SAM) is an image segmentation model for object
segmentation tasks [18]. SAM is not like other methods (e.g., YOLO) which
identify objects with labels; rather, SAM allows us to identify objects within
an image and, in our case, only observe a single “salient” object at a time in
a given frame from the robot’s camera. SAM is meant to act as a proxy for a
child’s capability to identify unlabeled objects. Figure 3 is an example of
how we use SAM to isolate the object through masking.

DINOv2 is a vision model we use to represent objects segmented from SAM [27].
DINOv2 is a good choice for this step because the model training is largely
unsupervised, similar to how children learn about the visual world more
directly than CLIP. We leverage the DINOv2 small model to encode the sub-image
from the SAM segment with a resulting feature vector of size 384. This vector
serves as the backbone representation of our unlabeled categories.

Figure 3: Using the SAM segments to mask a sub-image of the object in view.
This masked image is passed to the DINOv2 module.

###  III-B Modeling Curiosity: Explauto

The Explauto framework is the driving force behind the curiosity in our
system. It is a framework for studying, modeling, and simulating curiosity-
driven learning and exploration in both simulated and physical robotic
systems, developed by the Inria FLOWERS research team [3, 24, 28]. Explauto
comes with several sensorimotor systems available including simulated
environments and interfaces to real robots.444Based on Dynamixel actuators
through the Pypot library The Cozmo robot is not natively integrated with
Explauto, for example it does not have positional information of an arm end-
effector. While we do retrieve the camera output directly from Cozmo, which we
task Explauto to use, the sensory data we use in our model requires
postprocessing of the camera output through object detection and
representation.

The Explauto architecture consists of three core processing levels which we
elaborate on in the System Overview below. While the Explauto framework was
designed to support both motor and goal ‘babbling’ as well as inverse and
forward predictions, our implementation focuses only on motor babbling with
forward predictions. As such, in the following sections we remain focused on
the functionality relevant to our experiments, though it is important to note
that this only captures one aspect of the overall Explauto system.

Figure 4: With each iteration the Learning Potential for the region updates
dynamically. Regions with high learning potential have brighter colors.

The Sensorimotor Model learns mappings between the robot motor actions and the
sensory effect they produce. This is an active, unsupervised model which is
trained iteratively on the sensorimotor experience. This model is also used in
performing the forward predictions given a motor command. The sensorimotor
model maintains a dataset which includes a list of motor commands and sensory
effects from each execution as well as a KD Tree of the motor commands.555The
KD Tree implementation is from
https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.html
Explauto provides several sensorimotor models including locally weighted
linear regression (LWLR), which we use in our implementation. The LWLR
sensorimotor model predicts the sensory effect $s$ given a motor command $m$
by computing a linear regression of k nearest neighbors of m. The KD Tree
supports a nearest neighbor lookup, which allows for indexing the motor and
sensory lists to extract the X, Y values used in this linear regression.

The Interest Space is the area of exploration from which goals are sampled;
for our purposes we are using motor babbling strategies so the interest space
consists of four possible motor actions $M$: degrees of rotation (left or
right) and mm linear travel (forward or backward).

Figure 5: Results from an execution of our baseline setup. In the grid plot on
the left, the sub-regions (unlabeled categories) have a very loose alignment
with the objects in the space. The polar plot on the right shows that recent
motor goals continue to be sampled in both object and non-object areas with a
heavier concentration of activity at the a negative mm Linear Distance.

Figure 6: Results from an execution of Experiment 1. The grid plot to the left
shows the sub-regions continue to loosely align with the objects in the space
and with fewer sub-regions developing in non-object areas than the Baseline in
Figure 6. The polar plot to the right shows recent motor goals are more
frequently around the objects with fewer goals being sampled in the large
whitespace zone.

The Interest Model implements the active exploration process and it is what
approximates curiosity in our learning system. It is responsible for selecting
motor goals from the Interest Space which improve the predictions of the
Sensorimotor model. The sensory results associated with the execution of the
motor goal are added to the experiential history stored in the Sensorimotor
Model. At each time step, the sampling procedure selects a motor action and
the sensorimotor model attempts to predict the sensory outcome. A ‘measure of
competence’ is calculated using the actual and predicted sensory outcomes and
the competencies for a region are used to calculate a region’s learning
progress. In our experiments, we begin with the Euclidean distance calculation
then introduce a cosine distance calculation as a part of our contribution.

The sampling procedure we use for exploring the interest space is Learning
Progress Maximization. This sampling procedure discretizes the interest space
using a KD Tree like structure,666adapted from scipy.spatial.kdtree where each
split is along alternating axes. A region is split when the number of
observations exceeds the configured maximum per region. There are various
methods provided for determining the value along which to split the region, in
our experiments we use the provided Best Interest split introduced in Baranes
and Oudeyer [4] as well as a Variance of Cosine Similarity split we introduce
below. A visual representation of the discretized interest space as the result
of region splitting can be seen in Figure 10 from the Experiment section
below; a corresponding KD-subtree can be seen in Figure 7.

Figure 7:  KD sub-tree example; the top node splits using an angle threshold
whereas the bottom two nodes split on linear travel; each records information
about angle and linear travel bounds, variance of the vectors (euclidean or
cosine distance), and learning potential.

###  III-C Symbol Grounding: Words as Classifiers

We use the Words-as-Classifiers (WAC) model [16] in the Symbol Grounding phase
of our pipeline. WAC grounds an unlabeled category represented by its own
classifier trained on positive and negative examples of sensory data from the
region (negative samples being randomly sampled from other regions). WAC is
well suited to our task, where unlabeled category has at most 30 observations,
as it has been shown to perform well on small amounts of training data.

##  IV Experiments

In this section, we explain our experiments that build on prior work by
incremental extensions. We first establish that Explauto can work with our
framework and robot using YOLO+CLIP, then extend the curiosity model to make
use of high-dimensional, visual vectors then use category learning (i.e.,
region splitting as explained above) using cosine similarity and the YOLO+CLIP
then, SAM+DINOv2 pipeline. We first explain the procedure and metrics we use
for each experiment (as they are the same for all three), then we explain the
three experiments, then we compare the results.

Figure 8: The Experiment space used for our three experiments. The objects
were located at approximately +55 degrees and and +/-180 degrees. The robot is
initialized at at the center of the space facing 0 degrees.

#### Procedure

All experiments took place in a constrained environment with a limited
selection of objects; we use a toy cat and an elephant (see Figure 8).777Our
choice of using only two objects in a limited space for our experiments is
important: we systematically establish that each of our extensions to prior
work show improvement before we move onto the next experiment. Scaling our
pipeline to work in noisier environments with a greater degree of movements
and sensorimotor input is planned for the near future. To counteract the
gradual aggregate of positional error due to a 2 degree tolerance in Cozmo’s
rotation functionality, the robot re-centers itself after every fifth
action.888On occasion, it is necessary to manually adjust a poor re-centering
job which can be done during a brief pause before execution continues. In
order to support this automated re-centering, the elephant object has a symbol
on the front which signifies to Cozmo that it can be used for docking. While
the symbol has an impact on the CLIP or DINOv2 vector produced, we do not
expect it to have a detrimental impact to the category development process.

The experiment space measures 762x508 millimeters and has white floor and
walls to reduce background complexity. This space allows for the robot to
rotate up to $+/-$180 degrees and travel in a line up to $+/-$80mm. The
objects can be placed anywhere in the experiment space, though it is
beneficial for the object which doubles as a re-centering agent to be along
the X or Y axis. While object locations remained static throughout the
experiments outlined in this paper, we observed consistent category
development no matter where objects were placed in the space.

Each experiment consisted of 300 Action/Perception turns which lasted
approximately three hours. After each action/perception turn, Cozmo returned
to center before executing the next turn. For prediction, the Sensorimotor
Model LWLR retrieved neighbors from the KD Tree. The Interest Model parameters
had a max region size of 30, sampling epsilon (likelihood of random sampling)
of 0.1, and Learning Progress Window Size (lookback window for calculating
learning progress) of 10.

Figure 9: Results from an execution of Experiment 2.

Figure 10: Results from an execution of Experiment 3.

#### Metrics

Throughout our experiments, we observe where in the Interest Space each action
is sampled from, how each action impacts the learning progress for a region,
and the value along which the axis is split during region splitting. For every
action taken there is a target (predicted) and reached (actual) sensory
outcome, as defined in Explauto. These sensory outcomes are used in
calculating ‘competence’ which is used in calculating the region’s learning
progress which is used in determining the next action to sample. As each
experiment progresses, this cycle results in the development of sub-regions
throughout the space which we refer to as unlabeled categories. For each
experiment, we capture these values in figures to better understand the
category boundaries.

The grid plot (Figure 6, left) provides a view into the region splitting and
action sampling mechanisms. Each point represents an action taken, with the
degree of rotation along the x axis and mm linear distance traveled (post
rotation) along the y axis. In the grid plot to the left the sub-regions
(unlabeled categories) have a much closer alignment to the objects in the
space. The polar plot to the right shows recent motor goals are largely
concentrated in areas where an object will be perceived. There are three
shaded areas marked on the plot which roughly denote actions which will result
in an object being in Cozmo’s FOV. The boundary of each region is marked by a
white line with a number reflecting the unlabeled category in the bottom
right. As the experiment progresses, the color of each region dynamically
updates to reflect the current learning progress as seen in Figure 4. The
polar plot (Figure 6, right) provides a view into where in the Interest Space
more recent actions have been sampled. The darker the point denotes a more
recent action, which allows us to observe if and when the exploration begins
to narrow in on objects in the space; shaded areas marked to signify the
objects. Comparisons can easily be made visually using these plots.

In our final evaluation, we compare the categories produced by the pipelines
that use YOLO+CLIP and SAM+DINOv2 by comparing the accuracy of unlabeled WAC
classifiers to identify objects within a category.

#### Baseline Setup: Extending Explauto with Cozmo

We first establish that Explauto can be integrated into the Retico pipeline
using Cozmo as the physical robot (i.e., the ability to physically move as
well as visually observe its environment—an extension from prior work) while
changing as little as possible in the Explauto framework. As such, we used
YOLO only in determining the binary presence of an object. This aligned most
closely with the original Playground Experiment for Explauto. For a given
motor action of rotation or linear travel, the sensory output is 1 when an
object was detected and -1 otherwise. The results of this experiment are shown
in Figures 4 and 6. The figures show that the novel sensorimotor space and
Cozmo robot are usable with Explauto, setting the stage for our experiments.

###  IV-A Experiment 1: Sensing a high-dimensional, visual sensory space with
YOLO+CLIP

In this experiment, we go beyond prior work by increasing the richness of the
sensory space by incorporating the full CLIP+ feature vector as the sensory
output for a given motor action (prior work only was concerned with binary
presence of an object, not the visual representation of an object). However,
prior work only used very small vectors (smaller than size 10), whereas here
we are working with much bigger vectors that can vary in size depending on the
model (i.e., CLIP produces vectors of size 512, DINOv2 of size 384).
Therefore, to support the higher dimensional vectors and variability of their
sizes, we alter the Competence calculation from Euclidean Distance to a
measure of Cosine Distance formalized as

$1-\cos({\bf t},{\bf e})={{\bf t}{\bf e}\over\|{\bf t}\|\|{\bf
e}\|}=\frac{\sum_{i=1}^{n}{{\bf t}_{i}{\bf e}_{i}}}{\sqrt{\sum_{i=1}^{n}{({\bf
t}_{i})^{2}}}\sqrt{\sum_{i=1}^{n}{({\bf e}_{i})^{2}}}}$ (1)

We bound this distance between (0,1) by $(1-e^{2*-cos\\_dist})$.

The results of this experiment can be viewed in Figure 6. The grid plot shows
that the sub-regions align with the objects in space with fewer sub-regions in
non-object areas than the Baseline. Moreover, the polar plot shows that motor
goals coalesce around objects rather than whitespace.

###  IV-B Experiment 2: Region Splitting using Cosine Similarity with
YOLO+CLIP

In this experiment, we build on Experiment 1 but instead of splitting by the
value which maximizes the difference of learning potential in the two sub-
regions, we explore a new region splitting mechanism. This new region
splitting mechanism divides the regions such that the variance of the cosine
similarity is minimal on each side, encouraging region splitting along visual
object boundaries. Importantly, this experiment results in the robot learning
specific categories of object features and location from the CLIP+ vector—the
autonomous bottom-up category discovery that can later be symbolically
grounded into.

The results of this experiment can be viewed in Figure 10. In this plot, we
observe much closer alignment to the objects in the space, and the motor goals
are highly concentrated in areas where an object will be perceived.

###  IV-C Experiment 3: Region Splitting using Cosine Similarity with
SAM+DINOv2

Figure 11: Given a sensory exemplar from an unlabeled category, WAC is able to
predict the correct category. On the left is the photo used to produce the
sensory exemplar, on the right are images from the region of the unlabeled
category. The upper correctly predicts Category 27 and the bottom Category 30.
Figure 12: Given a sensory exemplar from an unlabeled category, WAC is not
able to predict the correct category but the images from the region of the
predicted unlabeled category are nearly identical to the sensory exemplar and
the categories shown in the associated grid plot cover a similar range of
objects (seen in 10. We consider this to be correct on a technicality and
liken it to human error. The upper row predicted Category 23 for an exemplar
from Category 22 and the bottom predicted Category 17 for an exemplar from 8.

Figure 13: Given a sensory exemplar from an unlabeled category, WAC is not
able to predict the correct category. The first row is an edge case where the
predicted unlabeled category 20, for an exemplar from 28, consists of both cat
and “whitespace” observations (“whitespace” is shown by images with SAM
segment overlays) and the region seen in the associated grid plot does not
cover a similar range of objects as 28 (seen in 10). The bottom row, which
predicted Category 17 for exemplar 20, captures the edge case where there was
a limitation in SAM and a cat in view was not successfully segmented and/or
isolated resulting in a false “whitespace” observation. sensory exemplar of
$\mathbf{[-1]*384}$.

Figure 14: SAM having successfully segmented the wall and floor, but missing
the cat. There is also a fair amount of segments which don’t appear meaningful
though may be due to light bounce or noise in the image.

In this experiment we build on Experiment 1 and 2 by replacing the object
detection and representation from YOLO+CLIP to SAM+DINOv2 to move closer to
our goal of learning in a noisier environment. Because SAM identifies many
segments for a given frame, it was necessary to implement a method to reduce
the noise which would occur if non-object segments were isolated and passed
down the pipeline despite an object being in view. Such noise results in
inconsistent exemplars for a region and have significant negative downstream
impact to the ability to correctly assess novelty. To prevent this noise we
only use an isolated segment if the average pixel color is above 191.25. This
threshold value was determined by observing the average pixel value for
segments consisting of only the whitespace background of our environment which
does have a reddish hue from the Cozmo camera. Additionally, as shown in 3 the
extracted image is no longer a tight crop around the object.

The results of this experiment can be viewed in Figure 10. The grid-plot shows
a similar precision in constructing category boundaries as seen in Experiment
2. The polar plot, however, deviates from the other experiments in that
actions were far less clustered in the observation space. Instead, actions
were sampled more broadly throughout the observation space, even late into the
execution, with a slightly heavier concentration in areas where an object will
be perceived.

###  IV-D Discussion

The grid plots from the Baseline and Experiment 1 both show a very loose
relationship between region and objects in the space, as expected, but they
both establish that the framework with the added functionality works in a
similar way as prior work. Notably, however, the region splitting
implementation introduced in Experiment 2 results in region development with
far more precision along object boundaries. Improving this relationship
between regions and objects in space allows us to treat the regions as
category boundaries and consequently leverage the sensory data associated with
each region as a representation of that category. The grid plot in Experiment
3 shows similar precision in category boundaries as those seen in Experiment
2.

The polar plots from each experiment show the points from more recent
observations as a darker shade of black. From this, we can determine that
Experiment 2 increasingly samples motor goals which result in the observation
of an object. This is promising behavior as we are looking to capture rich
sensory representations of the unlabeled categories and there is little to no
benefit in collecting more non-object data points than is necessary to
determine the object boundary. It is interesting to observe that Experiment 3
does not appear to have the same level of honing in on actions resulting in
object observation.

###  IV-E Final Evaluation: Symbol Grounding into Categories

We evaluate the quality of the unlabeled categories from Experiments 2 and 3
using WAC. The goal is for WAC to predict the category name given a sensory
exemplar from the Interest Space. Each experiment resulted in 16 categories
(making a random baseline of 6.3%); each category had a range of 8 and 30
samples, average of 18.75. We randomly sampled 80% from each for training and
20% for testing. To train WAC, we train a binary logistic regression
classifier for each category using observations from a specific category as
positive examples, and randomly sample observations from neighboring
categories as negative samples. We report the number of times WAC was able to
accurately predict a category both when using exemplars from Experiment 2 in
training and test, and separately, Experiment 3. This step is crucial because
if classifiers can distinguish between the categories, then that means they
are ready to have a word label attached to them, completing the word-level
symbol grounding for that category (though at this step, the process is not
interactive; we leave that for future work).

Table I shows the results using the categories from Experiment 2. There were 5
correct Category predictions which require no explanation. The rows marked SC
denote scenarios where the actual and predicted categories were visually
similar and we consider them to be “Similar Categories”. The rows marked YOLO
would appear incorrect at a glance, but the point of failure is YOLO rather
than Explauto.

Actual Category | Predicted Category | Manual Explanation
---|---|---
5 | 24 | SC
3 | 3 |
10 | 24 | SC
15 | 18 | SC
11 | 24 | SC
28 | 24 | YOLO
18 | 17 | SC
6 | 24 | SC
30 | 30 |
8 | 8 |
29 | 29 |
21 | 24 | YOLO
26 | 24 | YOLO
24 | 24 |
17 | 18 | SC
22 | 29 | SC

TABLE I: Results from WAC predictions trained on the Unlabeled Categories from
YOLO+CLIP in Experiment 2.

Table II shows the results using the categories from Experiment 3. There were
6 correct Category predictions which require no explanation. The rows marked
SAM would appear incorrect at a glance but the point of failure is SAM rather
than Explauto. The rows marked IPCF, we theorize, are incorrect due to “in
progress category formation”. Categories 28 and 20 show a split along
object/whitespace boundaries, category 28 due to natural object/whitespace
boundary and category 20 due to SAM failing to segment objects at a distance.

As evidenced by the manual explanation labels, evaluation with WAC is not yet
a fully automated process. This is further demonstrated by Figure 12 where a
sensory exemplar from a particular category is predicted as being from a
different category. Looking at the images of the exemplar vs those
constituting the predicted category, they are nearly identical. Moreover, the
grid plot from Figure 10 shows that the two categories are visually similar to
one another in the Interest Space, both overlapping where the similar features
in the exploration space, elephant (row 1) or whitespace (row 2). Because the
predicted vs actual categories are so similar, so much so that it can be
likened to human error, we consider this to be a technically correct
prediction. Figure 14 has examples of incorrect predictions, row 1 due to
underdeveloped categories and row 2 due to inconsistencies in SAM
segmentation.

Taken together, these results show that both pipelines work to allow the robot
to curiously explore its space and arrive at unsupervised object categories.
In the YOLO+CLIP pipeline, the object categories tended to be cleaner (as
expected), but the noisier and more theoretically pleasing SAM+DINOv2 pipeline
worked well; we conjecture that it simply needs more exposure to the objects
to learn more rigid categories. This also highlights the robustness of our
pipeline; theoretically, any kind of object detection model and object
representation model could be used because the cosine distance calculation is
agnostic to the vector size from the representation model.

Figure 11 captures examples where WAC was able to correctly predict the
unlabeled category. The first row correctly predicts a category which appears
to capture cat-ness at a nearer proximity to Cozmo, as seen in the grid plot
in Figure 10. The second row appears to capture elephant-ness up close, also
as seen in the grid plot.

Actual Category | Predicted Category | Manual Explanation
---|---|---
8 | 17 | SC
30 | 30 |
28 | 20 | IPCF
3 | 3 |
20 | 17 | SAM
23 | 17 | IPCF
13 | 13 |
29 | 17 | IPCF
17 | 17 |
10 | 10 |
16 | 17 | SC
5 | 17 | SC
7 | 17 | SC
27 | 27 |
22 | 23 | SC
15 | 17 | SC

TABLE II: Results from WAC predictions trained on the Unlabeled Categories
from SAM+DINOv2 in Experiment 3.

##  V Limitations

#### SAM

In the current implementation of Experiment 3 we are able to iterate through
the segments and isolate only those which captures the object in frame (as
opposed to a segment for wall, floor, or miscellaneous noise). As is, this can
not be scaled to a real-world scenario outside of the controlled experiment
space because the average pixel threshold is meaningless outside of a known
and consistent background. Working with unlabeled segments trained only using
image data is applicable to our long-term goals, however it does come with its
challenges. The foremost is that there are limited ways to ensure the same
object or feature within the exploration space for a given motor action is
consistently segmented, especially in an object rich space. Some of the
complexities of SAM segmentation can be seen in Figure 14.

#### Robot Movement Drift

For the purposes of our experiments, we always re-centered Cozmo from a
starting position so we could systematically compare our model’s decision as
to which action to take next. Cozmo unfortunately drifts when it turns (within
a tolerance of 2 degrees), requiring us to use the symbol (see Figure 3) to
re-orient Cozmo. This does not seem to interfere with our results, but ideally
this would not be required (though it is far better than manually re-centering
Cozmo after each step). This constraint is only applicable, however, to our
experimental setup; a setting where Cozmo can roam freely in a space would not
require strict re-centering.

##  VI Conclusion and Future Work

Through visual analysis of the output plots as well as our final evaluation
with WAC, our experiments illustrate the robot is able to autonomously explore
the space, constructing unsupervised, unlabeled categories in a bottom-up
progression. It is important to emphasize here that these categories have
never had a label manually assigned; rather, they develop bottom-up based on
the robot’s Sensorimotor experience and we simply used WAC to see if it could
distinguish between the (still unlabeled) categories. Categorization is a
crucial step of the language learning process and, more generally cognition
because “to cognize is to categorize” [13].

In future work, we will address the limitations explained above. We will also
include scenes with more objects, more complex objects, and a larger space for
Cozmo to explore. We will also include additional sensorimotor information as
well as increasing the actions that Cozmo can take.

## References

  * Aly et al. [2017] Amir Aly, Akira Taniguchi, and Tadahiro Taniguchi.  A generative framework for multimodal learning of spatial concepts and object categories: An unsupervised.  In _2017 Joint IEEE International Conference on Development and Learning and Epigenetic Robotics (ICDL-EpiRob)_ , pages 376–383. unknown, September 2017.
  * [2] Ali Ayub and Carter Fendley.  Few-Shot continual active learning by a robot.  In _36th Conference on Neural Information Processing Systems_.
  * Baranes and Oudeyer [2009] Adrien Baranes and Pierre-Yves Oudeyer.  R-IAC: Robust intrinsically motivated exploration and active learning.  _IEEE Trans. Auton. Ment. Dev._ , 1(3):155–169, October 2009.
  * Baranes and Oudeyer [2013] Adrien Baranes and Pierre-Yves Oudeyer.  Active learning of inverse models with intrinsically motivated goal exploration in robots.  _Rob. Auton. Syst._ , 61(1):49–73, January 2013.
  * Bisk et al. [2020] Yonatan Bisk, Ari Holtzman, Jesse Thomason, Jacob Andreas, Yoshua Bengio, Joyce Chai, Mirella Lapata, Angeliki Lazaridou, Jonathan May, Aleksandr Nisnevich, Nicolas Pinto, and Joseph Turian.  Experience grounds language.  _arXiv_ , 2020.
  * Bochkovskiy et al. [2020] Alexey Bochkovskiy, Chien-Yao Wang, and Hong-Yuan Mark Liao.  YOLOv4: Optimal speed and accuracy of object detection.  April 2020.
  * Bruner [1983] Jerome Bruner.  _Child’s Talk: Learning to Use Language_.  W. W. Norton & Company, 1983.
  * Cangelosi and Schlesinger [2015] Angelo Cangelosi and Matthew Schlesinger.  _Developmental robotics: From babies to robots_.  MIT press, 2015.
  * Craye et al. [2016] Celine Craye, David Filliat, and Jean-Francois Goudou.  RL-IAC: An exploration policy for online saliency learning on an autonomous mobile robot.  In _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, October 2016.
  * Eimas et al. [1987] Peter D Eimas, Joanne L Miller, and Peter W Jusczyk.  On infant speech perception and the acquisition of language.  In Stevan Harnad, editor, _Categorical perception: The groundwork of cognition , (pp_ , volume 599, pages 161–195. Cambridge University Press, x, New York, NY, US, 1987.
  * Ginsburg and Sylvia Opper [1988] Ginsburg and Herbert Sylvia Opper.  _Piaget’s Theory of Intellectual Development_.  International Psychotherapy Institute E-Books, 3rd edition, 1988.
  * Harnad [1990] Stevan Harnad.  The symbol grounding problem.  _Physica D_ , 42(1-3):335–346, 1990.
  * Harnad [2017] Stevan Harnad.  To cognize is to categorize: Cognition is categorization.  In _Handbook of Categorization in Cognitive Science_ , pages 21–54. 2017.
  * Johnson [2008] Mark Johnson.  _The meaning of the body: Aesthetics of human understanding_.  University of Chicago Press, 2008.
  * Kennington and Natouf [2022] Casey Kennington and Osama Natouf.  The symbol grounding problem re-framed as concreteness-abstractness learned through spoken interaction.  In _Proceedings of the 26th Workshop on the Semantics and Pragmatics of Dialogue - Full Papers_ , August 2022.
  * Kennington and Schlangen [2015] Casey Kennington and David Schlangen.  Simple learning and compositional application of perceptually grounded word meanings for incremental reference resolution.  In _Proceedings of the 53rd Annual Meeting of the Association for Computational Linguistics and the 7th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 292–301, Beijing, China, 2015. Association for Computational Linguistics.
  * Kennington et al. [2020] Casey Kennington, Daniele Moro, Lucas Marchand, Jake Carns, and David McNeill.  rrSDS: Towards a robot-ready spoken dialogue system.  In _Proceedings of the 21st Annual SIGdial Meeting on Discourse and Dialogue_ , Virtual, 2020. Association for Computational Linguistics.
  * [18] Alexander Kirillov, Eric Mintun, Nikhila Ravi, Hanzi Mao, Chloe Rolland, Laura Gustafson, Tete Xiao, Spencer Whitehead, Alexander C Berg, Wan-Yen Lo, Piotr Dollár, and Ross Girshick.  Segment anything.  _arXiv:2304. 02643_.
  * Lake et al. [2015] Brenden M Lake, Ruslan Salakhutdinov, and Joshua B Tenenbaum.  Human-level concept learning through probabilistic program induction.  _Science_ , 350(6266):1332–1338, December 2015.
  * Lomonaco and Maltoni [2017] Vincenzo Lomonaco and Davide Maltoni.  CORe50: a new dataset and benchmark for continuous object recognition.  In Sergey Levine, Vincent Vanhoucke, and Ken Goldberg, editors, _Proceedings of the 1st Annual Conference on Robot Learning_ , volume 78 of _Proceedings of Machine Learning Research_ , pages 17–26. PMLR, 2017.
  * [21] Christopher McClurg, Ali Ayub, Harsh Tyagi, Sarah M Rajtmajer, and Alan R Wagner.  Active class selection for Few-Shot Class-Incremental learning.  In _Conference on Lifelong Learning Agents_.
  * McCune [2008] Lorraine McCune.  _How Children Learn to Learn Language_.  Oxford University Press, 2008.
  * Michael and Möller [2019] Thilo Michael and Sebastian Möller.  ReTiCo: An open-source framework for modeling real-time conversations in spoken dialogue systems.  In _Tagungsband der 30. Konferenz Elektronische Sprachsignalverarbeitung 2019_ , ESSV, pages 134–140, Dresden, March 2019. TUDpress, Dresden.
  * Moulin-Frier et al. [2014] Clement Moulin-Frier, Pierre Rouanet, and Pierre-Yves Oudeyer.  Explauto: an open-source python library to study autonomous exploration in developmental robotics.  In _4th International Conference on Development and Learning and on Epigenetic Robotics_. IEEE, October 2014.
  * Muraki et al. [2022] Emiko J Muraki, Israa A Siddiqui, and Penny M Pexman.  Quantifying children’s sensorimotor experience: Child body-object interaction ratings for 3359 english words.  _Behav. Res. Methods_ , 54(6):2864–2877, December 2022.
  * Nie et al. [2023] Xiangli Nie, Zhiguang Deng, Mingdong He, Mingyu Fan, and Zheng Tang.  Online active continual learning for robotic lifelong object recognition.  _IEEE Trans Neural Netw Learn Syst_ , PP, September 2023.
  * Oquab et al. [2023] Maxime Oquab, Timothée Darcet, Theo Moutakanni, Huy V Vo, Marc Szafraniec, Vasil Khalidov, Pierre Fernandez, Daniel Haziza, Francisco Massa, Alaaeldin El-Nouby, Russell Howes, Po-Yao Huang, Hu Xu, Vasu Sharma, Shang-Wen Li, Wojciech Galuba, Mike Rabbat, Mido Assran, Nicolas Ballas, Gabriel Synnaeve, Ishan Misra, Herve Jegou, Julien Mairal, Patrick Labatut, Armand Joulin, and Piotr Bojanowski.  DINOv2: Learning robust visual features without supervision, 2023.
  * Oudeyer et al. [2005] P Oudeyer, V Hafner, and Andrew Whyte.  The playground experiment: Task-Independent development of a curious robot.  2005\.
  * Oudeyer and Kaplan [2006] Pierre-Yves Oudeyer and Frédéric Kaplan.  Discovering communication.  _Conn. Sci._ , 18(2):189–206, June 2006.
  * Oudeyer et al. [2007] Pierre-Yves Oudeyer, Frdric Kaplan, and Verena V Hafner.  Intrinsic motivation systems for autonomous mental development.  _IEEE Trans. Evol. Comput._ , 11(2):265–286, April 2007.
  * Radford et al. [2021] Alec Radford, Jong Wook Kim, Chris Hallacy, Aditya Ramesh, Gabriel Goh, Sandhini Agarwal, Girish Sastry, Amanda Askell, Pamela Mishkin, Jack Clark, Gretchen Krueger, and Ilya Sutskever.  Learning transferable visual models from natural language supervision.  February 2021.
  * Schlangen [2023] David Schlangen.  What a situated Language-Using agent must be able to do: A Top-Down analysis.  February 2023.
  * Smith and Gasser [2005] Linda Smith and Michael Gasser.  The development of embodied cognition: Six lessons from babies.  _Artif. Life_ , (11):13–29, 2005.
  * Squire and Levinson [2007] Kevin M Squire and Stephen E Levinson.  HMM-Based concept learning for a mobile robot.  _IEEE Trans. Evol. Comput._ , 11(2):199–212, April 2007.
  * Vincent-Lamarre et al. [2016] Philippe Vincent-Lamarre, Alexandre Blondin Massé, Marcos Lopes, Mélanie Lord, Odile Marcotte, and Stevan Harnad.  The latent structure of dictionaries.  _Top. Cogn. Sci._ , 8(3):625–659, July 2016.
  * Xu and Tenenbaum [2007] Fei Xu and Joshua B Tenenbaum.  Word learning as bayesian inference.  _Psychol. Rev._ , 114(2):245–272, April 2007.


# Recombination mechanism for $\bm{D^{0}}$-meson production and
$\bm{D^{0}\\!-\\!\overline{D^{0}}}$ production asymmetry in the LHCb
$\bm{p\\!+\\!\\!^{20}\\!N\\!e}$ fixed-target experiment

Rafał Maciuła<EMAIL_ADDRESS>Institute of Nuclear Physics, Polish
Academy of Sciences, ul. Radzikowskiego 152, PL-31-342 Kraków, Poland  Antoni
Szczurek<EMAIL_ADDRESS>Institute of Nuclear Physics, Polish
Academy of Sciences, ul. Radzikowskiego 152, PL-31-342 Kraków, Poland  College
of Natural Sciences, Institute of Physics, University of Rzeszów, ul. Pigonia
1, PL-35-310 Rzeszów, Poland

###### Abstract

We discuss production of neutral $D$ mesons in proton-proton collisions at the
LHC (fixed target mode) in the framework of the BJM recombination model. We
present rapidity and transverse momentum distributions of $D$ mesons and
compare the recombination contribution to the dominant gluon-gluon fusion
mechanism. Both the direct production, as dictated by the matrix element, and
fragmentation of the associated $c$ or $\bar{c}$ are included. The latter
mechanism generates $D$ mesons with smaller rapidities than those produced
directly. We calculate the $D^{0}+{\overline{D}^{0}}$ meson distributions
relevant for fixed target $p\\!+\\!\\!^{4}\\!H\\!e$ collisions at $\sqrt{s}$ =
86.6 GeV as well as for $p\\!+\\!\\!^{20}\\!N\\!e$ collisions at $\sqrt{s}$ =
69 GeV. The recombination component is consistent with the LHCb data and in
addition results in production asymmetry. The asymmetries in
$D^{0}\\!-\\!\overline{D^{0}}$ production as a function of rapidity and
transverse momentum are shown and the cancellation of terms for direct
production and associated $c/{\bar{c}}$ fragmentation is discussed.

##  I Introduction

The production mechanism of midrapidity $D$ mesons is relatively well known.
At low energy the quark-antiquark annihilation and gluon-gluon fusion must be
included in the leading order calculation. At higher energy the gluon-gluon
fusion is the dominant mechanism. We calculate the latter contribution in the
$k_{T}$-factorization approach. As already discussed in the literaturte such
an approach leads to a good description of single $D$ meson production data at
high energies, including RHIC Maciula:2015kea  and the LHC Maciula:2013wg ;
Maciula:2019izq , as well as of correlation observables for $D{\bar{D}}$-pair
Maciula:2013wg ; Karpishkov:2016hnx  and $DD$-pair production Maciula:2013kd ;
vanHameren:2014ava . Recently, it has been also applied for open charm meson
production at low energies within the fixed-target mode of the LHCb experiment
Maciula:2020cfy ; Maciula:2021orz .

On the other hand the mechanism of forward/backward production of $D$ mesons
is rather poorly understood. Two mechanisms are usually considered in this
context in phenomenological applications: (a) mechanism related to the knock-
out of preexisting intrinsic $c/\bar{c}$ from the nucleon BHPS1980 ;
Maciula:2020dxv ; Maciula:2021orz ; (b) mechanism of recombination of
$q/\bar{q}$ with $\bar{c}/c$ as proposed by Braaten-Jia-Mechen BJM2002a ;
BJM2002b ; BJM2002c . We shall call the first mechanism as intrinsic charm
(IC) mechanism and the second as BJM recombination.

As discussed e.g. in Refs. Maciula:2020dxv ; Maciula:2021orz  the Brodsky et
al. approach to the intrinsic charm BHPS1980  may lead to the production of
large rapidity $D$ mesons and consequently to large rapidity
neutrino/antineutrino production Maciula:2020dxv . The BHPS IC was included
recently into a global parton analysis Hou:2017khm . The original BHPS
formulation leads to $c_{IC}(x)=\bar{c}_{IC}(x)$. This mechanism was discussed
recently in the context of fixed target experiments Maciula:2021orz  as well
as of high-energy (anti)neutrinos observed recently by the IceCube laboratory
Goncalves:2021yvw  at the South Pole. There only upper limit for the BHPS
intrinsic charm can be obtained Goncalves:2021yvw  as there can be
extragalatic sources of very high (anti)neutrinos IceCube:2020wum . Another
mechanism of internal charm-anticharm component in the nucleon is related to
the meson cloud model Melnitchouk:1997ig ; Cazaroto:2013wy . In contrast to
the BHPS model it generates asymmetric intrinsic $c/{\bar{c}}$ distribution in
the nucleon. Similar model with light meson-baryon Fock components in the
nucleon was able to explain $\bar{d}-\bar{u}$ asymmetry in the nucleon, see
e.g. Holtmann:1996be , not possible to be explained in perturbative
calculations. As will be discussed here the $\bar{d}-\bar{u}$ asymmetry in the
nucleon, probably of the meson cloud origin, may have some influence on
$D-\bar{D}$ production asymmetries obtained in the recombination model. In our
calculations here we shall include distributions with ${\bar{d}}-{\bar{u}}$
asymmetry.

Here we rather concentrate on the recombination mechanism which may be
considered as an unwanted background for the large-$x$ intrinsic charm
contribution (which we will call here BHPS model for brevity). In this
mechanism there are two-body final states such as $D_{i}$ and $c$ or
${\bar{D}}_{i}$ and $\bar{c}$. We shall consider production of different
pseudoscalar and vector $D$ mesons as well as feed down due to decays of
vector $D^{*}$ mesons into measured pseudoscalar $D$ mesons. Somewhat
different model was used to explain asymmetries in $D$ meson production in
Rapp:2003wn .

The latter approach is based on the formalism proposed long ago by Das and Hwa
Das:1977cp . In this approach the transverse momentum of $D$ mesons is not
explicit in contrast to the BJM model. Stll another recombination approach was
presented in Berezhnoy:2000ji , but in the context of photoproduction. All the
previous studies of $D$ meson recombination concentrated on low-energy $pp$,
$\pi^{-}p$ or $\gamma p$ processes. According to our knowledge recombination
was not studied for higher energies. In Ref. Maciula:2017wov  a subleading
fragmentation ${\bar{d}}\to D^{+},d\to D^{-}$ was discussed to explain
observed asymmetry in $D^{+}/D^{-}$ production as observed in the LHCb
experiment LHCb:2012fb .

In principle, the asymmetry in $D$ meson production can be obtained also in
the Lund string approach Norrbin:2000zc . However, we shall not discuss it in
the present paper.

In the present paper we wish to estimate the contribution of the BJM mechanism
for fixed target experiment(s) at the LHC.

##  II Details of the model calculations

In the present study we take into consideration three different production
mechanisms of charm, including: a) the standard (and usually considered as a
leading) QCD mechanism of gluon-gluon fusion: $gg\to c\bar{c}$ ; b) the
mechanism driven by the intrinsic charm component of proton: $gc\to gc$; and
c) the recombination mechanism: $gq\to Dc$. Calculations of the first two
contributions are performed following our previous study in Ref.
Maciula:2021orz . Some of corresponding histograms shown there are repeated
here. What is completely new in the present study is the BJM recombination
component which was not considered in this context in the past.

###  II.1 The standard QCD mechanism for charm production

Figure 1:  A diagram of the standard QCD mechanism of charm production in the
$k_{T}$-factorization approach driven by the fusion of two off-shell gluons.

As was already mentioned here we follow the theoretical formalism for the
calculation of the $c\bar{c}$-pair production in the $k_{T}$-factorization
approach kTfactorization , as adopted and discussed in the context of the LHCb
fixed-target charm data in Ref. Maciula:2020cfy . In this framework the
transverse momenta $k_{t}$’s (or virtualities) of both partons entering the
hard process are taken into account, both in the matrix elements and in the
parton distribution functions. Emission of the initial state partons is
encoded in the transverse-momentum-dependent (unintegrated) PDFs (uPDFs). In
the case of charm flavour production the parton-level cross section is usually
calculated via the $2\to 2$ leading-order $g^{*}g^{*}\to c\bar{c}$ fusion
mechanism with off-shell initial state gluons that is the dominant process at
high energies (see Fig. 1). Even at lower energies as long as small transverse
momenta and not extremely backward/forward rapidities are considered the
$q^{*}\bar{q}^{*}\to c\bar{c}$ mechanism remains subleading. Then the hadron-
level differential cross section for the $c\bar{c}$-pair production, formally
at leading-order, reads:

$\displaystyle\frac{d\sigma(pp\to
c\bar{c}\,X)}{dy_{1}dy_{2}d^{2}p_{1,t}d^{2}p_{2,t}}$ $\displaystyle=$
$\displaystyle\int\frac{d^{2}k_{1,t}}{\pi}\frac{d^{2}k_{2,t}}{\pi}\frac{1}{16\pi^{2}(x_{1}x_{2}s)^{2}}\;\overline{|{\cal
M}^{\mathrm{off-shell}}_{g^{*}g^{*}\to c\bar{c}}|^{2}}$
$\displaystyle\times\;\delta^{2}\left(\vec{k}_{1,t}+\vec{k}_{2,t}-\vec{p}_{1,t}-\vec{p}_{2,t}\right)\;{\cal
F}_{g}(x_{1},k_{1,t}^{2},\mu_{F}^{2})\;{\cal
F}_{g}(x_{2},k_{2,t}^{2},\mu_{F}^{2})\;,$

where ${\cal F}_{g}(x_{1},k_{1,t}^{2},\mu_{F}^{2})$ and ${\cal
F}_{g}(x_{2},k_{2,t}^{2},\mu_{F}^{2})$ are the gluon uPDFs for both colliding
hadrons and ${\cal M}^{\mathrm{off-shell}}_{g^{*}g^{*}\to c\bar{c}}$ is the
off-shell matrix element for the hard subprocess. The gluon uPDF depends on
gluon longitudinal momentum fraction $x$, transverse momentum squared
$k_{t}^{2}$ of the gluons entering the hard process, and in general also on a
(factorization) scale of the hard process $\mu_{F}^{2}$. They must be
evaluated at longitudinal momentum fractions
$x_{1}=\frac{m_{1,t}}{\sqrt{s}}\exp(y_{1})+\frac{m_{2,t}}{\sqrt{s}}\exp(y_{2})$,
and
$x_{2}=\frac{m_{1,t}}{\sqrt{s}}\exp(-y_{1})+\frac{m_{2,t}}{\sqrt{s}}\exp(-y_{2})$,
where $m_{i,t}=\sqrt{p_{i,t}^{2}+m_{c}^{2}}$ is the quark/antiquark transverse
mass.

As we have carefully discussed in Ref. Maciula:2019izq , there is a direct
relation between a resummation present in uPDFs in the transverse momentum
dependent factorization and a parton shower in the collinear framework. In
most uPDF cases the off-shell gluon can be produced either from gluon or
quark, therefore, in the $k_{T}$-factorization all channels driven by
$gg,q\bar{q}$ and even by $qg$ initial states are open already at leading-
order (in contrast to the collinear factorization).

The kinematical configuration of the fixed-target LHCb experiment corresponds
to the region where in principle the Catani-Ciafaloni-Fiorani-Marchesini
(CCFM) CCFM  evolution equation is legitimate for any pQCD theoretical
calculations and could, in principle, be used to describe the dynamics behind
the mechanisms of e.g. open charm meson production. In the numerical
calculations below we follow the conclusions from Refs. Maciula:2020cfy ;
Maciula:2021orz  and apply the JH-2013-set2 gluon uPDFs Hautmann:2013tba  that
are determined from high-precision DIS measurements. As a default set in the
numerical calculations we take the renormalization scale
$\mu^{2}=\mu_{R}^{2}=\sum_{i=1}^{n}\frac{m^{2}_{it}}{n}$ (averaged transverse
mass of the given final state) and the charm quark mass $m_{c}=1.5$ GeV. The
strong-coupling constant $\alpha_{s}(\mu_{R}^{2})$ at next-to-next-to-leading-
order is taken from the CT14nnloIC PDF routines. The CCFM uPDFs here are taken
at a rather untypical value of the factorization scale
$\mu_{F}^{2}=M_{c\bar{c}}^{2}+P_{T}^{2}$, where $M_{c\bar{c}}$ and $P_{T}$ are
the $c\bar{c}$-invariant mass (or energy of the scattering subprocess) and the
transverse momentum of $c\bar{c}$-pair (or the incoming off-shell gluon pair).
This unusual definition has to be applied as a consequence of the CCFM
evolution algorithm Hautmann:2013tba .

###  II.2 The Intrinsic Charm induced component

Figure 2:  A diagrammatic representation of the intrinsic charm mechanism of
charm production within the hybrid model with the off-shell gluon and the on-
shell charm quark in the initial state.

The intrinsic charm contribution to charm production cross section (see Fig.
2) is obtained within the hybrid theoretical model discussed by us in detail
in Ref. Maciula:2020dxv . The LHCb fixed-target configuration allows to
explore the charm cross section in the backward rapidity direction where an
asymmetric kinematical configurations are selected. Thus in the basic $gc\to
gc$ reaction the gluon PDF and the intrinsic charm PDF are simultaneously
probed at different longitudinal momentum fractions - rather intermediate for
the gluon and large for the charm quark.

Within the asymmetric kinematic configuration $x_{1}\ll x_{2}$ the cross
section for the processes under consideration can be calculated in the so-
called hybrid factorization model motivated by the work in Ref. Deak:2009xt .
In this framework the small- or intermediate-$x$ gluon is taken to be off mass
shell and the differential cross section e.g. for $pp\to gcX$ via $g^{*}c\to
gc$ mechanism reads:

$\displaystyle d\sigma_{pp\to gcX}=\int d^{2}k_{t}\int\frac{dx_{1}}{x_{1}}\int
dx_{2}\;{\cal
F}_{g^{*}}(x_{1},k_{t}^{2},\mu^{2})\;c(x_{2},\mu^{2})\;d\hat{\sigma}_{g^{*}c\to
gc}\;,$ (2)

where ${\cal F}_{g^{*}}(x_{1},k_{t}^{2},\mu^{2})$ is the unintegrated gluon
distribution in one proton and $c(x_{2},\mu^{2})$ a collinear PDF in the
second one. The $d\hat{\sigma}_{g^{*}c\to gc}$ is the hard partonic cross
section obtained from a gauge invariant tree-level off-shell amplitude. A
derivation of the hybrid factorization from the dilute limit of the Color
Glass Condensate approach can be found e.g. in Ref. Dumitru:2005gt  (see also
Ref. Kotko:2015ura ). The relevant cross sections are calculated with the help
of the KaTie Monte Carlo generator vanHameren:2016kkz . There the initial
state quarks (including heavy quarks) can be treated as a massless partons
only.

Working with minijets (jets with transverse momentum of the order of a few
GeV) requires a phenomenologically motivated regularization of the cross
sections. Here we follow the minijet model Sjostrand:1987su  adopted e.g. in
Pythia Monte Carlo generator, where a special suppression factor is introduced
at the cross section level Sjostrand:2014zea :

$F(p_{t})=\frac{p_{t}^{2}}{p_{T0}^{2}+p_{t}^{2}}\;$ (3)

for each of the outgoing massless partons with transverse momentum $p_{t}$,
where $p_{T0}$ is a free parameter of the form factor that also enters as an
argument of the strong coupling constant $\alpha_{S}(p_{T0}^{2}+\mu_{R}^{2})$.

This suppression factor was originally proposed to remove singularity of
minijet cross sections in the collinear approach at leading-order. In the
hybrid model (or in the $k_{T}$-factorization) the leading-order cross
sections are finite as long as $k_{T}>0$, where $k_{T}$ is the transverse
momentum of the incident off-shell parton. Within this approach, a treatment
of the small-$k_{T}$ region in the construction of a given unintegrated parton
density is crucial. Different models of uPDFs may lead to different behaviour
of the cross section at small minijet transverse momenta but in any case the
cross sections should be finite. However, as it was shown in Ref.
Kotko:2016lej , the internal $k_{T}$ cannot give a minijet suppression
consistent with the minijet model and the related regularization seems to be
necessary even in this framework.

In the numerical calculations below, the intrinsic charm PDFs are taken at the
initial scale $m_{c}=1.3$ GeV, so the perturbative charm contribution is
intentionally not taken into account. In the numerical calculations below we
apply the grid of the intrinsic charm distribution from the CT14nnloIC PDF
Hou:2017khm  that corresponds to the BHPS model BHPS1980 .

###  II.3 Recombination model and charm production

Figure 3: Generic leading-order diagrams for $D$ meson production via the BJM
recombination model.

The underlying mechanism of the BJM recombination is illustrated in Fig. 3.
Differential cross section for production of $Dc$ final state reads:

$\displaystyle\frac{d\sigma}{dy_{1}dy_{2}d^{2}p_{t}}=\frac{1}{16\pi^{2}{\hat{s}}^{2}}$
$\displaystyle[x_{1}q_{1}(x_{1},\mu^{2})\,x_{2}g_{2}(x_{2},\mu^{2})\overline{|{\cal
M}_{qg\to Dc}(s,t,u)|^{2}}$ (4) $\displaystyle+$ $\displaystyle
x_{1}g_{1}(x_{1},\mu^{2})\,x_{2}q_{2}(x_{2},\mu^{2})\overline{|{\cal M}_{gq\to
Dc}(s,t,u)|^{2}}]\,.$

Above $y_{1}$ is rapidity of the $D$ meson and $y_{2}$ rapidity of the
associated $c$ or $\bar{c}$. The fragmentation of the latter will be discussed
below.

The matrix element squared in (4) can be writted as

$\overline{|{\cal M}_{qg\to Dc}(s,t,u)|^{2}}=\overline{|{\cal
M}_{qg\to({\bar{c}}q)^{n}c}|^{2}}\cdot\rho\;,$ (5)

where $n$ enumerates quantum numbers of the ${\bar{c}}q$ system
$n\equiv^{2J+1}L$. $\rho$ can be interpreted as a probability to form real
meson. For illustration as our default set we shall take $\rho$ = 0.1, but the
precise number should be adjusted to experimental data. For the discussion of
the parameter see e.g. Refs. BJM2002b ; BJM2002c  and references therein. The
asymmetries observed in photoproduction can be explained with $\rho$ = 0.15
BJM2002c .

The explicit form of the matrix element squared can be found in BJM2002a  for
pseudoscalar and vector meson production for color singlet and color octet
meson-like states.

Similar formula can be written for production of ${\bar{D}}{\bar{c}}$. Then
the quark distribution is replaced by the antiquark distribution.

In the following we include only color singlet $(q{\bar{c}})^{n}$ or
$({\bar{q}}c)^{n}$ components.

The factorization scale in the calculation is taken as:

$\mu^{2}=p_{t}^{2}+\frac{m_{t,D}^{2}+m_{t,c}^{2}}{2}\;.$ (6)

Within the recombination mechanism we include fragmentation of $c$-quarks or
$\bar{c}$-antiquarks accompanying directly produced $D$-mesons or
$\overline{D}$-antimesons, e.g.:

$\displaystyle
d\sigma[qg\to\overline{D}_{\mathrm{direct}}+D_{\mathrm{frag.}}]$
$\displaystyle=$ $\displaystyle d\sigma[qg\to\overline{D}+c]\otimes
F^{\mathrm{frag.}}_{c\to D}\;,$ (7)

where $F^{\mathrm{frag.}}_{c\to D}$ is the relevant fragmentation function.
How the convolution $\otimes$ is understood is explained in Maciula:2019iak .

In explicit calculations of $D$ meson production discussed in the next section
we shall consider proton - isoscalar nucleus scattering. For relatively light
nuclei as ${}^{20}\\!N\\!e$ the nuclear effects can be neglected. Then the
$p\\!+\\!\\!^{20}\\!N\\!e$ scattering is a superposition of $p+p$ and $p+n$
scatterings. The quark/antiquark distributions in the neutron are obtained
from those in the proton by assuming isospin symmetry.

We shall discuss in the present paper also the asymmetry in production of
$D^{0}$ meson and ${\overline{D}}^{0}$ antimeson. The asymmetry is defined as:

$A_{p}=\frac{d\sigma^{D^{0}}\\!\\!/d\xi-d\sigma^{\overline{D}^{0}}\\!\\!/d\xi}{d\sigma^{D^{0}}\\!\\!/d\xi+d\sigma^{\overline{D}^{0}}\\!\\!/d\xi}\;,$
(8)

where $\xi$ represents single variable ($y$ or $p_{t}$) or even a pair of
variables ($(y,p_{t})$).

Only a part of the pseudoscalar $D$ mesons is directly produced. A second part
originates from vector meson decays. The vector $D$ mesons promptly and
dominantly decay to pseudoscalar mesons:

$\displaystyle D^{*0}$ $\displaystyle\to$ $\displaystyle
D^{0}\pi(0.619),D^{0}\gamma(0.381),$ (9) $\displaystyle D^{*+}$
$\displaystyle\to$ $\displaystyle
D^{0}\pi^{+}(0.677),D^{+}\pi^{0}(0.307),D^{+}\gamma(0.0016),$ (10)
$\displaystyle D_{s}^{*+}$ $\displaystyle\to$ $\displaystyle
D_{s}^{+}\gamma(0.935),D_{s}^{+}\pi^{0}(0.058),D_{s}^{+}e^{+}e^{-}(6.710^{-3}),$
(11)

###  II.4 Hadronization of charm quarks

The transition of charm quarks to open charm mesons is done in the framework
of the independent parton fragmentation picture (see e.g. Refs.
Maciula:2019iak ) where the inclusive distributions of open charm meson can be
obtained through a convolution of inclusive distributions of charm
quarks/antiquarks and $c\to D$ fragmentation functions. Here we follow exactly
the method devoted to rather low energies, which was applied by us in our
previous study reported in Ref. Maciula:2021orz . According to this approach
we assume that the $D$-meson is emitted in the direction of parent
$c$-quark/antiquark, i.e. $\eta_{D}=\eta_{c}$ (the same pseudorapidities or
polar angles) and the $z$-scaling variable is defined with the light-cone
momentum i.e. $p^{+}_{c}=\frac{p^{+}_{D}}{z}$ where $p^{+}=E+p$. In numerical
calculations we take the Peterson fragmentation function Peterson:1982ak  with
$\varepsilon=0.05$, often used in the context of hadronization of heavy
flavours. Then, the hadronic cross section is normalized by the relevant charm
fragmentation fractions for a given type of $D$ meson Lisovyi:2015uqa . In the
numerical calculations below for $c\to D^{0}$ meson transition we take the
fragmentation probability $\mathrm{P}_{c\to D}=61\%$.

##  III Numerical results

Figure 4:  The rapidity (left) and transverse momentum (right) distributions
of $D^{0}$ meson (plus $\overline{D^{0}}$ antimeson) for $p+^{4}\\!\\!H\\!e$
collisions at $\sqrt{s}=86.6$ GeV together with the LHCb data Aaij:2018ogq .
Here four different contributions to charm meson production are are shown
separately, including the standard $g^{*}g^{*}\to c\bar{c}$ mechanism
(dotted), the intrinsic charm contribution (double-dot-dashed) and two
recombination components (dashed and dash-dotted). The solid histograms
correspond to the sum of all considered mechanisms. Details are specified in
the figure.

Let us now start presentation of numerical results with a discussion of the
LHCb fixed-target data on neutral open charm meson
($D^{0}+{\bar{D}}^{0}$)production in $p\\!\\!+\\!\\!^{4}\\!He$ collisions at
$\sqrt{s}=86.6$ GeV. In Fig. 4 we show the rapidity (left panel) and the
transverse momentum (right panel) distributions for the $D^{0}$-meson. Here we
plot separately four different components that lead to charm production: the
standard $g^{*}g^{*}\to c\bar{c}$ mechanism (dotted), the intrinsic charm
mechanism with the intrinsic charm probability $P_{IC}=1\%$ (double-dot-
dashed) and the two recombination contributions - the direct one (dashed) and
the one from associated fragmentation (dash-dotted). The two former
contributions are taken from our previous study of Ref. Maciula:2021orz  and
the recombination components are new. Here we assume the recombination
probability $\rho=10\%$. The solid histograms correspond to the sum of the
four contributions. The predicted distributions are compared to the LHCb data
taken from Ref. Aaij:2018ogq .

We observe that the recombination components may play an important role in
understanding the LHCb data in the region of backward rapidities. There,
according to our model the intrinsic charm contribution is the dominant one,
however, the recombination components are only slightly smaller, depending of
course on the recombination probability parameter. Both recombination
contributions are concentrated only at small meson transverse momenta and are
negligible at larger $p_{T}$’s. On contrary, the intrinsic charm component
seems to be necessary to describe the large $p_{T}$ region. Therefore, in the
present analysis we shall assume that the presence of the recombination
components should not change the $P_{IC}$ parameter in order not to weaken the
agreement between our model and the meson $p_{T}$ spectrum measured by the
LHCb. In this scenario one could use the most backward rapidity bin of the
distribution in meson rapidity to extract an upper limit for the recombination
probability parameter $\rho$. For the intrinsic charm probability $P_{IC}=1\%$
there is still a room for recombination contribution. Our study shows that the
recombination component with the probability below 5% seems to be not excluded
by the LHCb data. Within our model its inclusion does not significantly weaken
a very good description of the most backward LHCb data point where the cross
section is dominated by the intrinsic charm component. Therefore $\rho=5\%$
for the recombination mechanism might be considered here as the rough
estimation of its upper limit. The dependence of our predictions on the value
of the $\rho$ parameter is shown explicitly in Fig. 4.

Figure 5:  The same as in Fig. 4 but here different histograms correspond to a
total prediction of our model for three different values of the recombination
probability $\rho=15\%$ (dashed), $\rho=10\%$ (solid), and $\rho=5\%$ (dash-
dotted). Details are specified in the figure.

In Fig. 6 we show predictions of our model similar to those presented in Fig.
4 but here we calculate the cross sections for $p+^{20}\\!\\!N\\!e$ collisions
at $\sqrt{s}=69$ GeV within the LHCb experimental acceptance. However, the
basic conclusions are the same as in the case of the $p+^{4}\\!\\!H\\!e$
collisions at $\sqrt{s}=86.6$ GeV discussed above.

Figure 6:  The rapidity (left) and transverse momentum (right) distributions
of $D^{0}$ meson (plus $\overline{D^{0}}$ antimeson) for $p+^{20}\\!\\!N\\!e$
collisions at $\sqrt{s}=69$ GeV within the LHCb experimental acceptance. Here
four different contributions to charm meson production are shown separately,
including the standard $g*g*\to c\bar{c}$ mechanism (dotted), the intrinsic
charm contribution (double-dot-dashed) and two recombination components
(dashed and dash-dotted). The solid histograms correspond to the sum of all
considered mechanisms. Details are specified in the figure.

Finally, we discuss production asymmetry for $D^{0}$-meson and
$\overline{D}^{0}$-antimeson that can be also measured by the LHCb experiment.
The asymmetry is defined as:
$A_{p}=\frac{D^{0}-\overline{D}^{0}}{D^{0}+\overline{D}^{0}}$. According to
our present model only the recombination components may lead to the introduced
production asymmetry. Therefore, we expect different values of the asymmetry
for different values of the recombination parameter $\rho$. Future data for
the asymmetry may therefore be also used to extract the probability for the
recombination mechanism. The predictions of our model for the asymmetry
$A_{p}$ as a function of meson/antimeson rapidity (left panel) and transverse
momentum (right panel) are shown in Fig. 7. Here different histograms
correspond to different values of the $\rho$ parameter. We get smaller
production asymmetries for lower recombination probabilities $\rho$. The
asymmetry depends on the meson rapidity and transverse momentum. We predict
the production asymmetry in the most backward rapidity bin probed by the LHCb
fixed-target experiments as follows:

$\displaystyle A_{p}=-0.096\;\;\;\mathrm{for}\;\;\;\rho=5\%,$ (13)
$\displaystyle A_{p}=-0.160\;\;\;\mathrm{for}\;\;\;\rho=10\%,$ $\displaystyle
A_{p}=-0.208\;\;\;\mathrm{for}\;\;\;\rho=15\%,$

The predicted numbers can be confronted in future with the corresponding LHCb
measurement.

The asymmetry we observed in the backward rapidity region is a negative
number. It changes to a positive value when approaching midrapidities. The
sign of the asymmetry is a result of interplay between the direct and the
$c\to D$ fragmentation components of the recombination mechanism. The direct
mechanism leads to a negative asymmetry while the $c\to D$ fragmentation
component results in a positive one (see Fig. 8). The total asymmetry is the
result of the cancellation between the two components and takes a negative
sign in regions where the direct recombination component dominates over its
fragmentation counterpart.

Figure 7:  The production asymmetry $A_{p}$ for $D^{0}$-meson and
$\overline{D}^{0}$-antimeson as a function of rapidity (left) and transverse
momentum (right) for $p+^{20}\\!\\!N\\!e$ collisions at $\sqrt{s}=69$ GeV
within the LHCb experimental acceptance. Here three different values of the
recombination probability $\rho=15\%$ (dashed), $\rho=10\%$ (solid), and
$\rho=5\%$ (dash-dotted) are used. Details are specified in the figure.

Figure 8:  The production asymmetry $A_{p}$ for $D^{0}$-meson and
$\overline{D}^{0}$-antimeson as a function of rapidity (left) and transverse
momentum (right) for $p+^{20}\\!\\!N\\!e$ collisions at $\sqrt{s}=69$ GeV
within the LHCb experimental acceptance. Here the recombination probability
$\rho=10\%$ is used. The asymmetry related with the direct only (dashed) and
the $c\to D$ fragmentation only (dash-dotted) components of the recombination
mechanism are shown separately. Details are specified in the figure.

##  IV Conclusions

In the present paper we have discussed production of $D^{0}+{\bar{D}}^{0}$
mesons in fixed target $p\\!+\\!\\!^{4}\\!H\\!e$ and
$p\\!+\\!\\!^{20}\\!N\\!e$ experiments at the LHC. In addition to the
previously included gluon-gluon fusion and intrinsic charm contributions we
have included also recombination mechanism as proposed some time ago by BJM.
The latter contributions has characteristics similar to that of the intrinsic
charm.

Two contributions are included in the recombination mechanism: the directly
produced $D^{0}+{\bar{D}}^{0}$ mesons as well as $D^{0}+{\bar{D}}^{0}$
produced in fragmentation of associated $c$ or ${\bar{c}}$. The latter has
been done with the help of phenomenological fragmentation functions. In our
analysis we have used the Peterson fragmentation functions with well known
$\epsilon$ parameter.

The recombination model has one free parameter - probability of the D meson
production. This parameter can be in principle fitted to experimental data.
Our fit shows that values of the order of 5 % - 15 % seems relevant for the
fit to the cross section data for the $p\\!+\\!\\!^{4}\\!H\\!e$ scattering.

A better extraction of this parameter can be obtained from the fit to
$D^{0}-{\bar{D}}^{0}$ asymmetry. The recombination model leads to such an
asymmetry. We have discussed such an asymmetry for the
$p\\!+\\!\\!^{20}\\!N\\!e$ collisions. A future experimental data maurice
could be used to extract the $\rho$ parameter. This would better constrain
also contributions to the cross section.

Acknowledgments

We are indebted to Emilie Maurice for discussion of the LHCb
$p\\!+\\!\\!^{20}\\!N\\!e$ experimental analysis in the context of charm
production. This study was supported by the Polish National Science Center
grant UMO-2018/31/B/ST2/03537 and by the Center for Innovation and Transfer of
Natural Sciences and Engineering Knowledge in Rzeszów.

## References

  * (1) R. Maciuła, A. Szczurek and M. Łuszczak, Phys. Rev. D 92, no.5, 054006 (2015).
  * (2) R. Maciuła and A. Szczurek, Phys. Rev. D 87, no.9, 094022 (2013)
  * (3) R. Maciuła and A. Szczurek, Phys. Rev. D 100, no.5, 054001 (2019).
  * (4) A. Karpishkov, V. Saleev and A. Shipilova, [erratum: Phys. Rev. D 94, no.11, 114012 (2016)].
  * (5) R. Maciuła and A. Szczurek, Phys. Rev. D 87, no.7, 074039 (2013).
  * (6) A. van Hameren, R. Maciuła and A. Szczurek, Phys. Rev. D 89, no.9, 094019 (2014).
  * (7) R. Maciuła, Phys. Rev. D 102, no.1, 014028 (2020).
  * (8) R. Maciuła and A. Szczurek, Phys. Rev. D 105, no.1, 014001 (2022).
  * (9) S.J. Brodsky, P. Hoyer, C. Peterson and N. Sakai, Phys. Lett. B93, 451 (1980).
  * (10) R. Maciuła and A. Szczurek, J. High Energy Phys. 10, 135 (2020).
  * (11) E. Braaten, Y. Jia and T. Mehen, Phys. Rev. D66, 034003 (2002).
  * (12) E. Braaten, Y. Jia and T. Mehen, Phys. Rev. Lett. 89, 122002 (2002).
  * (13) E. Braaten, Y. Jia and T. Mehen, Phys. Rev. D66, 014003 (2002).
  * (14) T. J. Hou, S. Dulat, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, C. Schmidt, J. Winter, K. Xie and C. P. Yuan, J. High Energy Phys. 02, 059 (2018).
  * (15) V. P. Goncalves, R. Maciuła and A. Szczurek, [arXiv:2103.05503 [hep-ph]].
  * (16) R. Abbasi et al. [IceCube], Phys. Rev. D 104, 022002 (2021).
  * (17) W. Melnitchouk and A. W. Thomas, Phys. Lett. B 414, 134-139 (1997).
  * (18) E. R. Cazaroto, V. P. Goncalves, F. S. Navarra and M. Nielsen, Phys. Lett. B 724, 108-114 (2013).
  * (19) H. Holtmann, A. Szczurek and J. Speth, Nucl. Phys. A 596, 631-669 (1996).
  * (20) R. Rapp and E. V. Shuryak, Phys. Rev. D 67, 074036 (2003).
  * (21) K. P. Das and R. C. Hwa, Phys. Lett. B 68, 459 (1977) [erratum: Phys. Lett. B 73, 504 (1978)].
  * (22) A. V. Berezhnoy, V. V. Kiselev and A. K. Likhoded, Phys. Rev. D 62, 074013 (2000).
  * (23) R. Maciuła and A. Szczurek, Phys. Rev. D 97, no.7, 074001 (2018).
  * (24) R. Aaij et al. [LHCb], Phys. Lett. B 718, 902-909 (2013).
  * (25) E. Norrbin and T. Sjostrand, Eur. Phys. J. C 17, 137-161 (2000).
  * (26) S. Catani, M. Ciafaloni and F. Hautmann, Phys. Lett. B242 (1990) 97; Nucl. Phys. B366 (1991) 135; Phys. Lett. B307 (1993) 147.
J.C. Collins and R.K. Ellis, Nucl. Phys. B360, 3 (1991).
L.V. Gribov, E.M. Levin, and M.G. Ryskin, Phys. Rep. 100, 1 (1983);
E.M. Levin, M.G. Ryskin, Yu.M. Shabelsky and A.G. Shuvaev, Sov. J. Nucl. Phys.
53, 657 (1991).

  * (27) M. Ciafaloni, Nucl. Phys. B 296, 49 (1988);
S. Catani, F. Fiorani, G. Marchesini, Phys. Lett. B 234, 339 (1990);
S. Catani, F. Fiorani, G. Marchesini, Nucl. Phys. B 336, 18 (1990);
G. Marchesini, Nucl. Phys. B 445, 49 (1995).

  * (28) F. Hautmann and H. Jung, Nucl. Phys. B 883, 1-19 (2014).
  * (29) M. Deak, F. Hautmann, H. Jung and K. Kutak, JHEP 09, 121 (2009).
  * (30) A. Dumitru, A. Hayashigaki and J. Jalilian-Marian, Nucl. Phys. A 765, 464-482 (2006).
  * (31) P. Kotko, K. Kutak, C. Marquet, E. Petreska, S. Sapeta and A. van Hameren, JHEP 09, 106 (2015).
  * (32) A. van Hameren, Comput. Phys. Commun. 224, 371-380 (2018).
  * (33) T. Sjostrand and M. van Zijl, Phys. Rev. D 36, 2019 (1987).
  * (34) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen and P. Z. Skands, Comput. Phys. Commun. 191, 159-177 (2015).
  * (35) P. Kotko, A. M. Stasto and M. Strikman, Phys. Rev. D 95, no.5, 054009 (2017).
  * (36) T. J. Hou, S. Dulat, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, C. Schmidt, J. Winter, K. Xie and C. P. Yuan, JHEP 02, 059 (2018).
  * (37) R. Maciuła and A. Szczurek, J. Phys. G 47, no.3, 035001 (2020).
  * (38) C. Peterson, D. Schlatter, I. Schmitt and P. M. Zerwas, Phys. Rev. D 27, 105 (1983).
  * (39) M. Lisovyi, A. Verbytskyi and O. Zenaiev, Eur. Phys. J. C 76, no.7, 397 (2016).
  * (40) M. Cacciari, M. Greco and P. Nason, J. High Energy Phys. 05, 007 (1998).
  * (41) R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 122, no. 13, 132002 (2019).
  * (42) Emilie Maurice, a private communication.