Mathieu Dugré et al

Mathieu Dugré, Department of Computer Science and Software Engineering,
Concordia University, Quebec, Montreal, Canada,

# Performance comparison of Dask and Apache Spark on HPC systems for
Neuroimaging

Mathieu Dugré  Valérie Hayot-Sasson  Tristan Glatard  Department of Computer
Science and Software Engineering, Concordia University, Quebec, Canada
<EMAIL_ADDRESS>
###### Abstract

[Summary] The general increase in data size and data sharing motivates the
adoption of Big Data strategies in several scientific disciplines. However,
while several options are available, no particular guidelines exist for
selecting a Big Data engine. In this paper, we compare the runtime performance
of two popular Big Data engines with Python APIs, Apache Spark, and Dask, in
processing neuroimaging pipelines. Our experiments use three synthetic
neuroimaging applications to process the $606\text{\,}\mathrm{GiB}$ BigBrain
image and an actual pipeline to process data from thousands of anatomical
images. We benchmark these applications on a dedicated HPC cluster running the
Lustre file system while using varying combinations of the number of nodes,
file size, and task duration. Our results show that although there are slight
differences between Dask and Spark, the performance of the engines is
comparable for data-intensive applications. However, Spark requires more
memory than Dask, which can lead to slower runtime depending on configuration
and infrastructure. In general, the limiting factor was the data transfer
time. While both engines are suitable for neuroimaging, more efforts need to
be put to reduce the data transfer time and the memory footprint of
applications.

###### keywords:

Performance, Big Data, Dask, Spark, Neuroimaging, HPC

††articletype: WORKS19 Special Issue on Workflows in Support of Large-Scale
Science GEs: David Walker, Chase Wu

This is the accepted version of the following article: Dugré, M., Hayot-
Sasson, V., & Glatard, T. (2023). Performance comparison of Dask and Apache
Spark on HPC systems for neuroimaging. Concurrency and Computation: Practice &
Experience, 35(21). https://doi.org/10.1002/cpe.7635, which has been published
in final form at [https://onlinelibrary.wiley.com/doi/10.1002/cpe.7635]. This
article may be used for noncommercial purposes in accordance with the Wiley
Self-Archiving Policy [http://www.wileyauthors.com/self].

##  1 Introduction

The rise in data sharing coupled with improved data collection technologies
leads neuroimaging into the Big Data era1, 2, 3. While standard neuroimaging
workflow engines, such as Nipype4, are well rounded to process standard
compute-intensive pipelines, they lack Big Data strategies (i.e., in-memory
computing, data locality, and lazy-evaluation) to improve the performance of
increasingly prevalent data-intensive pipelines. A previous study5 noted the
importance of in-memory computing and data locality to improve the performance
of data-intensive pipelines. This work studies how performance benchmarks of
neuroimaging applications generalize across Big Data engines in an HPC
environment. With the rising need for Big Data engines in neuroimaging and
none commonly used in the field, evaluating the performance of different
frameworks is essential. In particular, we study in-depth the performance
difference between Dask6 and Apache Spark7 for their suitability to process
neuroimaging pipelines. Our main contribution is to assess whether Spark or
Dask has any substantial performance advantage in processing neuroimaging
pipelines in HPC systems.

Spark and Dask both provide in-memory computing, data locality, and lazy-
evaluation, typical for Big Data engines. The scheduler for both engines
operates dynamically, which benefits applications with runtime unknown ahead
of time6. They also provide rich high-level APIs and support a variety of
schedulers and deployment environments, such as Mesos8, YARN9, Kubernetes, and
HPC clusters. Although sharing similarities, these engines have differences.

First and foremost, Spark is written in Scala while Dask is in Python. The
popularity of Python in the scientific community arguably gives an edge to
Dask due to the serialization cost between Python and Scala. However, if not
careful, the Python Global Interpreter Lock (GIL) can significantly reduce the
parallelism of an application in some cases. The difference in programming
languages also provides different benefits. On the one hand, as part of the
Scipy ecosystem, Dask provides almost transparent integration with APIs such
as Numpy arrays, Pandas dataframe, or RAPIDS GPU accelerated code framework.
On the other hand, Spark’s Java, Scala, R, and Python APIs efficiently
parallelize pipelines with minimal performance loss due to the Java Virtual
Machine. Our study is restrained to performance, although we recognize that
other factors affect the choice of an engine in practice.

While laptops or workstations can be sufficient for some applications,
neuroimaging pipelines often require larger infrastructures to process
datasets in the scale of several hundreds of $\text{\,}\mathrm{GiB}$ to few
$\text{\,}\mathrm{TiB}$. This paper focuses on the performance of Big Data
engines in HPCs environments for neuroimaging. We installed the Slurm
scheduler and Lustre file system 10 on a dedicated cluster to mimic a typical
HPC system. Lustre is a reference file system for HPC environments, used in
over 60 TOP100 supercomputers 11, offering parallel file accesses.

In a previous study12 we found no significant performance difference between
Apache Spark and Dask for data-intensive applications. However, using the
network file system (NFS) led to an I/O bottleneck limiting application
performance. Another study 13 also reported network and I/O bottlenecks when
concurrency increased. This same study found that the startup overhead time
for Spark was more prominent than for Dask. On the contrary, the work in 14
claims that the startup overhead is more significant for Dask than for Spark.
This contradiction suggests that the startup overhead might vary depending on
cluster configuration. The latter study processed a $100\text{\,}\mathrm{GB}$
dataset using a neuroimaging pipeline and found Dask to be up to 14% faster
than Spark due to “more efficient pipelining" and serialization time to
Python. On the contrary, the work in 15 shows that with large datasets, Spark
provides better speedup than Dask. However, that study used the Dask Bag API,
which has performance limitations compared to other Dask APIs 16.

The above studies compared the performance of the engine at a high level. In
comparison, our study provides a detailed analysis of the performance
differences and similarities between Dask and Spark for data-intensive
applications executed on HPC systems.

The following section introduces the Lustre file system and the Big Data
engines. Next, we present the design of our benchmarking experiments. We
consider two types of neuroimaging pipelines: high-resolution imaging and
functional MRI studies for a cohort of subjects. The pipelines we chose
involve different compute and I/O patterns, e.g., map-reduce or map-only,
requiring partial data in memory or the whole dataset. Further sections
present our results, discussion, and conclusion.

##  2 Background

###  2.1 Lustre

The Lustre file system is a parallel file system of reference in HPC
environments. Lustre offers a POSIX-compliant interface with the file system.
Installed in some of the most powerful supercomputers, the design of Lustre
allows scaling to thousands of nodes, petabytes of storage, and hundreds of
gigabits per second of I/O throughput.

At the core of Lustre are the management servers (MGSs) and metadata server
(MDSs) in charge of the file operations and metadata, the object storage
targets (OSTs) store file data, and the object storage servers (OSSs) handle
the OSTs. Lustre clients make requests to the MGSs to access files on the
OSTs. To handle failures in the system, the MGSs, MDSs, and OSTs have failover
capabilities. InfiniBand, Ethernet, or a combination of both can be used to
interconnect these different components.

While other filesystems such as HDFS or cloud storage could be considered, we
chose Lustre to represent an HPC environment best. HDFS uses permanent local
storage on compute node, which is usually not feasible in HPC. Moreover, HPC
nodes rarely have access to the internet, limiting cloud storage use.

###  2.2 Dask

Dask is a Python-based Big Data engine with growing popularity in the
scientific Python ecosystem. Dask was designed with data locality and in-
memory computing in mind to mitigate the data transfer bottleneck in Big Data
workflows. Data locality, popularized by Map-Reduce 17, schedules tasks where
the data reside. In-memory computing minimizes the overhead of transferring
data to disks by keeping data in memory. Dask uses lazy evaluation to reduce
unnecessary communication and computation. The engine builds a dynamic graph
before execution, allowing it to determine which task to compute. Dask
workflows can further reduce data transfers by leveraging multithreading
whenever the Python GIL does not restrict it. It achieves fault tolerance by
recording data lineage, the sequence of operations used to modify the initial
data.

Dask offers five data structures: Array, Bag, DataFrame, Delayed, and Futures.
Arrays offer a clone of NumPy API for distributed processing of large arrays.
Bags are a distributed collection of Python objects that offers a programming
abstraction similar to PyToolz. Dataframes are a parallel composition of
Pandas Dataframes used to process a large amount of tabular data. Dask Delayed
offers an API for distributing arbitrary functions that do not fit the above
frameworks. Lastly, Dask Futures can also execute arbitrary functions;
however, it launches computation immediately rather than lazily. Dask allows
users to install only required components making it lightweight.

In Dask, a scheduler decides where and when to execute tasks using the Dask
graph. API operations generate multiple fine-coarse tasks in the computation
graph, allowing a more straightforward representation of complex algorithms.

The Dask engine is compatible with multiple distributed schedulers, including
YARN and Mesos. It also provides its scheduler, Dask Distributed scheduler. We
chose to use Dask Distributed scheduler to keep the environments similar
between Dask and Spark.

In the Dask Distributed scheduler, a dask-scheduler process administrates the
resource provided by dask-workers in the cluster. The scheduler receives jobs
from clients and assigns tasks to available workers using a global FIFO
(First-In-First-Out) job scheduling policy. However, the Dask scheduler
attempts to complete immediate task dependencies using a LIFO (Last-in-First-
Out) policy to reduce the memory footprint of the workers.

Dask offers multiple ways to deploy a cluster, including, but not limited to,
SSH configs, Kubernetes, SLURM, PBS. For our experiments, we used the Dask
SLURM cluster API.

###  2.3 Apache Spark

Apache Spark is a widely-used general-purpose Big Data engine. Like Dask, it
aims at reducing data transfer costs by incorporating data locality, in-memory
computing, and lazy evaluation.

Spark offers three options to schedule jobs: Spark Standalone, Mesos, and
YARN. Spark Standalone is a simple built-in scheduler. YARN is mainly used to
schedule Hadoop-based workflows, while Mesos can be used for various
workflows. We limit our focus to Spark Standalone scheduler, as researchers
are likely to execute their workflows in an HPC environment where neither YARN
nor Mesos is usually available.

In the Spark Standalone scheduler, a leader (a.k.a master) coordinates the
resource provisioned by workers in the cluster. A driver process receives jobs
from clients and requests workers from the leader. Jobs are divided into
stages to execute onto workers. Each operation in a stage is represented by a
high-level task in the computation graph. Like Dask, Spark Standalone
scheduler uses a LIFO policy to schedule tasks. Spark Standalone has two
execution modes: (1) the client mode, where the driver process runs in a
dedicated process, and (2) the cluster mode, where the driver runs within a
worker process. Our experiments use the client mode since cluster mode is not
available in PySpark.

Spark’s primary data structure is the Resilient Distributed Dataset (RDD)18, a
fault-tolerant, parallel collection of data elements. RDDs are the basis of
other Spark data structures: Datasets and DataFrames. Datasets are similar to
RDDs but benefit additional performance by leveraging Spark SQL’s optimized
execution engine. DataFrames are Datasets organized into named columns and are
used to process tabular data. While the DataFrame API is available in all
supported languages, Datasets are limited to Scala and Java.

Python is a standard programming language in the scientific community,
offering numerous data processing libraries. While serialization from Python
to Java, an operation required when using Spark’s Python API, creates
overhead, we found it minimal in a previous study 12. We focus on PySpark API
to have a similar environment between the different engines and its
suitability for neuroimaging.

##  3 Methods

###  3.1 Infrastructure

We used the “slashbin" cluster at Concordia University. The cluster has eight
compute nodes, four storage nodes, one login node, and one Lustre metadata
node. Each compute node has two 16-core Intel(R) Xeon(R) Gold 6130 CPU @
2.10GHz, $250\text{\,}\mathrm{GiB}$ of RAM, $126\text{\,}\mathrm{GiB}$ of
tmpfs, CentOS 8.1 and Linux kernel 4.18.0-240.1.1.el8_lustre.x86_64. The
compute and storage nodes are inter-connected by
$10\text{\,}\mathrm{Gbit}\text{\,}{\mathrm{s}}^{-1}$ Ethernet link through a
Mellanox MSN2410 switch. The login node, hosting an NFS server, is connected
to the compute nodes through the same switch, with
$10\text{\,}\mathrm{Gbit}\text{\,}{\mathrm{s}}^{-1}$ Ethernet.

We configured Lustre to have 1 MGS and 4 OSS. Each OSS has 11 OSTs (disks) of
$8.8\text{\,}\mathrm{TiB}$. The OSTs are configured to have a maximum of
$1000\text{\,}\mathrm{MiB}$ of dirty data written into the client page cache.
Both the MDT and OSTs are configured to a maximum of 64 concurrent RPCs in
flight. A timeout of $100\text{\,}\mathrm{s}$ is set for Lustre, and OSTs
checksum are disabled. Lustre 2.14.0-1.el8 is installed on both the Lustre
clients and servers. We used the default inode configuration for Lustre.

While our infrastructure is not of the scale of TOP500 HPC systems, we can
assure that no external processes impacted our measurements which would be
hard to achieve on a production system. Our system uses similar hardware as
Compute Canada’s Béluga cluster — a representative cluster available to
researchers in Canada — and replicates as closely as possible its
configuration. This lets us compare Spark and Dask in an isolated environement
with the condition of a realistic allocation of resources, for neuroimaging
applications, in an HPC system.

We configure Dask to have eight worker processes per node, each with eight
threads. This is the default recommended configuration for the Dask SLURM
cluster configuration 19. Similarly, we configured Spark to have eight
executors per node (hereafter called workers for consistency with Dask), each
with eight threads. However, when trying different configurations for the
workers, we did not observe performance differences for our benchmark.
Therefore we chose the recommended configuration from Dask. Each worker was
allocated $31.5\text{\,}\mathrm{GiB}$ of memory. For Spark, the JVM heap space
uses 10% of workers’ memory, and the location of the log is set on an NFS-
mounted directory. The worker workspace for both engines is located on Lustre
due to the limited fast storage space on the worker nodes. Other Dask
configurations are left to default. A new Dask or Spark cluster was started
and torn down for each experiment. Dask 2021.7.0 and Spark 3.1.2 were used.

###  3.2 Dataset

We used BigBrain20, a 3-D image of the human brain with voxel intensities
ranging from 0 to 65,535. We converted the blocks into the NIfTI format, a
popular format in neuroimaging. We left the NIfTI blocks uncompressed,
resulting in a total data size of $606\text{\,}\mathrm{GiB}$. To evaluate the
effect of block size, we split the image into 1000, 2500, and 5000 files of
$606\text{\,}\mathrm{MiB}$, $242.4\text{\,}\mathrm{MiB}$, and
$121.2\text{\,}\mathrm{MiB}$, respectively. We used the sam 21 library to
split the image.

We also used the dataset provided by the Consortium for Reliability and
Reproducibility (CoRR) 22, freely available on DataLad 23. The entire dataset
is $379.83\text{\,}\mathrm{GiB}$, containing anatomical, diffusion, and
functional images of 1,654 subjects acquired in 35 sites. We used all 3,491
anatomical images, representing $26.67\text{\,}\mathrm{GiB}$ overall
($7.82\text{\,}\mathrm{MiB}$ per image on average).

We chose these dataset as they represent common dataset sizes for neuroimaging
imaging studies and are open-access.

###  3.3 Applications

This section describes the applications used to benchmark Dask and Spark. We
use three simple synthetic applications to have a deep understanding of the
underlying behavior of the applications with simple I/O patterns, namely
Increment, Multi-Increment, and Histogram. Also, two neuroimaging applications
are used to study more realistic applications, namely Kmeans and BIDS App
Example24. We chose these neuroimaging data processing applications to study
the effects of various Big Data characteristics. Table 1 depicts the Big Data
characteristics used in each application. In this table, Map and Reduce
represent the Map-Reduce paradigm operations, and the Shuffle column
represents the amount of data shuffled between stages if any. The column
Dataset in-memory depicts whether the whole dataset is required to be in
memory for processing. Then Compute intensive and Data intensive shows if the
applications are compute or data intensive. Lastly, the Containerized column
shows applications wrapped in a container.

Table 1: Big Data characteristics of the applications | Map | Reduce |  | Data
---
shuffling
| Dataset
---
in-memory
| Compute
---
intensive
| Data
---
intensive
Containerized
Increment | X |  |  | Partial |  | X |
Histogram | X | X | Minimal | Partial |  | X |
Kmeans | X |  | X | Whole | X |  |
BIDS App Example | X | X | Minimal | Partial | X |  | X

To minimize the performance difference from implementation decisions, the code
base of the applications is common for both engines. Therefore, the
performance observed comes directly from the engine scheduling and their APIs.
The scripts used for our experiments and the results generated are available
on GitHub at: https://github.com/mathdugre/paper-big-data-engines/.

####  3.3.1 Increment

We adapted the increment application used in 5. This synthetic application
reads blocks of the BigBrain from Lustre and simulates computation by sleeping
for a specified period (Fig. 2). To simulate intermediate results, we repeat a
$1\text{\,}\mathrm{s}$ sleep delay during each incrementation. We prevent data
caching of the blocks by incrementing their voxels value by one after each
sleep operation. Finally, we write the resulting NIfTI image back to Lustre.
This application allows us to study the engines when their inputs are
processed independently. The map-only pattern mimics the processing of
multiple independent subjects in parallel.

####  3.3.2 Histogram

As our second application, we calculate the histogram of the BigBrain image.
The application reads the BigBrain blocks from Lustre, calculates each
intensity’s frequency, and then writes the aggregated result back on Lustre
(Fig. 2). This map-reduce application has a very high read over write ratio.
Moreover, this application requires shuffling, albeit of a limited amount of
data. The amount of inter-worker communication is in-between the increment and
multi-increment applications.

Figure 1: Task graph for Incrementation with 3 iterations and 3 BigBrain
blocks.

Figure 2: Task graph for Histogram with 3 BigBrain blocks.

####  3.3.3 Kmeans

For our third application, we apply Kmeans clustering to the voxel intensities
of the BigBrain image. We set the number of clusters to 3. The application
starts by reading the image blocks, combining all voxels in a 1-D array, and
choosing initial centroids using the min, max, and intermediate values. It
assigns each voxel to its closest centroid and updates each centroid with the
average of its assigned voxels. It repeats the assignment and updates steps
for a configurable number of iterations. Finally, the voxels of the image
blocks are classified and written back to the file system. Updating the
centroids involves substantial data communication between the workers. The
computation to find the closest centroids for each voxels scales with the
number of centroids. It can be computed in parallel at the cost of higher
memory or sequentially to reduce memory usage. We decided to compute the
nearest centroids in parallel since we use only three centroids. However, due
to resources limitation on our infrastructure, the Kmeans experiments are
performed with only half of the BigBrain dataset ($303\text{\,}\mathrm{GiB}$).
We chose this application as it involves substantial data shuffling between
workers, and requires to keep all the data in memory.

####  3.3.4 BIDS App Example

Our fourth application is the BIDS App Example, a neuroimaging pipeline to
measure brain volume from MRIs. For this application, we use the CoRR dataset.
The application extracts the brain volume of each participant then computes
the average for each group of participants. Unlike the other applications, the
BIDS App Example is a command-line executed in a Docker image (bids/example on
DockerHub). We converted the Docker image to a Singularity image for use in
HPC environments, using docker2singularity We chose this application for its
real-life use case and the fact that it is containerized unlike our other
applications.

###  3.4 Experiments

Table 2 shows the parameters that we varied throughout the experiments. We
varied (1) the number of nodes to assess the scalability of the engines’
scheduler, (2) the number of files to segment BigBrain to measure the effect
of different I/O patterns and parallelization degrees, and (3) the number of
iterations to evaluate the effect of the number of tasks. It should be noted
that increasing the number of iterations for a given sleep delay also
increases the total compute time. The choice of parameters for our experiments
was restrained by the datasets and resources available from the cluster. The
unrestrained parameters were chosen to have an expected linear scaling to
facilitate the interpretation of the results.

To avoid potential external biases due to caching, background processes and,
network load, we ran the applications in randomized order and cleared the page
cache between every experiment. Moreover, we ran each benchmark ten times.

For each execution, we measured the time spent in the different functions of
the application to read, process, and write data. We calculate the overhead
for each CPU thread by subtracting the total time spent in the reading,
writing, and processing functions to the makespan of the application. Summing
those results gives the total overhead for the application.

Table 2: Parameters for the experiments | Increment | Histogram | Kmeans | BIDS App Example
---|---|---|---|---
Dataset size [$\text{\,}\mathrm{GiB}$] | 606 | 303 | 26.67
# of Nodes | 2, 4, 8
Number of file | 1000, 2500, 5000 | 5000 | 3491
# of Iterations | 1, 5, 25 | - | 1, 3 | -

##  4 Results

###  4.1 Increment

Figure 3 depicts the total execution time spent in each code segment of the
Increment application. The bars are broken down into Overhead, Read, Write,
and other application-specific functions. The mean of each bar is calculated
from ten repetitions, with the error bar representing the standard deviation.
As expected, the compute time, from Increment, remains the same when we vary
the number of nodes. However, the Read and Write time increases with the
number of nodes for both engines. Due to the large amount of worker threads
accessing the Lustre file system concurrently, the disk bandiwdth or network
are saturated, which results in an I/O bottleneck.

From Figure 3, we also observe that Dask’s overhead is higher than Spark’s.
While the overhead of both engines increases with the number of nodes, Dask’s
overhead is significantly worse when scaling. Although Dask has a higher
overhead, the total time for both engines is overall similar. This is due to a
tradeoff between disk bandwidth and compute/overhead time. On the one hand,
when multiple threads access Lustre concurrently they may saturate the shared
disk bandwidth, resulting in longer I/O time. On the other hand, when most
threads are either idle or performing computations, the contention on the disk
bandwidth is reduced, resulting in shorter I/O time. The only significant
difference in total time is when two nodes are used. This results in an
approximately 3% total runtime difference.

Figure 3: Increment: varying nodes – 5000 files, 5 iterations

Figure 4 shows the detailed execution time for the Increment application while
varying the number of iterations. As expected, the overhead for both engines
increases with the number of iterations. However, overhead for Dask increases
significantly faster than Spark, resulting in a 7% makespan difference for 25
iterations.

Figure 4 also shows that the total I/O time decreases as total computation
time increases. This is explained by an I/O bottleneck on the Lustre file
system. When threads are busy with computations, it reduces the contention on
Lustre, thus improving I/O speed for individual accesses.

Figure 4: Increment: varying iterations – 4 nodes, 5000 files

Figure 5 shows an increase in Increment time. This is because the application
performs an incrementation for each file, thus increasing the computation
time. Similarly to Figure 4, increasing the total compute time of the
application reduces the total I/O time. Again, this is explained by the I/O
tradeoff between compute time and contention of disk bandwidth.

From Figure 5 we also observe that Dask has significantly higher overhead than
Spark. However, the total execution time of the application is equivalent due
to the tradeoff between I/O and compute/overhead time.

Figure 5: Increment: varying # of files – 4 nodes, 5 iterations

###  4.2 Histogram

Figure 6 depicts the execution time for the Histogram application when varying
the number of nodes. Overall, the execution time for Spark is significantly
lower than for Dask by up to 5%, with the difference coming from a shorter
execution time for Calculate_histogram. This is unexpected as the code base
for both engines is the same, with only the calling APIs differing. It is most
likely due to the use of Dask Bags, which are slower than the other Dask APIs
16. As with previous experiments, the overhead for Spark is less and scales
better, although compensated by a longer total I/O time.

Figure 6: Histogram: varying nodes – 5000 files

Figure 7 shows the execution time for the Histogram when varying the number of
files. The results are similar to varying the number of nodes with Spark
having faster compute time for Calculate_histogram. This results in a
significant makespan difference with Spark being up to 3% faster.

Figure 7: Histogram: varying # of files – 4 nodes

###  4.3 Kmeans

When varying the number of workers, we observed that Dask is multifold faster
than Spark, resulting in a $1.79\times$ to $2.71\times$ speed-up ratio on the
makespan. The problem we observed was twofold. Firstly, the process ID of the
workers’ threads is not reported correctly and leading to falsely excessive
overhead for Spark. Secondly, the Spark workers were restarting due to a lack
of memory, leading to task recomputation and having to reload data into
memory.

We increase the worker memory by reducing the number of worker processes per
compute node since all available memory from the compute node is already
allocated to the pool of workers. When reducing the number of processes from
64 to 32 for Spark and Dask (effectively doubling the memory of each process),
the difference in execution time between Spark and Dask was reduced
significantly. With only 32 processes and more memory per worker, Spark no
longer recomputed tasks for the Kmeans application. We suppose that given more
memory per thread/process, Dask and Spark would then perform similarly.

In a follow-up experiment, we further optimized our Kmeans implementations for
both Dask and Spark to minimize the size of data kept in memory. Figure 8 and
Figure 9 show the Gantt chart for the execution of this new experiments for
Dask and Spark. The x-axis shows the time in second while the y-axis depicts
the workers, with one worker thread ID (Dask) and process ID (Spark) per line.
We represent each task on the Gantt chart with its associated function’s
color. We again observe that Spark reports three times more processes than
Dask threads. It is unclear what is the cause for reporting extra Spark
processes. However, we found that both Dask and Spark use the same amount of
maximum concurrent processes or threads. Additionally, there was no
significant difference in the average amount of concurrent processes or
threads between Dask and Spark.

Figure 8 and Figure 9 show that Dask has a significantly shorther makespan
than Spark, $338\text{\,}\mathrm{s}$ and $409\text{\,}\mathrm{s}$
respectively. The main differences in makespan for Spark appear from a longer
initial load of the data and a longer scheduling delay between the data load
and the initial computation.

Figure 8: Kmeans: Dask with 8 nodes, 5000 files, and 3 iterations

Figure 9: Kmeans: Spark with 8 nodes, 5000 files, and 3 iterations

For brevity we omitted the figures when varying the number of iterations as
the results are similar to varying the number of workers.

###  4.4 BIDS App Example

Figure 10 shows the total execution time for the BIDS App Example application
when varying the number of nodes. We observe that Dask has a slightly lower
overhead, thank Spark, although not significant. However, the makespan for
Dask is between 7.5% to 14.5% lower than for Spark. The main difference comes
from an increased amount of stagger tasks for Spark, thus increasing the
overhead of the application. Moreover, although minimal, a shorter cluster
deployment time for Dask reinforces the difference in overhead time.

Figure 10: BIDS App Example: varying nodes

##  5 Discussion

###  5.1 I/O bottleneck

We used Lustre as shared file system as it is a de facto solution in HPC
systems. The performance of our Lustre installation degrades when a large
number of concurrent tasks attempt to perform I/O to the shared file system,
as shown in Figures 3 and 6. In a previous study 12, we observed nearly no
improvement when increasing the number of nodes to process an application
while using an NFS. In this paper, the performance scales when increasing the
number of nodes, although not linearly due to an I/O bottleneck. The system
was designed such that the filesystem has a similar speed to the network.
Upgrading to more performant networking and storage components would reduce
the bottleneck when increasing the number of workers.

As in our previous study 12, we observe that the I/O bottleneck has a lower
impact on more compute-intensive applications, such as Kmeans and BIDS App
Example. However, it significantly affects data-intensive applications, such
as Increment and Histogram.

These results suggest that increasing the parallelism of the file system and
speed of the network is critical when scaling the compute resources to process
a data-intensive application. Otherwise, the compute resource might be
underutilized.

###  5.2 Tradeoff between overhead and data transfer

Overall, except for Kmeans, the overhead for Dask is larger than for Spark.
However, the total execution time for both engines is comparable. This almost
exact compensation between the engine overhead and the data transfer time (see
for example Figure 3) is not surprising. When there is high contention on the
bandwidth of the file system, and in the absence of trashing, desynchronizing
the file transfers generally does not reduce the total data transfer time. To
take an extreme example, the makespan of $n$ concurrent file transfer of size
$F$ from a file system with a bandwidth $B$ would be $nF/B$, as the $n$
transfer would share the bandwidth $B$, which is equivalent to sequentially
transferring the $n$ files. In a previous study 12, we also observed this
tradeoff between overhead and data transfer when the file system was
saturated. Exploring strategies that take advantage of this tradeoff could be
interesting to reduce the application makespan.

###  5.3 Memory management

Our initial Kmeans results showed a significant difference between Dask and
Spark. The difference in performance came from Spark worker restarting due to
running out of memory; thus, Spark needed to recompute several tasks. This is
unexpected, given that Dask and Spark were configured similarly and were
allocated the same amount of resources. We think those issues arose from a
combination of (1) the reserved memory for the JVM heap space, (2) the
serialization of Python to Java kept in memory, (3) Spark using processes
compared to Dask using threads, thus limiting shared memory consumption, and
(4) different policies to persist data to disk when running out of memory.

We optimized the memory usage in our Kmeans implementation for Dask and Spark
and observed that the Spark workers no longer failed with this improved
implementation. Nevertheless, Figure 8 and Figure 9 depicts that a substantial
difference in makespan remained between Spark and Dask with this
implementation, showing that the previously-based gap was not only due to
worker failures.

###  5.4 Scheduling overhead

Overall, our results show that Dask has a slightly longer startup overhead
than Spark. However, we found that this did not lead directly to the makespan
difference. While both engines showed similar overhead for applications with
low inter-worker communication, Spark had significantly more scheduling
overhead than Dask for the Kmeans applications. Spark was considerably slower
to schedule the first set of computations after reading the data on the
workers. Moreover, Spark read all the data before starting computation while
Dask was still reading data alongside doing computation. This difference in
scheduling is beneficial for data-intensive applications as it reduces the
risk of I/O bottleneck, as explained in Section 5.2.

###  5.5 Comparison to other studies

Overall, except for Kmeans, our experiments did not show substantial
performance differences between Spark and Dask. In contradiction with the
experiments of 14, we found that Dask had a larger overhead than Spark,
especially when scaling the number of workers. However, this was compensated
by the previously mentioned overhead and data transfer tradeoff. Since our
focus is on neuroimaging applications, we conclude that the performance of
Dask and Spark is equivalent.

In this study, we observed that Dask has a larger startup overhead than Spark.
These results are consistent with a previous study 14. However, other studies
found the opposite 13 or that the difference was negligible 12. Those apparent
contractions most likely come from differences in infrastructure and
configuration of the engines.

Our results from the Histogram experiments show that the Dask Bags API was
slower than Spark. The experiments from 15 also found that Dask Bags performed
slower than Spark. We conclude that when using Dask, the other Dask APIs
should be favored unless Dask Bags are strictly required.

Our results showed that, for data-intensive applications, increasing the
number of workers should only be done if the bandwidth of the file system is
not saturated. Two other studies support this 12, 13. Therefore, it is
critical to ensure enough bandwidth to the file system before scaling the
compute resources for data-intensive applications. Otherwise, those extra
resources might be underused.

Another study 25 developed a parametrized benchmark framework to evaluate the
performance of different engines, including Spark and Dask. On the one hand,
the framework reduces the implementation effort and helps with comparing
frameworks fairly. On the other hand, it is hard to characterize neuroimaging
realistically applications using their framework since the runtime of these
applications can vary greatly depending on the subject process, and the
framework lacks kernels with a variable amount of compute-intensity. Due to
these limitations, we decided to implement the benchmarked applications
directly into Spark and Dask.

###  5.6 Recommendations

The results from our experiments show that the performance difference between
Dask and Spark are minimal, in general. The main differences we observed are
qualitative. We found that Dask integration with the scientific Python
ecosystem made it easier to implement neuroimaging applications. Spark
configuration required to set multiple environment variables. In comparison,
Dask is ready-to-go after an installation through PIP, the Python package
manager. The Dask dashboard provides more information than Spark’s, making the
debugging of Dask applications easier and faster. Therefore, for the
neuroimaging field, we recommend the use of Dask over Spark. However, in other
fields with non-Python-based ecosystem, both Dask and Spark may be equal
candidates as a Big Data engine.

While the engines showed similar performance, they both come with their
caveats that can drastically affect their performance. For Dask, the user
should be aware of the Python GIL limitations. For Python code, multiples
worker processes should be used while for code releasing the GIL, a balance
between processes and thread is recommended. For Spark, we experienced severe
performance issues due to workers running out of memory. Users should
carefully consider their Spark worker configurations to prevent similar
issues.

###  5.7 Limitations

To make a fair comparison between Spark and Dask, we set both engine
configurations to be equivalent to the best of our capabilities, which comes
with some limitations. Firstly, both engines use a shared implementation for
the algorithms. This is done to keep the theoretical runtime equivalent;
however, the chosen implementation might not be optimal for either of the
engines. Secondly, we decided to use Lustre as a file system for both Spark
and Dask. Spark could benefit from using HDFS as a file system, but HDFS is
uncommon in HPC. Finally, our Lustre filesystem suffers from an I/O bottleneck
coming from either the network or disk bandwidth. As explained in section 5.2,
the data transfer time is tightly correlated to the overhead of the
application. However, while reducing this bottleneck should improve the
performance of both engines similarly, it can be hard to predict the actual
behavior.

While both Dask and Spark support various high-performance network strategies
such as UCX, we focused on the standpoint of neuroscientists that would
commonly use the frameworks mostly out of the box. With higher network
throughput, the overall performance of the engines could differ. Further work
would be required to compare Dask and Spark using different high-performance
network protocols.

Dask and Spark both support GPUs; however, we did not use applications that
leverage them in our work. Although several neuroimaging applications use
GPUs, only a few preprocessing pipelines use them. Therefore, we omitted GPUs
from this study because preprocessing comprises the majority of computing time
for neuroimaging workflow. Moreover, the few preprocessing applications using
GPUs are not commonly used due to their recent development.

Due to the limited amount of benchmarks available to evaluate the performance
of neuroimaging applications, we use our own collection of applications.
Further efforts would be needed to create a credible and in-depth performance
benchmark specific to the neuroimaging field.

##  6 Conclusion

We presented a detailed comparison of two Big Data engines, Apache Spark and
Dask. We studied the engines using three synthetic neuroimaging applications
and a neuroimaging pipeline to calculate brain volume. Overall, our results
showed no substantial difference between the engines. Although Spark overhead
was generally lower than Dask, it was significantly higher when applications
incur inter-worker communication. For one application, our results showed that
Spark used more memory, causing issues with its workers resulting in a
significantly faster runtime for Dask. However, using an infrastructure
offering more memory per worker would most likely solve this issue.

These results suggest that future research should focus on strategies to:

  * •

Limit the impact of data transfers for data-intensive applications

  * •

Reduce the footprint of application memory

  * •

Manage the memory of workers more efficiently

## Acknowledgment

Mathieu Dugré was partially funded by a Alexander Graham Bell scholarship from
the Natural Sciences and Engineering Research Council of Canada and a Master’s
(B1X) Research Scholarships. The Canada Research Chairs program partially
funded this work. The Canada Foundation for Innovation funded the HPC cluster.

## Data Availability

### BigBrain

The data that support the findings of this study are available from the
corresponding author, M. Dugré, upon reasonable request. A lower resolution of
$40\text{\,}\mathrm{\SIUnitSymbolMicro m}$ is available publicly by FTP at:
https://ftp.bigbrainproject.org/bigbrain-
ftp/BigBrainRelease.2015/3D_Blocks/40um/nii

### CoRR

The data that support the findings of this study are openly available in
DataLad at:
https://datasets.datalad.org/?dir=/corr/RawDataBIDS

### Code

The code that support the findings of this study is openly available in GitHub
at:
https://github.com/big-data-lab-team/paper-big-data-engines/tree/CPE

## Conflict of Interests

The authors declare no conflict of interests.

## References

  * 1 Alfaro-Almagro F, Jenkinson M, Bangerter NK, et al. Image processing and Quality Control for the first 10,000 brain imaging datasets from UK Biobank. NeuroImage 2018; 166: 400-424.  doi: https://doi.org/10.1016/j.neuroimage.2017.10.034
  * 2 Van Horn JD, Toga AW. Human neuroimaging as a "Big Data" science. Brain imaging and behavior 2014; 8(2): 323–331.
  * 3 CONP . Canadian Open Neuroscience Platform (CONP) Portal. https://portal.conp.ca; .
  * 4 Gorgolewski K, Burns C, Madison C, et al. Nipype: A Flexible, Lightweight and Extensible Neuroimaging Data Processing Framework in Python. Frontiers in Neuroinformatics 2011; 5: 13.  doi: 10.3389/fninf.2011.00013
  * 5 Hayot-Sasson V, Brown ST, Glatard T. Performance Evaluation of Big Data Processing Strategies for Neuroimaging. In: 2019 19th IEEE/ACM International Symposium on Cluster, Cloud and Grid Computing (CCGRID). ; 2019: 449-458
  * 6 Rocklin M. Dask: Parallel Computation with Blocked algorithms and Task Scheduling. In: Proceedings of the 14th Python in Science Conference. ; 2015: 126 - 132
  * 7 Zaharia M, Xin RS, Wendell P, et al. Apache Spark: A Unified Engine for Big Data Processing. Commun. ACM 2016; 59(11): 56–65.  doi: 10.1145/2934664
  * 8 Hindman B, Konwinski A, Zaharia M, et al. Mesos: A Platform for $\\{$Fine-Grained$\\}$ Resource Sharing in the Data Center. In: 8th USENIX Symposium on Networked Systems Design and Implementation (NSDI 11). ; 2011.
  * 9 Vavilapalli VK, Murthy AC, Douglas C, et al. Apache Hadoop YARN: Yet another resource negotiator. In: ACM. ; 2013: 5.
  * 10 Braam PJ. The Lustre Storage Architecture. ArXiv 2019; abs/1903.01955.
  * 11 OpenSFS . Lustre® File System, Version 2.4 Released. https://www.opensfs.org/press-releases/lustre-file-system-version-2-4-released/; .
  * 12 Dugré M, Hayot-Sasson V, Glatard T. A Performance Comparison of Dask and Apache Spark for Data-Intensive Neuroimaging Pipelines. In: 2019 IEEE/ACM Workflows in Support of Large-Scale Science (WORKS). ; 2019: 40-49
  * 13 Chantzialexiou G, Luckow A, Jha S. Pilot-Streaming: A Stream Processing Framework for High-Performance Computing. In: 2018 IEEE 14th International Conference on e-Science (e-Science). ; 2018: 177-188
  * 14 Mehta P, AlSayyad Y, Dorkenwald S, et al. Comparative evaluation of big-data systems on scientific image analytics workloads. Proceedings of the VLDB Endowment 2017; 10(11): 1226-1237.  doi: 10.14778/3137628.3137634
  * 15 Paraskevakos I, Luckow A, Khoshlessan M, et al. Task-Parallel Analysis of Molecular Dynamics Trajectories. In: ICPP 2018. Proceedings of the 47th International Conference on Parallel Processing. Association for Computing Machinery; 2018; New York, NY, USA
  * 16 Dask Bag Documentation: Known Limitations. https://docs.dask.org/en/latest/bag.html#known-limitations; .  Accessed: 2022-08-17.
  * 17 Dean J, Ghemawat S. MapReduce: simplified data processing on large clusters. Communications of the ACM 2008; 51(1): 107–113.
  * 18 Zaharia M, Chowdhury M, Das T, et al. Resilient Distributed Datasets: A Fault-tolerant Abstraction for In-memory Cluster Computing. In: NSDI’12. Proceedings of the 9th USENIX Conference on Networked Systems Design and Implementation. USENIX Association; 2012; Berkeley, CA, USA: 2–2.
  * 19 Dask SLURM Cluster Documentation. https://jobqueue.dask.org/en/latest/generated/dask_jobqueue.SLURMCluster.html#dask-jobqueue-slurmcluster; .  Accessed: 2022-08-23.
  * 20 Amunts K, Lepage C, Borgeat L, et al. BigBrain: An Ultrahigh-Resolution 3D Human Brain Model. Science 2013; 340(6139): 1472–1475.  doi: 10.1126/science.1235381
  * 21 SAM Library. https://github.com/big-data-lab-team/sam; .  Accessed: 2022-09-12.
  * 22 Zuo XN, Anderson JS, Bellec P, et al. An open science resource for establishing reliability and reproducibility in functional connectomics. Scientific data 2014; 1: 140049.
  * 23 Datalad dataset. https://datasets.datalad.org/?dir=/corr/RawDataBIDS; .  Accessed: 2022-09-12.
  * 24 Gorgolewski KJ, Alfaro-Almagro F, Auer T, et al. BIDS apps: Improving ease of use, accessibility, and reproducibility of neuroimaging data analysis methods. PLoS computational biology 2017; 13(3): e1005209.
  * 25 Slaughter E, Wu W, Fu Y, et al. Task Bench: A Parameterized Benchmark for Evaluating Parallel Runtime Performance. In: SC20: International Conference for High Performance Computing, Networking, Storage and Analysis. ; 2020: 1-15


# Gluon Sivers function from forward exclusive $\chi_{c1}$ photoproduction on
unpolarized protons

Sanjin Benić  Department of Physics, Faculty of Science, University of Zagreb,
Bijenička c. 32, 10000 Zagreb, Croatia  Adrian Dumitru  Department of Natural
Sciences, Baruch College, CUNY, 17 Lexington Avenue, New York, NY 10010, USA
The Graduate School and University Center, The City University of New York,
365 Fifth Avenue, New York, NY 10016, USA  Leszek Motyka  Jagiellonian
University, Institute of Theoretical Physics, Łojasiewicza 11, 30-348 Kraków,
Poland  Tomasz Stebel  Jagiellonian University, Institute of Theoretical
Physics, Łojasiewicza 11, 30-348 Kraków, Poland

###### Abstract

Exclusive production of a $\chi_{c1}$ axial vector quarkonia in photon-proton
scattering at high energies requires a $C$-odd $t$-channel exchange. In the
limit of vanishing momentum transfer this occurs either via the exchange of a
photon, the Primakoff process, where the spin of the proton does not change.
For axial-vector meson production, as a consequence of the Landau-Yang
theorem, the Primakoff cross section is finite as $t\to 0$. Alternatively, a
$C$-odd spin dependent Odderon can be exchanged, which involves a spin flip of
the proton. The resulting cross section is related to the square of the gluon
Sivers function or its $\boldsymbol{k}_{\perp}^{2}$-moment. Using two models
for the gluon Sivers function from the literature we compute the ratio of
Sivers to Primakoff cross sections and the angular coefficient
$\lambda_{\theta}$ governing the angular distribution of the $\chi_{c1}\to
J/\psi+\gamma$ decay as functions of $x$. We point out that these observables
constrain the magnitude of the gluon Sivers function at small $x$ which could
be accessed in electron-proton scattering and ultraperipheral proton-proton
and nucleus-proton collisions.

††preprint: ZTF-EP-24-06

Introduction: In our previous work [1] we have computed the exclusive cross
section of axial vector $\chi_{c1}$ quarkonia:
$\gamma^{*}(q)p(P)\to\chi_{c1}(\Delta)p(P^{\prime})$. This arises from tri-
gluon (or, Odderon) exchange with the proton, interfering with the photon
exchange (the Primakoff process). The focus was on the so-called spin-
independent Odderon [2] which is important at finite momentum transfer
$t=(P-P^{\prime})^{2}\approx-\boldsymbol{\Delta}_{\perp}^{2}$. For $t\to 0$,
on the other hand, the corresponding cross section vanishes due to odd parity
of the spin independent Odderon amplitude under
$\boldsymbol{\Delta}_{\perp}\to-\boldsymbol{\Delta}_{\perp}$.

In near-forward kinematics the Primakoff process becomes the major background
to Odderon searches. Nevertheless, the $\chi_{c1}$ is special because, thanks
to the Landau-Yang theorem [3, 4], the near-forward Coulomb tail ($\sim
1/|t|$) of the Primakoff cross section, in the photoproduction
$q^{2}=-Q^{2}\to 0$ limit111We require here the scale hierarchy
$Q^{2}\ll|t|\lesssim 1/R^{2}_{p}$ where $R_{p}\simeq 1$ fm is the radius of
the proton., is screened [1], and the cross section remains finite. In the
non-relativistic (NR) $m_{c}\to\infty$ limit we simply have [1, 5]

$\lim_{t\to 0}\frac{\mathrm{d}\sigma_{\rm Prim}}{\mathrm{d}|t|}=\frac{3\pi
q_{c}^{4}\alpha^{3}N_{c}|R^{\prime}(0)|^{2}|F_{1}(0)|^{2}}{m_{c}^{9}}\,.$ (1)

Here $R^{\prime}(0)$ is the derivative of the non-relativistic (NR) wave
function at the origin, $m_{c}$ is the charm quark mass, $\alpha$ is the
electromagnetic fine structure constant, with $q_{c}=2/3$ the quark fractional
charge, and $N_{c}=3$ denotes the number of colors. The Dirac form-factor
$F_{1}(t)$ is evaluated at $t=0$. Note that this result is also independent of
the collision energy $W^{2}=(q+P)^{2}$ in the eikonal limit, $W^{2}\gg|t|$.

It has been pointed out in [2] (see also [6]) that the near-forward region of
the unpolarized $\gamma^{*}p$ cross section provides an opportunity to
constrain the spin-dependent Odderon [2, 7], that is, the gluon Sivers
function. The gluon Sivers function represents an azimuthally asymmetric
transverse momentum distribution (TMD) of gluons in the proton that is
expected in collisions off a transversely polarized proton [8, 9, 10]. This
function represents a key element in our understanding of the structure of the
proton; to date, it is rather poorly constrained, as will be clear from fig. 2
below. The determination of the gluon Sivers function at high energies
represents one of the key “golden measurements” to be performed at the
upcoming electron-ion collider (EIC) [11, 12, 13].

The traditional strategy for measuring the gluon Sivers function is through
single spin asymmetries (SSAs) in single-inclusive $ep^{\uparrow}$ [14, 15,
16, 17, 18, 19] or $pp^{\uparrow}$ collisions [20, 21, 22] with heavy quark
production. Considering instead exclusive $C$-even meson production, a finite
contribution arises to the unpolarized $ep$ cross section from the gluon
Sivers function when $t\to 0$ [2]. The simple reason is the possibility of a
proton helicity flip as associated with the Sivers function; this is
illustrated by the lower diagram in Fig 1. We argue that measuring axial
quarkonia such as $\chi_{c1}$ for the same near-forward kinematics has a
unique advantage over other mesons ($\pi^{0}$ or $\eta$, but also $a_{2}$,
$f_{2}$,…) or even other quarkonia ($\eta_{c}$, $\chi_{c0}$, $\chi_{c2}$,…).
In these cases the Primakoff cross section blows up at $t\to 0$, and it will
always dominate the Sivers contribution. By contrast, for axial quarkonia such
as $\chi_{c1}$, the Primakoff process gives a finite contribution (1).
Depending on the relative magnitudes of the Sivers vs. Primakoff channels
there is potentially a chance to become directly sensitive to the gluon Sivers
function, without the need of subtracting the Primakoff background at all.

Figure 1: Representative diagrams illustrating the Primakoff (top) and the
Sivers (bottom) contributions to the exclusive
$\gamma(q)p(P)\to\chi_{c}(\Delta)p(P^{\prime})$ process. The arrows indicate
the helicity transfer in the respective channels.

We propose another promising opportunity to constrain the gluon Sivers
function via the polarization of $\chi_{c1}$ in the Sivers and the Primakoff
production channels, respectively. The argument is as follows. Since the
photon exchange at high energies cannot flip the proton helicity (at $t\to 0$
and $Q^{2}\to 0$), conservation of helicity dictates that the $\chi_{c1}$ from
the Primakoff channel must have transverse polarization. By contrast, the spin
dependent Odderon flips the proton helicity and so the produced $\chi_{c1}$
must have longitudinal polarization. This way, the polarization of $\chi_{c1}$
is directly sensitive to the gluon Sivers function. The different polarization
states can be accessed through the angular distribution of the $\chi_{c1}$
decay products for which we provide an example below. The helicity transfer
from the initial photon to the final $\chi_{c1}$ is illustrated for both
channels in fig. 1, with red arrows representing photon or $\chi_{c1}$
polarization vectors, and the blue arrows stand for proton helicity.

Exclusive $\gamma p\to\chi_{c1}p$ cross section: The helicity dependent
Odderon exchange amplitude from eq. (5) in ref. [1] is

$\begin{split}\mathcal{O}_{\lambda_{p}\lambda_{p}^{\prime}}&(\boldsymbol{k}_{\perp},\boldsymbol{\Delta}_{\perp})=\int\mathrm{d}^{2}\boldsymbol{r}_{\perp}\mathrm{d}^{2}\boldsymbol{b}_{\perp}\mathrm{e}^{-\mathrm{i}\boldsymbol{k}_{\perp}\cdot\boldsymbol{r}_{\perp}}\mathrm{e}^{-\mathrm{i}\boldsymbol{\Delta}_{\perp}\cdot\boldsymbol{b}_{\perp}}\\\
&\times\frac{1}{N_{c}}\frac{{\rm Im}\left\\{\langle
P^{\prime}\lambda^{\prime}_{p}|{\rm
tr}\left[V^{\dagger}(\boldsymbol{x}_{\perp})V(\boldsymbol{y}_{\perp})\right]|P\lambda_{p}\rangle\right\\}}{\langle
P\lambda_{p}|P\lambda_{p}\rangle}\,.\end{split}$ (2)

Here
$\boldsymbol{x}_{\perp}=\boldsymbol{b}_{\perp}+\frac{\boldsymbol{r}_{\perp}}{2}$,
$\boldsymbol{y}_{\perp}=\boldsymbol{b}_{\perp}-\frac{\boldsymbol{r}_{\perp}}{2}$
are transverse coordinates, and $V(\boldsymbol{x}_{\perp})$ are Wilson lines
describing scattering of the $c$ quarks off the gluon fields in the proton in
the eikonal (high-energy) limit. $\lambda_{p}$ ($\lambda_{p}^{\prime}$) are
the helicities of the incoming (outgoing) proton. The general decomposition of
the above off-forward proton matrix elements classifies into three different
odderons [2] in the eikonal limit. Focusing on forward kinematics the only
non-zero contribution arises from the spin dependent Odderon [16, 2]

$\begin{split}\mathcal{O}_{\lambda_{p}\lambda^{\prime}_{p}}&(\boldsymbol{k}_{\perp},0)=\lambda_{p}\delta_{\lambda_{p},-\lambda^{\prime}_{p}}(\boldsymbol{k}_{\perp}\times\boldsymbol{\epsilon}_{\perp}^{\lambda_{p}})\\\
&\times\frac{(2\pi)^{3}\mathrm{i}g^{2}}{2\sqrt{2}N_{c}}\frac{xf_{1T}^{\perp
g}(x,\boldsymbol{k}_{\perp})}{\boldsymbol{k}_{\perp}^{2}M_{p}}\,,\end{split}$
(3)

where
$\boldsymbol{\epsilon}_{\perp}^{\lambda}=(-\lambda,-\mathrm{i})/\sqrt{2}$ and
$f_{1T}^{\perp g}(x,\boldsymbol{k}_{\perp})$ is the gluon Sivers function that
is associated with the spin flip of the proton. $M_{p}$ is the proton mass and
$g$ is the strong coupling constant. Finally, $x=q\cdot(P-P^{\prime})/q\cdot
P$ denotes the usual momentum fraction in diffractive processes.

With eq. (3) for
$\mathcal{O}_{\lambda_{p}\lambda^{\prime}_{p}}(\boldsymbol{k}_{\perp},0)$ we
now take the photoproduction, $Q^{2}\to 0$, as well as the forward,
$\boldsymbol{\Delta}_{\perp}\to 0$, limits of the $\chi_{c1}$ production
amplitude of ref. [1]. As a consequence of helicity conservation only the
amplitude for longitudinally polarized $\chi_{c1}$ survives,
$\lambda_{\chi}=0$. It is given by [1]

$\begin{split}&\mathcal{M}_{\lambda_{p}\lambda_{\gamma},\lambda^{\prime}_{p}\lambda_{\chi}}=\frac{8\mathrm{i}N_{c}eq_{c}}{M_{\chi}}\delta_{\lambda_{\chi}0}\\\
&\times\int\frac{\mathrm{d}^{2}\boldsymbol{k}_{\perp}}{(2\pi)^{2}}\mathcal{O}_{\lambda_{p}\lambda_{p}^{\prime}}(\boldsymbol{k}_{\perp},0)\int_{0}^{1}\frac{\mathrm{d}z}{4\pi}\int\frac{\mathrm{d}^{2}\boldsymbol{l}_{\perp}}{(2\pi)^{2}}\frac{\phi_{\chi}(\boldsymbol{l}_{\perp},z)}{z\bar{z}}\\\
&\times\frac{(\boldsymbol{l}_{\perp}^{2}+m_{c}^{2})(\boldsymbol{\epsilon}_{\perp}^{\lambda_{\gamma}}\times\boldsymbol{l}_{\perp})+\boldsymbol{l}_{\perp}^{2}(\boldsymbol{\epsilon}_{\perp}^{\lambda_{\gamma}}\times\boldsymbol{k}_{\perp})}{(\boldsymbol{l}_{\perp}+\boldsymbol{k}_{\perp})^{2}+m_{c}^{2}}\,,\end{split}$
(4)

where $\lambda_{\gamma}$ is the incoming photon helicity,
$\phi_{\chi}(\boldsymbol{l}_{\perp},z)$ is the $\chi_{c1}$ light-cone wave
function, and $M_{\chi}$ its mass. Parameterizing the overall helicity
structure as
$\mathcal{M}_{\lambda_{p}\lambda_{\gamma},\lambda^{\prime}_{p}\lambda_{\chi}}=\lambda_{p}\delta_{\lambda_{p},-\lambda^{\prime}_{p}}\delta_{\lambda_{\gamma},-\lambda_{p}}\delta_{\lambda_{\chi}0}\mathcal{M}_{\rm
Siv}$, the cross section is

$\lim_{t\to 0}\frac{\mathrm{d}\sigma_{\rm
Siv}}{\mathrm{d}|t|}=\frac{1}{32\pi}|\mathcal{M}_{\rm Siv}|^{2}\,.$ (5)

Taking the collinear limit on the proton side ($\boldsymbol{k}_{\perp}^{2}\ll
m_{c}^{2}$) and the NR limit ($\boldsymbol{l}_{\perp}^{2}\ll m_{c}^{2}$), the
cross section takes a remarkably simple form

$\lim_{t\to 0}\frac{\mathrm{d}\sigma_{\rm
Siv}}{\mathrm{d}|t|}=\frac{3\pi^{3}q_{c}^{2}\alpha\alpha_{S}^{2}M_{p}^{2}|R^{\prime}(0)|^{2}|xf_{1T}^{\perp(1)g}(x)|^{2}}{N_{c}m_{c}^{11}}\,,$
(6)

where $xf_{1T}^{\perp(1)g}(x)$ is the $\boldsymbol{k}_{\perp}^{2}$ moment of
the gluon Sivers function222The connection to the so-called $d$-type twist-3
gluonic Qiu-Sterman functions [23] was established in [16] at small-$x$.

$xf_{1T}^{\perp(1)g}(x,\mu^{2})\equiv\int^{\mu^{2}}\mathrm{d}^{2}\boldsymbol{k}_{\perp}\frac{\boldsymbol{k}_{\perp}^{2}}{2M_{p}^{2}}xf_{1T}^{\perp
g}(x,\boldsymbol{k}_{\perp})\,,$ (7)

with $\mu^{2}$ the factorization scale. Eq. (6) parallels the classic result
for the exclusive $J/\psi$ production cross section at $t\to 0$ which is
proportional to the square of the unpolarized gluon distribution function
[24]. Comparing to the Primakoff cross section in (1), eq. (6) is suppressed
by two additional powers of $m_{c}$ due to the twist-3 nature of the Sivers
function.

Figure 2: Ratio of the $\gamma p\to\chi_{c1}p$ cross sections in the Sivers to
Primakoff channels, in the $t\to 0$ limit, as a function of diffractive $x$.
These processes involve $C$-odd tri-gluon (Odderon) and photon exchanges,
respectively.

In practice only the sum of the Sivers and the Primakoff cross sections can be
measured. Nevertheless, we argue that the two contributions can be
disentangled as follows. Firstly, thanks to the spin flip nature of the Sivers
function there is no interference of the amplitudes. Secondly, at high energy
only the gluon Sivers function is expected to exhibit $x$ dependence; recall
that (1) is completely independent of $x$. Whether the $x$ dependence in the
total cross section would be significant or not depends, of course, on the
relative magnitudes of the Primakoff and Sivers components. It is useful to
consider the ratio

$r\equiv\left(\frac{\mathrm{d}\sigma_{\rm
Siv}}{\mathrm{d}|t|}\Big{/}\frac{\mathrm{d}\sigma_{\rm
Prim}}{\mathrm{d}|t|}\right)_{t=0}\,,$ (8)

where the Sivers cross section is in general given by eqs. (4) and (5), and
the Primakoff cross section follows from eqs. (33) and (43) in [1]. The NR
limit is again very illustrative

$r\to\frac{4\pi^{2}}{q_{c}^{2}N_{c}^{2}}\frac{\alpha_{S}^{2}}{\alpha^{2}}\frac{M^{2}_{p}}{M^{2}_{\chi}}|xf_{1T}^{\perp(1)g}(x)|^{2}\,,$
(9)

with $R^{\prime}(0)$ cancelling in the ratio! In general, we can expect the
uncertainty in modelling the non-perturbative $\chi_{c1}$ wave function to
partially cancel in $r$. Here $M_{\chi}=2m_{c}$ is the mass of $\chi_{c1}$ in
the NR limit.

In fig. 2 we show numerical results for (9) using two models for
$xf_{1T}^{\perp g}(x,\boldsymbol{k}_{\perp})$. The model labelled “SIDIS2” and
shown with the dashed blue curve was obtained from a fit to transverse single
spin asymmetries measured at RHIC [21]. Its first Sivers moment is [21]

$|f_{1T}^{\perp(1)g}(x)|=\sqrt{\frac{\mathrm{e}}{2}}\rho^{2}\sqrt{\frac{1-\rho}{\rho}}\frac{\sqrt{\langle
k_{\perp}^{2}\rangle}}{M_{p}}\mathcal{N}_{g}(x)g(x)\,,$ (10)

where
$\mathcal{N}_{g}(x)=N_{g}x^{\alpha}(1-x)^{\beta}(\alpha+\beta)^{\alpha+\beta}/(\alpha^{\alpha}\beta^{\beta})$
with $N_{g}=0.05$, $\alpha=0.8$, $\beta=1.4$, $\rho=0.576$, $\langle
k_{\perp}^{2}\rangle=0.25$ GeV2, and $g(x)$ is the unpolarized gluon PDF. The
label “5% positivity” (full magenta curve) represents a model [18] where the
magnitude of the Sivers function is set at $5\%$ of the function that
saturates the positivity bound within the model of [20]:

$\begin{split}|f_{1T}^{\perp(1)g}(x)|&=0.05\times\frac{\sigma}{M_{p}}\left[1+\rho\mathrm{e}^{\rho}{\rm
Ei}(-\rho)\right]g(x)\,,\end{split}$ (11)

where ${\rm Ei}(\rho)$ is the exponential integral, with
$\rho=\sigma^{2}/\langle k_{\perp}^{2}\rangle$. In the numerical computation
$\sigma^{2}=\langle k_{\perp}^{2}\rangle$ and $\sigma=0.8$ GeV [20]. In both
models the factorization scale in $g(x,\mu^{2})$ is taken as
$\mu^{2}=4m_{c}^{2}$.

The difference between the two models in fig. 2 illustrates the large
uncertainty in our current knowledge of $f_{1T}^{\perp
g}(x,\boldsymbol{k}_{\perp})$. The numerical values of $r$ further indicate a
potentially large contribution from the Sivers component to the total cross
section, emphasizing that a measurement of the exclusive $\gamma
p\to\chi_{c1}p$ cross section, in particular of its $x$ dependence, can be
used to constrain the gluon Sivers function. A computation based on eqs. (4)
and (5) gives a factor of 2.0 (3.0) lower Sivers cross section for the SIDIS2
(5% positivity) models. For the Primakoff cross section the analogous ratio is
1.12.

It is interesting to contrast the above two phenomenological models to a
small-x model for the spin dependent Odderon [16, 25, 26] which follows from a
group theory constraint. The trace properties of the special unitary matrix in
(2) restrict the magnitude of the Odderon to
$|\mathcal{O}(\boldsymbol{r}_{\perp})|\leq\mathcal{N}(\boldsymbol{r}_{\perp})^{3/2}/3$,
in the weak field limit where the $C$-even (Pomeron) amplitude
$\mathcal{N}(\boldsymbol{r}_{\perp})$ is small. The form and magnitude of
$\mathcal{N}(\boldsymbol{r}_{\perp})$ at the initial condition $x_{0}\simeq
0.01$ of the small-$x$ dipole model [27, 28, 29, 16] are fairly well
constrained [30, 31, 32, 33, 34]. Relating
$\mathcal{N}(\boldsymbol{r}_{\perp})$ (in the double-log regime) to the
leading-twist collinear gluon density $xg(x)$ [35, 30] we obtain that this
“maximal Odderon” corresponds to the Sivers moment

$|f_{1T}^{\perp(1)g}(x)|=\frac{1}{3}\frac{Q_{S}(x)}{M_{p}}g(x)~{}.$ (12)

$Q_{S}(x)$ denotes the intrinsic transverse momentum scale of the proton at
$x$. With $Q_{S}(x_{0})=0.25-0.5$ GeV, the upper limit for
$|f_{1T}^{\perp(1)g}(x_{0})|$ reaches $9-18\%$ of the unpolarized gluon
density. Comparing to (11), but with 100% of the positivity bound (at the TMD
level), $f^{\perp(1)g}_{1T}(x)$ is about 35% of $g(x)$. Interestingly, the
group theory bound, valid at small-$x$, is more stringent than the positivity
bound.

The above models do not represent solutions of the small-$x$ evolution
equations. Rather, they should be understood as models of the initial
condition. In contrast to the Pomeron, the small-$x$ linear evolution of the
Odderon [36] leads to a roughly constant energy dependence, $\sim 1/x^{0}$,
and suppressed by non-linear unitarity corrections [27, 37, 29, 25, 26, 38].

At this point it is useful to provide an estimate of the magnitude of the
Primakoff cross section for the EIC but also for the ultra-peripheral
collisions (UPCs) at RHIC and the LHC. Taking $|R^{\prime}(0)|^{2}=0.1296$
GeV5 [39] we obtain from (1) that $(\mathrm{d}\sigma_{\rm
Prim}/\mathrm{d}|t|)_{t=0}\approx 0.69$ pb/GeV2. For $pA$ UPCs we convolute
the $\gamma-p$ cross section with the photon flux factor [40] and obtain
$0.05$ nb/GeV2 ($0.79$ nb/GeV2) at $\sqrt{s}=200$ GeV ($5.02$ TeV) RHIC (LHC)
energy.

Figure 3: Angular coefficient $\lambda_{\theta}$ as a function of $x$ based on
different models of the gluon Sivers function (see text).

Angular distributions in $\chi_{c1}\to J/\psi\gamma$ decay: As already
mentioned, the Primakoff process produces transversely polarized $\chi_{c1}$’s
while those from the Sivers channel are longitudinally polarized. In order to
distinguish the polarization we focus on the dominant decay channel,
$\chi_{c1}\to J/\psi+\gamma$, assuming a pure electric dipole (E1) transition,
for simplicity. The angular distribution of $J/\psi$’s is [41, 42, 43]

$W(\theta,\phi)\propto\frac{N}{3+\lambda_{\theta}}\big{(}1+\lambda_{\theta}\cos^{2}\theta+\dots\big{)}\,.$
(13)

The ellipses in (13) correspond to four additional angular coefficients
$\lambda_{i}$ that drop out by explicit calculation. The polarization angles
$\theta$ and $\phi$ are frame-dependent; we assume the Collins-Soper frame.
Since our calculation [1] has been performed in the $\gamma p$ center-of-mass
frame, with the helicity quantized along the $z$-direction, we have confirmed
that the general choice of the Collins-Soper polarization vectors [44]
collapses to that in [1] when $t\to 0$. The simple reason is that the Collins-
Soper frame and the $\gamma p$ center-of-mass frame are related by a
longitudinal boost when $t\to 0$. All the angular coefficients $\lambda_{i}$
can be expressed in terms of the spin density matrix
$\rho_{\lambda^{\prime}_{\chi},\lambda_{\chi}}\propto\sum_{\lambda_{\gamma}\lambda_{p}\lambda^{\prime}_{p}}\mathcal{M}_{\lambda_{p}\lambda_{\gamma},\lambda_{p}^{\prime}\lambda^{\prime}_{\chi}}^{*}\mathcal{M}_{\lambda_{p}\lambda_{\gamma},\lambda_{p}^{\prime}\lambda_{\chi}}$.
The $\lambda_{\theta}$ coefficient is given by [41, 42, 43]

$\lambda_{\theta}=\frac{3\rho_{0,0}-N}{3N-\rho_{0,0}}\,,$ (14)

where $N=\rho_{1,1}+\rho_{0,0}+\rho_{-1,-1}$. The remaining angular
coefficients are all proportional to off-diagional entries of the density
matrix (unlike $\lambda_{\theta}$) which are zero in both the Sivers and the
Primakoff channels. The diagonal entries yield

$\begin{split}&\rho_{0,0}\propto 16\pi\lim_{|t|\to
0}\frac{\mathrm{d}\sigma_{\rm Siv}}{\mathrm{d}|t|}\,,\\\
&\rho_{1,1}=\rho_{-1,-1}\propto 8\pi\lim_{|t|\to 0}\frac{\mathrm{d}\sigma_{\rm
Prim}}{\mathrm{d}|t|}\,,\end{split}$ (15)

with a common proportionality factor that cancels in (14). The result can be
conveniently written in terms of the ratio $r$ from (8)

$\lambda_{\theta}=\frac{2r-1}{2r+3}\,.$ (16)

In the NR limit $\lambda_{\theta}$ will be completely independent of
$|R^{\prime}(0)|$ and thus sensitive only to the gluon Sivers function. The
two interesting cases of (16) are $r\to 0$ (Primakoff contribution only with
transverse $\chi_{c1}$) and $1/r\to 0$ (Sivers contribution only with
longitudinal $\chi_{c1}$), where $\lambda_{\theta}\to-1/3$ and
$\lambda_{\theta}\to 1$, respectively. Hence, even the sign of
$\lambda_{\theta}$ could indicate the presence of the spin dependent Odderon.

In fig. 3 we show $\lambda_{\theta}(x)$ for the above-mentioned models for the
gluon Sivers functions. For the “SIDIS2” model the Sivers channel is small
(see fig. 2) and drops with $x$ ($r\ll 1$) and so $\lambda_{\theta}$ is
numerically almost at the pure Primakoff value of $-1/3$. The “5% positivity”
model shows the opposite behavior, with a slightly negative $\lambda_{\theta}$
at $x\leq 0.1$ which increases towards the pure Sivers value of $+1$ at $x\sim
10^{-3}$.

Discussion and conclusion: We suggest that the measurement of the exclusive
$\chi_{c1}$ production cross section and polarization in the forward limit of
high-energy unpolarized $\gamma-p$ collisions can provide stringent
constraints on the gluon Sivers function $f_{1T}^{\perp
g}(x,\boldsymbol{k}_{\perp})$. In fact, the cross section (6), being
proportional to the square of the first moment of the Sivers function in the
NR limit, is potentially even more sensitive than the traditional SSAs with a
linear dependence on $f_{1T}^{\perp g}(x,\boldsymbol{k}_{\perp})$. The gluon
Sivers function is a key element of our understanding of the proton and, to
date, is very poorly known. Based on the connection between the spin dependent
Odderon and $f_{1T}^{\perp g}(x,\boldsymbol{k}_{\perp})$ we compute the cross
section and find a non-zero contribution at $t\to 0$. Thanks to the Landau-
Yang theorem the background from the Primakoff channel is finite. Some current
models for the gluon Sivers function suggest that the two channels may have
cross sections of similar magnitude. By examining angular distributions of the
$\chi_{c1}$ decay products we find that the $\cos^{2}\theta$ angular
coefficient in the Primakoff and Sivers channels, has opposite sign. In
conjunction with similar measurements of $C$-odd $t$-channel exchanges at high
momentum transfer [1] the processes discussed here have great potential to
provide fundamental insight into the structure of the proton.

Acknowledgments: We thank Y. Hatta for useful comments on the manuscript. The
work of S. B. is supported by the Croatian Science Foundation (HRZZ) no. 5332
(UIP-2019-04). A. D. acknowledges support by the DOE Office of Nuclear Physics
through Grant DE-SC0002307. T. S. kindly acknowledges the support of the
Polish National Science Center (NCN) Grant No. 2021/43/D/ST2/03375. S. B.
acknowledges the hospitality of the EIC Theory Institute at the Brookhaven
National Laboratory in March 2024 where this work was initiated.

## References

  * Benić  _et al._ [2024] S. Benić, A. Dumitru, A. Kaushik, L. Motyka, and T. Stebel, Photon-Odderon interference in exclusive $\chi_{c}$ charmonium production at the Electron-Ion Collider, (2024), arXiv:2402.19134 [hep-ph] .
  * Boussarie  _et al._ [2020] R. Boussarie, Y. Hatta, L. Szymanowski, and S. Wallon, Probing the Gluon Sivers Function with an Unpolarized Target: GTMD Distributions and the Odderons, Phys. Rev. Lett. 124, 172501 (2020), arXiv:1912.08182 [hep-ph] .
  * Landau [1948] L. D. Landau, On the angular momentum of a system of two photons, Dokl. Akad. Nauk SSSR 60, 207 (1948).
  * Yang [1950] C.-N. Yang, Selection Rules for the Dematerialization of a Particle Into Two Photons, Phys. Rev. 77, 242 (1950).
  * Jia  _et al._ [2023] Y. Jia, Z. Mo, J. Pan, and J.-Y. Zhang, Photoproduction of C-even quarkonia at the EIC and EicC, Phys. Rev. D 108, 016015 (2023), arXiv:2207.14171 [hep-ph] .
  * Ma [2003] J. P. Ma, Diffractive photoproduction of eta(c), Nucl. Phys. A 727, 333 (2003), arXiv:hep-ph/0301155 .
  * Boer  _et al._ [2016] D. Boer, M. G. Echevarria, P. Mulders, and J. Zhou, Single spin asymmetries from a single Wilson loop, Phys. Rev. Lett. 116, 122001 (2016), arXiv:1511.03485 [hep-ph] .
  * Burkardt [2004] M. Burkardt, Sivers mechanism for gluons, Phys. Rev. D 69, 091501 (2004), arXiv:hep-ph/0402014 .
  * Meissner  _et al._ [2007] S. Meissner, A. Metz, and K. Goeke, Relations between generalized and transverse momentum dependent parton distributions, Phys. Rev. D 76, 034002 (2007), arXiv:hep-ph/0703176 .
  * Boer  _et al._ [2015] D. Boer, C. Lorcé, C. Pisano, and J. Zhou, The gluon Sivers distribution: status and future prospects, Adv. High Energy Phys. 2015, 371396 (2015), arXiv:1504.04332 [hep-ph] .
  * Pro [2020] _Proceedings, Probing Nucleons and Nuclei in High Energy Collisions: Dedicated to the Physics of the Electron Ion Collider: Seattle (WA), United States, October 1 - November 16, 2018_ (WSP, 2020) arXiv:2002.12333 [hep-ph] .
  * Abdul Khalek  _et al._ [2022] R. Abdul Khalek _et al._ , Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report, Nucl. Phys. A 1026, 122447 (2022), arXiv:2103.05419 [physics.ins-det] .
  * Achenbach  _et al._ [2024] P. Achenbach _et al._ , The present and future of QCD, Nucl. Phys. A 1047, 122874 (2024), arXiv:2303.02579 [hep-ph] .
  * Kang and Qiu [2008] Z.-B. Kang and J.-W. Qiu, Single transverse-spin asymmetry for D-meson production in semi-inclusive deep inelastic scattering, Phys. Rev. D 78, 034005 (2008), arXiv:0806.1970 [hep-ph] .
  * Beppu  _et al._ [2010] H. Beppu, Y. Koike, K. Tanaka, and S. Yoshida, Contribution of Twist-3 Multi-Gluon Correlation Functions to Single Spin Asymmetry in Semi-Inclusive Deep Inelastic Scattering, Phys. Rev. D 82, 054005 (2010), arXiv:1007.2034 [hep-ph] .
  * Zhou [2014] J. Zhou, Transverse single spin asymmetries at small x and the anomalous magnetic moment, Phys. Rev. D 89, 074050 (2014), arXiv:1308.5912 [hep-ph] .
  * Godbole  _et al._ [2018] R. M. Godbole, A. Kaushik, and A. Misra, Transverse single-spin asymmetry in the low-virtuality leptoproduction of open charm as a probe of the gluon Sivers function, Phys. Rev. D 97, 076001 (2018), arXiv:1709.03074 [hep-ph] .
  * Zheng  _et al._ [2018] L. Zheng, E. C. Aschenauer, J. H. Lee, B.-W. Xiao, and Z.-B. Yin, Accessing the gluon Sivers function at a future electron-ion collider, Phys. Rev. D 98, 034011 (2018), arXiv:1805.05290 [hep-ph] .
  * Chen  _et al._ [2023] L. Chen, H. Xing, and S. Yoshida, The twist-3 gluon contribution to Sivers asymmetry in $J/\psi$ production in semi-inclusive deep inelastic scattering, Phys. Rev. D 108, 054021 (2023), arXiv:2306.12647 [hep-ph] .
  * Anselmino  _et al._ [2004] M. Anselmino, M. Boglione, U. D’Alesio, E. Leader, and F. Murgia, Accessing Sivers gluon distribution via transverse single spin asymmetries in p(transv. polarized) p —$>$ D X processes at RHIC, Phys. Rev. D 70, 074025 (2004), arXiv:hep-ph/0407100 .
  * D’Alesio  _et al._ [2015] U. D’Alesio, F. Murgia, and C. Pisano, Towards a first estimate of the gluon Sivers function from AN data in pp collisions at RHIC, JHEP 09, 119, arXiv:1506.03078 [hep-ph] .
  * Kato  _et al._ [2024] N. Kato, L. Maxia, and C. Pisano, Spin asymmetries for $C$-even quarkonium production as a probe of gluon distributions, (2024), arXiv:2403.20017 [hep-ph] .
  * Ji [1992] X.-D. Ji, Gluon correlations in the transversely polarized nucleon, Phys. Lett. B 289, 137 (1992).
  * Ryskin  _et al._ [1997] M. G. Ryskin, R. G. Roberts, A. D. Martin, and E. M. Levin, Diffractive J / psi photoproduction as a probe of the gluon density, Z. Phys. C 76, 231 (1997), arXiv:hep-ph/9511228 .
  * Lappi  _et al._ [2016] T. Lappi, A. Ramnath, K. Rummukainen, and H. Weigert, JIMWLK evolution of the odderon, Phys. Rev. D 94, 054014 (2016), arXiv:1606.00551 [hep-ph] .
  * Yao  _et al._ [2019] X. Yao, Y. Hagiwara, and Y. Hatta, Computing the gluon Sivers function at small-$x$, Phys. Lett. B 790, 361 (2019), arXiv:1812.03959 [hep-ph] .
  * Kovchegov  _et al._ [2004] Y. V. Kovchegov, L. Szymanowski, and S. Wallon, Perturbative odderon in the dipole model, Phys. Lett. B 586, 267 (2004), arXiv:hep-ph/0309281 .
  * Jeon and Venugopalan [2005] S. Jeon and R. Venugopalan, A Classical Odderon in QCD at high energies, Phys. Rev. D 71, 125003 (2005), arXiv:hep-ph/0503219 .
  * Hatta  _et al._ [2005] Y. Hatta, E. Iancu, K. Itakura, and L. McLerran, Odderon in the color glass condensate, Nucl. Phys. A 760, 172 (2005), arXiv:hep-ph/0501171 .
  * Kowalski  _et al._ [2006] H. Kowalski, L. Motyka, and G. Watt, Exclusive diffractive processes at HERA within the dipole picture, Phys. Rev. D 74, 074016 (2006), arXiv:hep-ph/0606272 .
  * Albacete  _et al._ [2009] J. L. Albacete, N. Armesto, J. G. Milhano, and C. A. Salgado, Non-linear QCD meets data: A Global analysis of lepton-proton scattering with running coupling BK evolution, Phys. Rev. D 80, 034031 (2009), arXiv:0902.1112 [hep-ph] .
  * Lappi and Mäntysaari [2013] T. Lappi and H. Mäntysaari, Single inclusive particle production at high energy from HERA data to proton-nucleus collisions, Phys. Rev. D 88, 114020 (2013), arXiv:1309.6963 [hep-ph] .
  * Dumitru  _et al._ [2023] A. Dumitru, H. Mäntysaari, and R. Paatelainen, High-energy dipole scattering amplitude from evolution of low-energy proton light-cone wave functions, Phys. Rev. D 107, 114024 (2023), arXiv:2303.16339 [hep-ph] .
  * Casuga  _et al._ [2024] C. Casuga, M. Karhunen, and H. Mäntysaari, Inferring the initial condition for the Balitsky-Kovchegov equation, Phys. Rev. D 109, 054018 (2024), arXiv:2311.10491 [hep-ph] .
  * Bartels  _et al._ [2002] J. Bartels, K. J. Golec-Biernat, and H. Kowalski, A modification of the saturation model: DGLAP evolution, Phys. Rev. D 66, 014001 (2002), arXiv:hep-ph/0203258 .
  * Bartels  _et al._ [2000] J. Bartels, L. N. Lipatov, and G. P. Vacca, A New odderon solution in perturbative QCD, Phys. Lett. B 477, 178 (2000), arXiv:hep-ph/9912423 .
  * Motyka [2006] L. Motyka, Nonlinear evolution of pomeron and odderon in momentum space, Phys. Lett. B 637, 185 (2006), arXiv:hep-ph/0509270 .
  * Benić  _et al._ [2023] S. Benić, D. Horvatić, A. Kaushik, and E. A. Vivoda, Exclusive $\eta$c production from small-x evolved Odderon at an electron-ion collider, Phys. Rev. D 108, 074005 (2023), arXiv:2306.10626 [hep-ph] .
  * Eichten and Quigg [2019] E. J. Eichten and C. Quigg, Quarkonium wave functions at the origin: an update, (2019), arXiv:1904.11542 [hep-ph] .
  * Baltz [2008] A. J. Baltz, The Physics of Ultraperipheral Collisions at the LHC, Phys. Rept. 458, 1 (2008), arXiv:0706.3356 [nucl-ex] .
  * Faccioli  _et al._ [2010] P. Faccioli, C. Lourenco, and J. Seixas, A New approach to quarkonium polarization studies, Phys. Rev. D 81, 111502 (2010), arXiv:1005.2855 [hep-ph] .
  * Faccioli  _et al._ [2011] P. Faccioli, C. Lourenco, J. Seixas, and H. K. Wohri, Determination of $\chi_{c}$ and $\chi_{b}$ polarizations from dilepton angular distributions in radiative decays, Phys. Rev. D 83, 096001 (2011), arXiv:1103.4882 [hep-ph] .
  * Shao and Chao [2014] H.-S. Shao and K.-T. Chao, Spin correlations in polarizations of P-wave charmonia $\chi_{cJ}$ and impact on $J/\psi$ polarization, Phys. Rev. D 90, 014002 (2014), arXiv:1209.4610 [hep-ph] .
  * Beneke  _et al._ [1998] M. Beneke, M. Kramer, and M. Vanttinen, Inelastic photoproduction of polarized J / psi, Phys. Rev. D 57, 4258 (1998), arXiv:hep-ph/9709376 .