# A Novel GAN Approach to Augment Limited Tabular Data for Short-Term
Substance Use Prediction

Nguyen Thach1 Patrick Habecker1 Bergen Johnston1 Lillianna Cervantes1 Anika
Eisenbraun1 Alex Mason1 Kimberly Tyler1 Bilal Khan2&Hau Chan1 1University of
Nebraska-Lincoln
2Lehigh University<EMAIL_ADDRESS>{phabecker2, bergenjohnston,
lcervantes4, aeisenbraun3, amason19, ktyler2<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Substance use is a global issue that negatively impacts millions of persons
who use drugs (PWUDs). In practice, identifying vulnerable PWUDs for efficient
allocation of appropriate resources is challenging due to their complex use
patterns (e.g., their tendency to change usage within months) and the high
acquisition costs for collecting PWUD-focused substance use data. Thus, there
has been a paucity of machine learning models for accurately predicting short-
term substance use behaviors of PWUDs. In this paper, using longitudinal
survey data of 258 PWUDs in the U.S. Great Plains collected by our team, we
design a novel GAN that deals with high-dimensional low-sample-size tabular
data and survey skip logic to augment existing data to improve classification
models’ prediction on (A) whether the PWUDs would increase usage and (B) at
which ordinal frequency they would use a particular drug within the next 12
months. Our evaluation results show that, when trained on augmented data from
our proposed GAN, the classification models improve their predictive
performance (AUROC) by up to 13.4% in Problem (A) and 15.8% in Problem (B) for
usage of marijuana, meth, amphetamines, and cocaine, which outperform state-
of-the-art generative models.

##  1 Introduction

Substance use can create short- and long-term negative consequences for
persons who use drugs (PWUDs) HHS (2023); Newcomb and Bentler (1989); NIDA
(2020). These consequences include mental illness, HIV/AIDS, hepatitis, drug
overdose, and death. Within the U.S. alone, an estimated 161.8 million people
aged 12 or older used a substance (out of which 40.0 million used an illicit
drug) in the past month before being interviewed in 2021 SAMHSA (2022).
Furthermore, according to NIDA (2023), substance-involved overdose deaths,
including those related to illicit drugs and prescription opioids, continue to
increase over the years, with 106,699 deaths in 2021 compared to 91,799 (+16%)
in 2020 and 70,630 (+51%) in 2019. This alarming trend also applies globally,
with the estimated number of PWUDs in the past 12 months reaching 296 million
in 2021 from 240 million in 2011 (out of which 39.5 million and 27.3 million
had drug use disorders, respectively) United Nations (2023b).

In response to these large-scale negative impacts, various public and private
organizations around the globe have prompted initiatives for preventing and
reducing substance use at both population and individual levels. Prominent
examples are the United Nation’s Sustainability Goal 3: “Strengthen the
prevention and treatment of substance abuse” United Nations (2023a) and U.S.
Department of Health and Human Services (HHS)’s Healthy People 2030: “Reduce
misuse of drugs and alcohol” ODPHP (2020).

Generally, approaches toward these initiatives focus on designing and
deploying intervention and outreach programs/resources (e.g., rehabs and
consulting services) for PWUDs, with the main goal of reducing and eliminating
their usage of certain substances Ray et al. (2020); Colledge-Frisby et al.
(2023). While these programs and resources have shown to be effective to some
extent Ouimette et al. (1997); Tanner-Smith et al. (2016), they often require
volunteer participation from PWUDs, who face the difficulties and reluctances
of (self-)evaluating and (self-)determining whether they want or need help
Russell et al. (2021); Wilson and Brown (2023). Even when PWUDs agree to use
these programs/resources, they may have already experienced prior harms such
as overdose and mental illness SAMHSA (2021); Andersson et al. (2023).
Therefore, it is important to prevent harm from occurring in the first place
by carefully identifying PWUDs at the highest risk (i.e., those who are prone
to drastically increase usage of some drug) and allocating them appropriate
resources to reduce or eliminate potential harms.

Unfortunately, forecasting individual substance use behavior is challenging
due to its complex patterns (e.g., drug use frequency and co-use of multiple
drugs) and tendency to change over time at short timescales (i.e., within
months) Karamouzian et al. (2022); Linden-Carmichael et al. (2022); Lorvick et
al. (2023), as well as the lack of appropriate data and models for predicting
short-term future substance use (see Existing Efforts and Limitations and
additional contents in the full paper attached in the supplementary material).

Thus, our goal is to design accurate predictive models for modeling short-term
drug usage (i.e., within months) from PWUDs to aid healthcare agencies, local
communities, policymakers, and other stakeholders in the efficient allocation
of resources to PWUDs who need the most help. An effective predictive model
has the potential to improve the well-being of millions of PWUDs and hampers
the ongoing rapid growth of drug use prevalence and drug overdose death rates.

##### Our Approach and Associated Challenges.

To address the above shortcomings of substance use data and models, we formed
a collaboration between computer scientists and domain experts in substance
use research, including social scientists, intervention specialists, and
survey interviewers from University of Nebraska-Lincoln and Rural Drug
Addiction Research Center. With IRB approval, our team recruited a local
sample of 258 PWUDs in the Great Plains of the U.S. from which longitudinal
survey data were collected and stored in tabular format (detailed below in
Section 2.1). Despite its high value and relevance to our stated goal, the
sample size is small due to various challenges during data collection (mainly
disruptions due to COVID-19 and the lack of visible harm reduction movements
and punitive laws in the region), As a result, our preliminary predictive
models trained on currently available data achieved subpar performance with
respect to the baseline (demonstrated in our full paper). These models seem to
overfit due to the limited training examples ($<$200) with respect to the
number of (preprocessed) features ($>$600).

One effective way to tackle overfitting on small datasets is via data
augmentation Shorten and Khoshgoftaar (2019) i.e., creating synthetic samples
based on real data to increase its sample size. _Deep generative models_ such
as the popular Generative Adversarial Network (GAN) Goodfellow et al. (2014)
provide powerful tools for this purpose due to their flexibility in
representing complex and high-dimensional data distributions. In recent years,
these models have been extended to generate tabular data, especially
healthcare records, using specialized architectures (e.g., GOGGLE Liu et al.
(2022)) and/or novel training algorithms (CTGAN Xu et al. (2019)), which in
turn achieved promising performance across different evaluation criteria.

However, these models were benchmarked mostly on moderate to large datasets
that typically contain at least thousands of training examples. In rare cases
when being applied on small datasets with less than 1000 samples, e.g.,
Diabetes (768) and Breast (569) Street et al. (1993) as reported respectively
in the work for GOGGLE and TabDDPM Kotelnikov et al. (2023), the proposed
models and its selected comparatives often either fail to or narrowly
outperform simple baselines such as Bayesian networks Pearl (2011) and SMOTE
Chawla et al. (2002) in terms of their considered evaluation metrics. The
benchmark datasets for state-of-the-art models also contain no more than 200
features, which is substantially smaller than the size of our feature set
($>$600 after preprocessing). Additionally, the survey used for collecting our
tabular data contains _skip logic_ , a commonly employed functionality in
survey design for social science applications Fowler (1995); Dillman et al.
(2014), which has not been addressed by the aforementioned state-of-the-art
models.

##### Our Contributions.

In this paper, we design a novel specialized generative model to augment our
tabularized survey data with skip logic in order to improve the predictive
performance of our classification models, which ultimately predict the
following two short-term substance use behaviors: for a given PWUD and a
certain drug, (A) whether they would increase its usage and (B) at which
frequency (on a pre-defined ordinal scale) they would use it within the next
12 months. We summarize our contributions as follows:

  1. (I)

We believe ours is the first work that addresses _skip logic_ from surveys in
tabular data generation. We demonstrate its practical value by showing both
conceptually and empirically how enforcing constraints from skip logic (or
_skip constraints_ for brevity) positively affects the training of our
generative model. We are also one of the first to investigate the real-world
feasibility of deep generative models in settings where the number of features
exceed the sample size i.e., high-dimension low-sample-size (HDLSS). (See
Related Work in our full paper.)

  2. (II)

We design a novel GAN that deals with small tabular data containing $<$258
samples and 210 features (209 plus one target variable). Specifically, we
leverage CTGAN Xu et al. (2019), a well-known tabular GAN, by incorporating an
auxiliary classifier Park et al. (2018); Zhao et al. (2022) within its
architecture to generate high-quality samples conditioned on the corresponding
target variable. Since the transformed feature space has high dimensionality
(over 600), we also embed a global feature selection mechanism while training
the auxiliary classifier’s by employing the novel approach from Margeloiu et
al. (2023) that was shown to perform well on even smaller and higher
dimensional data than ours. Finally, to enforce the skip constraints stated in
(i), we take advantage of CTGAN’s built-in conditional vector.

  3. (III)

We implement and train the proposed GAN and use it to augment our (training)
data. The augmented data is then used to train binary and multiclass
classification models for predictive problems (A) and (B), respectively, for
each of the following drugs: marijuana, methamphetamine, amphetamines, and
cocaine. Our experimental results show that the average Area under the
Receiver Operating Characteristic curve (AUROC) evaluated on multiple distinct
sets of test data is improved by up to 13.4% in Problem (A) and 15.8% in
Problem (B) when the data is augmented, which is significantly higher than
what yielded using state-of-the-art generative models.

##  2 Problem Description

###  2.1 Background

##### Our Tabular Data.

The recruitment of PWUDs started in 2019 under the respondent-driven sampling
scheme Heckathorn (2014) in the Great Plains of the U.S. and has continued to
the present time. Enrolled PWUDs were followed up within 4–12 months after
their initial visit and took the same survey as before. In the survey, each
PWUD answered questions on a computer regarding their individual attributes,
including the drug use behavior of 18 different drugs. Use frequencies of
considered drugs (e.g., marijuana, cocaine, amphetamines, and methamphetamine)
were inquired on an ordinal scale (1–8) of {never, less than once a month,
once a month, once a week, 2–6 times a week, once a day, 2–3 times a day, 4 or
more times a day}. Collectively, the responses to these questions form the
(raw) features stored in a 2-D table.There are 258 samples (each represented
as a row) and 151 features (each represented as a column) in total. (See
supplementary material for the survey questions.) After being preprocessed
(detailed in Section 4.1), the table contains 209 features (plus the target
variable for teaching the desired classification models): 2 continuous and 207
categorical.

##### Definition 1 (Skip Logic).

In survey design, skip logic is a set of automated navigational rules that
allows respondents to skip to relevant questions depending on their prior
answers Couper (2008). Each of the rules is defined as a _skip constraint_ ,
which restricts the possible values of a subset of features $A$ in accordance
with the value of some feature $a$. We say the skip constraint on the _chain_
$A$ is _imposed_ by $a$ and denote this as $a\rightarrow A$. For the example
in Figure 1, the skip constraint on $A=$ {TB4} is imposed by $a=$ TB3 i.e.,
TB3 $\rightarrow$ {TB4}, in which TB4 is _omissible_ when TB3=“No”.

TB3 In the past 6 months, have you smoked cigarettes?

  1. No (0) — SKIP TO TB5

  2. Yes (1)

TB4 How many cigarettes do you usually smoke in a day?

  1. Not at all (1)

  2. Less than 1 cigarette a day (2)

  3. 1-5 cigarettes a day (3)

  4. Half a pack a day (4)

  5. A pack or more a day (5)

Figure 1: Skip Logic: If respondents answer “No” in TB3, they will
automatically be directed to TB5 without being asked on TB4.

##### Definition 2 (Tabular Data Generation).

Following the notations from Xu et al. (2019), given a table
$\mathbf{T}_{real}$ that is partitioned row-wise into training set
$\mathbf{T}_{train}$ and test set $\mathbf{T}_{test}$, the task involves
training a data generator $G$ on $\mathbf{T}_{train}$ and then independently
sampling rows using the learned $G$ to generate a synthetic table
$\mathbf{T}_{syn}$ such that $|\mathbf{T}_{syn}|=|\mathbf{T}_{train}|$ with
similar probability mass function for some target variable, $y$ (defined in
Section 2.2) (in order for us to fairly evaluate the efficacy of $G$).

Note that a capable $G$ can satisfy the latter requirement without affecting
the overall quality of the generated synthetic samples—that is, the features
in each generated row should be consistent with the label in that row. For our
setting, $\mathbf{T}_{real}$ and its derivatives consist of $N_{c}=2$
continuous features and $N_{d}=207+1$ categorical features.

##### Definition 3 (GAN and Its Extensions).

GANs are _deep_ generative models that have recently found success in modeling
tabular data Borisov et al. (2022) in addition to images and text. A GAN
typically consists of two separate networks: a generator $\mathcal{G}$ that
maps a noise distribution (typically Gaussian) to the data distribution and a
discriminator $\mathcal{D}$ that estimates the probability an input sample
came from the data distribution. The learning process is defined as an
adversarial game between $\mathcal{G}$ and $\mathcal{D}$ in which
$\mathcal{G}$ attempts to consistently fool $\mathcal{D}$ Goodfellow et al.
(2014).

To stabilize the training of GANs, Arjovsky et al. (2017) introduce WGAN that
provides a meaningful loss, the Wasserstein distance, for quantifying the
difference between the generated and real data distributions. Using the value
function from WGAN-GP Gulrajani et al. (2017) (an improvement of WGAN), Xu et
al. (2019) design CTGAN that aims to tackle class imbalance in categorical
features of tabular data by modifying $\mathcal{G}$ to additionally take a
vector as input. This so-called _conditional vector_ represents a certain
class/category of some categorical feature and is used to condition both the
generated samples and the real training samples. The model can thus
efficiently learn proper conditional distributions for each feature. We
further describe the main components of CTGAN in Section 3.2.

###  2.2 Problem Formulation

In this work, we focus on _data augmentation_ by designing a novel GAN to
generate high-quality synthetic samples that resemble our tabular data with
small sample size and skip logic incorporated. The augmented data is then used
to train classifiers for predicting PWUDs’ usage of a given drug within the
next 12 months, including (A) whether they would increase its usage and (B) at
which ordinal frequency they would use it. The former is a binary
classification problem wherein PWUDs exhibiting an increase and non-increase
(i.e., decrease or unchanged) in usage belong to the positive and negative
class, respectively; the latter is a multiclass classification problem in
which PWUDs are labeled according to their ordinal usage of the corresponding
drug within 12 months after. We mainly concern with predicting usage of the
most prevalent drugs (i.e., used by at least half of the PWUDs in our tabular
data), which include marijuana, methamphetamine (or meth for brevity),
amphetamines, and cocaine.

There are two unique properties of our tabularized survey data that impede
existing generative and predictive models from achieving desirable
performance: being high-dimension low-sample-size (HDLSS) and the presence of
skip logic (whose impact is demonstrated in Challenges and Observations of the
full paper).

##  3 Our Proposed GAN

###  3.1 Overview

We focus on GAN architectures due to their prevalence in tabular data
generation literature Borisov et al. (2022) and their convenient properties
(e.g., flexibility for conditional generation) that help us approach the
discussed challenges systematically. Overall, we extend the popular CTGAN
(reviewed in Section 3.2) as follows:

  1. i)

Due to the large number of (categorical) features, it is difficult for the
conditional generator $\mathcal{G}$ to learn to generate samples that are
conditioned on a particular feature. Therefore, during the training of
$\mathcal{G}$, we raise the production of synthetic samples in order to ensure
adequate training for conditional generation on a wide range of features.
Furthermore, we prioritize the generation of synthetic samples that are
conditioned on the empirical distributions of the target variable.

  2. ii)

To enhance the quality of synthetic samples _conditioned on_ the target
variable $y$, we incorporate an _auxiliary classifier_ $\mathcal{C}$ Park et
al. (2018) into the architecture of CTGAN for learning the correlation between
$y$ and other features. That is, $\mathcal{C}$ is trained to predict the label
of an input (synthetic or real) sample with $y$ removed.

  3. iii)

Within $\mathcal{C}$, which is originally a multilayer perceptron (MLP), we
add two auxiliary networks that together perform global feature selection and
compute the weights of the first hidden layer of $\mathcal{C}$ Margeloiu et
al. (2023). The objectives of this addition are twofold: facilitate effective
learning by $\mathcal{C}$ on HDLSS data and provide a mean to approximate the
feature importance scores.

  4. iv)

Lastly, we enforce skip logic on synthetic samples during training by
leveraging CTGAN’s conditional vector.

###  3.2 Conditional Generator

We first revisit the two main components of CTGAN, the _conditional generator_
and the _training-by-sampling_ procedure. Formally, the conditional generator
takes a vector $cond$ as additional input, which represents the condition
($D_{i^{*}}=k^{*}$) for some category $k^{*}$ of the $i^{*}$th categorical
feature $D_{i^{*}}\in\\{D_{1},\ldots,D_{N_{d}}\\}$ (∗ denotes the selected
feature/category in $cond$). It is worth noting that each condition only
involves _one_ category of _one_ categorical feature and does not involve any
continuous features. The sampling of $cond$ is twofold: sampling for the
$i^{*}$th categorical feature and sampling for the $k^{*}$th category of that
feature. Let $\oplus$ denote vector concatenation e.g., given
$\mathbf{x}_{1}=[0,0,0]$ and $\mathbf{x}_{2}=[1,0]$,
$\mathbf{x}_{1}\oplus\mathbf{x}_{2}=[0,0,0,1,0]$. If each categorical feature
$D_{i}$ is encoded as a one-hot vector
$\mathbf{d}_{i}=[\mathbf{d}_{i}^{(1)},\ldots,\mathbf{d}_{i}^{(k)},\ldots,\mathbf{d}_{i}^{(|D_{i}|)}]$,
where $|D_{i}|$ is the number of categories in $D_{i}$ and
$\mathbf{d}_{i}^{(k)}\in\\{0,1\\}$, then $cond$ is defined as
$cond=\mathbf{m}_{1}\oplus\ldots\oplus\mathbf{m}_{N_{d}}$ where
$\mathbf{m}_{i}=[\mathbf{m}_{i}^{(1)},\ldots,\mathbf{m}_{i}^{(k)},\ldots,\mathbf{m}_{i}^{(|D_{i}|)}]$
is the _mask_ vector of zeros associated with $\mathbf{d}_{i}$ with
$\mathbf{m}_{i}^{(k)}=1$ at $i=i^{*}$ and $k=k^{*}$. To ensure the generator
produces samples in accordance with the given conditions, a cross-entropy term
measuring the difference between $\mathbf{m}_{i^{*}}$ from the input $cond$
and the generated (denoted by ^) one-hot feature $\mathbf{\hat{d}}_{i^{*}}$ is
added to its loss.

To help the model evenly explore all possible categories in categorical
features, a procedure for sampling the $cond$ vector, termed _training-by-
sampling_ , is employed in CTGAN as follows: randomly choose a categorical
feature $D_{i^{*}}$ with uniform probability; construct the probability mass
function across the categories available in $D_{i^{*}}$ by taking the
logarithm of their frequencies in that feature (with respect to all training
examples in $\mathbf{T}_{train}$); then sample a category $k^{*}$ accordingly
and calculate the corresponding $\mathbf{m}_{i^{*}}$ and $cond$. Afterward,
the $cond$ vector is used to condition both synthetic and real training
samples in order for the discriminator to properly estimate the (Wasserstein)
distance between the learned and real conditional distributions
$P_{\mathcal{G}}(\text{row}|cond)$ and $P(\text{row}|cond)$, respectively.

We design our GAN to place more emphasis on generating samples conditioned on
$y$ by leveraging the conditional generator (discussed in the following
paragraph) as well as other known techniques (discussed in subsequent
subsections).

##### Proper Conditional Generation in HDLSS Setting.

In each iteration of CTGAN, $cond$ is sampled twice during the respective
training of the discriminator $\mathcal{D}$ and the generator $\mathcal{G}$,
with sample size equal the specified batch size for both. For our HDLSS
setting with $|\mathbf{T}_{train}|<200$ and over 600 columns/available
conditions to consider, even when the batch size is maximally set to
$|\mathbf{T}_{train}|$, any condition ($D_{i^{*}}=k^{*}$) would either be
missing or inadequately sampled in the minibatch, and hence it would be
impossible for $\mathcal{G}$ to properly learn to produce samples that
preserve the input conditions. Therefore, we increase the sample size of
$cond$ and hence the number of synthetic samples to generate during the
training of $\mathcal{G}$ by a factor of $q>1$ with respect to the batch size
($q=1$ in CTGAN). Moreover, instead of randomly sampling $D_{i^{*}}$ while
constructing $cond$, we want $\mathcal{G}$ to sample certain features more
frequently than others in order to prioritize learning their conditional
distributions, particularly for the target variable $y$ since our main purpose
of generating synthetic samples is to help classification models improve their
prediction on $y$. Therefore, we dedicate a portion of $cond$ to conditions
solely for $y$. We discuss on how we sample the remaining $D_{i^{*}}$’s in
Section 3.4.

###  3.3 Auxiliary Classifier

Incorporating an auxiliary classifier $\mathcal{C}$ into GAN architecture has
been shown to improve conditional generation quality for both image and
tabular data Odena et al. (2017); Park et al. (2018). $\mathcal{C}$ takes
either a synthetic or a real sample having its label removed as input and aims
to predict that label. Its loss, termed _classification loss_ Park et al.
(2018), quantifies the discrepancy between the label of a real sample and the
label predicted by $\mathcal{C}$ for that same sample, which is formulated as
binary and categorical cross entropy for Problems (A) and (B), respectively.
The addition of $\mathcal{C}$ also introduces an extra loss term into the loss
function of $\mathcal{G}$, which has the same form as the classification loss
but concerns synthetic samples instead. We refer to this loss as _downstream
loss_ as in Zhao et al. (2022). Altogether, $\mathcal{C}$ is trained to learn
the actual correlation between the label and the features, then teach
$\mathcal{G}$ how to generate realistic samples accordingly.

The synthetic samples that are fed to $\mathcal{C}$ are conditioned via the
$cond$ vector that is sampled earlier during the training of $\mathcal{G}$ (we
iteratively train
$\mathcal{D}\blacktriangleright\mathcal{G}\blacktriangleright\mathcal{C}$).
Since we are mainly concerned with learning the conditional distribution for
the target variable $P_{\mathcal{G}}(\text{row}|y)$, we only feed synthetic
samples that were conditioned on $y$ to $\mathcal{C}$ when computing the
downstream loss. The proper conditional generation step in Section 3.2 ensures
that we have sufficient amount of such samples for training $\mathcal{C}$.

###  3.4 Learning Important Features

The auxiliary classifier $\mathcal{C}$ was originally proposed to be an MLP
having the same architecture as $\mathcal{D}$ Park et al. (2018). On HDLSS
data, however, $\mathcal{C}$ is very likely to overfit, especially during the
first few iterations and epochs when $\mathcal{C}$ encounters few training
(either real or synthetic) examples, which in turn would negatively impact
$\mathcal{G}$. Recently, Margeloiu et al. (2023) propose a way to overcome
overfitting on HDLSS tabular classification problem by adding two auxiliary
networks before the first hidden layer of some classification neural network
in order to reduce the number of its learnable parameters and simultaneously
perform global feature selection. We generalize this idea by integrating these
networks into $\mathcal{C}$ to prevent it from overfitting. Furthermore, by
leveraging the global importance scores
$\bm{s}=[s_{1},\ldots,s_{N}]\in(0,1)^{N}$ (higher indicates greater
importance) for $N$ features, which are learned by one of the auxiliary
networks, we can inform the conditional generator of important categorical
features during training-by-sampling. More specifically, while sampling the
$cond$ vector, we sample features $D_{i}$ (excluding the target variable $y$)
each with probability proportional to its important score $s_{i}$. Hence,
conditions for features that have significant effects on the prediction of $y$
are sampled more frequently, allowing $\mathcal{G}$ to prioritize generating
synthetic samples conditioned on those features.

###  3.5 Enforcing Skip Logic

We first provide an intuitive explanation of how enforcing skip logic on
synthetic samples can benefit the training of our GAN. Recall that CTGAN
attempts to minimize the Wasserstein distance between
$P_{\mathcal{G}}(\text{row}|cond)$ and $P(\text{row}|cond)$, where $cond$
represents some condition ($D_{i^{*}}=k^{*}$) as defined earlier in Section
3.2. Let $cond_{1}$ and $cond_{2}$ be two distinct samples of $cond$ with
different sampled features e.g., $D_{1}$ ($i^{*}=1$) and $D_{2}$ ($i^{*}=2$),
respectively (and hence different sampled categories). With the presence of
skip logic, if $D_{1}$ and $D_{2}$ both belong to the same chain imposed by
another feature $D_{3}$ and the corresponding skip constraint
$D_{3}\rightarrow\\{D_{1},D_{2}\\}$ is enforced such that both $D_{1},D_{2}$
are omissible, then $cond_{1}$ and $cond_{2}$ must be redefined to match the
same conditional vector, $cond^{*}$, that satisfies such constraint i.e.,
($D_{1}=[\text{BLANK}]$) AND ($D_{2}=[\text{BLANK}]$). It follows that on real
data, $P(\text{row}|cond_{1})=P(\text{row}|cond_{2})=P(\text{row}|cond^{*})$.
This implies that the corresponding synthetic samples associated with either
condition should follow the same distribution
$P_{\mathcal{G}}(\text{row}|cond^{*})$. Therefore, enforcing skip logic by
inferring $cond^{*}$ effectively reduces the search space for
$P_{\mathcal{G}}$, which leads to more efficient and stable learning and hence
more consistency in the quality of the generated samples. We empirically
demonstrate this claim in Section 4.2.

Existing methods for enforcing column-wise constraints require either creating
customized transformation functions coupled with validity check or using
reject sampling Patki et al. (2016), both of which are ad hoc and highly
inefficient for large number of constrained columns. Instead, we leverage the
$cond$ vector to enforce our constraints. Formally, recall from Section 3.2
that $cond$ is constructed as $\bigoplus_{i=1}^{N_{d}}\mathbf{m}_{i}$ with
$\mathbf{m}_{i^{*}}^{(k^{*})}$ set to 1 and all other entries set to 0 for
representing the condition ($D_{i^{*}}=k^{*}$). Let us assume the
corresponding categorical feature $D_{i^{*}}$ from $cond$ imposes some skip
constraint $\kappa$ (e.g., $D_{i^{*}}=$ TB3 and $k^{*}=$ “No”) on a chain of
features $\\{D_{i^{\prime}}\\}_{i^{\prime}\in M}$, where
$M\subseteq\\{1,\ldots,N_{d}\\}$ such that $|M|$ is the size of such chain.
Let $k^{\prime}_{i^{\prime}}\in\\{1,\ldots,|D_{i^{\prime}}|\\}$ be a specific
category that each $D_{i^{\prime}}$ takes in accordance with $\kappa$ (e.g.,
$D_{i^{\prime}}=$ TB4 and $k^{\prime}_{i^{\prime}}=$ [BLANK] under
TB3$=$“No”). Then, we can define $\kappa$ as the extension of a condition as
follows:

$\displaystyle\begin{split}\kappa&=\left[D_{i^{*}}\rightarrow\\{D_{i^{\prime}}\\}_{i^{\prime}\in
M}\right]\\\ &=\left[(D_{i^{*}}=k^{*})\Rightarrow\bigwedge_{i^{\prime}\in
M}(D_{i^{\prime}}=k^{\prime}_{i^{\prime}})\right].\end{split}$ (1)

Therefore, whenever applicable, we can restrict
$cond\xrightarrow{\kappa}\zeta(cond)$ by reconstructing the individual mask
vectors of $cond$ (with a slight abuse of notation) as

$\mathbf{m}_{i}^{(k)}\xrightarrow{\kappa}\zeta\left(\mathbf{m}_{i}^{(k)}\right)=\begin{cases}1&\parbox[t]{65.04034pt}{if
$i\in\\{i^{*},i^{\prime}\\}$ and
$k\in\\{k^{*},k^{\prime}_{i^{\prime}}\\}$,}\\\ 0&\text{otherwise}.\end{cases}$
(2)

Note that $\kappa$ is only defined for a fixed set of
$\\{k^{*},k^{\prime}_{i^{\prime}}\\}$ and hence can take various forms in
practice. For instance, the following is an exhaustive list of valid
expressions for $\kappa=\text{TB3}\rightarrow\\{\text{TB4}\\}$:

  * $\blacktriangleright$

TB3$=$“No” $\Rightarrow$ TB4$=$[BLANK] (omissible)

  * $\blacktriangleright$

TB3$=$“Yes” $\Rightarrow$ TB4$=$“Not at all”

  * $\blacktriangleright$

TB3$=$“Yes” $\Rightarrow$ TB4$=$“Less than 1 cigarettes a day”

  * $\blacktriangleright$

TB3$=$“Yes” $\Rightarrow$ TB4$=$“1-5 cigarettes a day”

  * $\blacktriangleright$

TB3$=$“Yes” $\Rightarrow$ TB4$=$“Half a pack a day”

  * $\blacktriangleright$

TB3$=$“Yes” $\Rightarrow$ TB4$=$“A pack or more a day”.

###  3.6 The Complete Model

The complete training procedure via minibatch stochastic gradient descent
(SGD) is summarized in Algorithm 1. $q$ was previously defined in Section 3.2.
We denote the conditions for $y$ as $cond^{(y)}$. $\omega\in[0,1]$ is the
ratio for controlling the prevalence of such conditions in some sample of
$cond$. Let $\mathcal{I}$ be the probability mass function across the
categorical features (excluding $y$) wherein the probability for selecting a
feature $D_{i^{*}}$ is defined as its normalized feature importance score:
$s_{i^{*}}/\sum_{i=1}^{N_{d}-1}s_{i}$ ($-1$ for excluding $y$). The remaining
conditions in $cond$ are sampled following the training-by-sampling procedure
but with features selected according to $\mathcal{I}$ (the categories in each
categorical feature are still sampled according to their log probabilities in
$\mathbf{T}_{train}$ as before), for which we express as
$cond_{j}\sim\mathcal{I}$ with a slight abuse of notation. The loss for
$\mathcal{D}$, $L^{\mathcal{D}}$, is defined similarly as in CTGAN i.e., WGAN-
GP loss. The loss function for $\mathcal{G}$ is
$L^{\mathcal{G}}=L^{\mathcal{G}}_{orig}+L^{\mathcal{G}}_{dstream}$ where
$L^{\mathcal{G}}_{orig}$ is the original loss function for $\mathcal{G}$ in
CTGAN and $L^{\mathcal{G}}_{dstream}$ is the downstream loss defined in
Section 3.3 along with the classification loss $L^{\mathcal{C}}$.

At each iteration of an epoch, the training sequence is as follows: train the
discriminator (lines 3-6), train the generator (lines 7-10), and train the
auxiliary classifier while simultaneously updating $\mathcal{I}$ according to
$\bm{s}$ (lines 11-13).

Algorithm 1 Training Our Proposed GAN

1:Discriminator $\mathcal{D}$, generator $\mathcal{G}$, auxiliary classifier
$\mathcal{C}$, constants $q$ and $\omega$

2:Initialize $\mathcal{I}$ uniformly

3:for each minibatch $B$ of equal size from $\mathbf{T}_{train}$ do

4: Sample noise $Z$ and $cond$ (for training $\mathcal{D}$)
$Z=[z_{1},\ldots,z_{|B|}]$ with $z_{j}\sim\mathcal{N}(0,1)$
$cond=[cond_{1},\ldots,cond_{(1-\omega)*|B|}]$ with $cond_{j}\sim\mathcal{I}$

5: $cond\leftarrow\zeta(cond)\oplus[cond^{(y)}]_{\times\omega|B|}$

6: Sample $|B|$ rows from $\mathbf{T}_{train}$ such that each row
$\mathbf{r}_{j}$ satisfies $cond_{j}$

7: Perform SGD update of $\mathcal{D}$ with
$L^{\mathcal{D}}(\mathcal{G}(Z,cond),B)$

8: Sample noise $Z$ and $cond$ (for training $\mathcal{G}$ and $\mathcal{C}$)
$Z=[z_{1},\ldots,z_{q*|B|}]$ with $z_{j}\sim\mathcal{N}(0,1)$
$cond=[cond_{1},\ldots,cond_{q*(1-\omega)*|B|}]$ with
$cond_{j}\sim\mathcal{I}$

9: $cond\leftarrow\zeta(cond)\oplus[cond^{(y)}]_{\times q\omega|B|}$

10: Perform SGD update of $\mathcal{G}$ with
$L^{\mathcal{G}}_{orig}(\mathcal{G}(Z,cond))$

11: Perform SGD update of $\mathcal{G}$ with $\quad
L^{\mathcal{G}}_{dstream}\left(\mathcal{C}\left(\\{\mathcal{G}(Z,cond_{j})\\}_{cond_{j}\in
cond^{(y)}}\right)\right)$

12: Obtain $\bm{s}$ by computing $\mathcal{C}(B)$

13: Update $\mathcal{I}$ according to $\bm{s}$

14: Perform SGD update of $\mathcal{C}$ with $L^{\mathcal{C}}(\mathcal{C}(B))$

15:end for

16:Trained $\mathcal{G}$

##  4 Experiments

###  4.1 Methodology

####  4.1.1 Implementation Details

All experiments were conducted using PyTorch 1.13.1, CUDA 11.7, and scikit-
learn 1.3.2. Our implementation of the proposed
GAN111https://github.com/AnonyMouse3005/HDLSS-GAN is based on CTGAN’s. We
refer readers to our full paper for further implementation details,
particularly on our data preprocessing step.

##### Hyperparameters.

For all considered generative models, unless otherwise stated, we adopt the
same specifications as in the cited original work. We use a batch size of
$|B|=30$ to train each model, For CTGAN, we set the pac size Lin et al. (2018)
to 3. $\mathcal{C}$ is a 2-layer MLP with (256, 256) neurons and either a
sigmoid (for Problem A) or a softmax (for Problem B) activation in the last
layer. Each of the two attached auxiliary networks prior to the first hidden
layer of $\mathcal{C}$ is a 4-layer MLP (256, 256, 256, 256) with a tanh and a
sigmoid activation in the last layer, respectively. We manually tuned $q$ and
$\omega$ (via 10-fold cross validation on $\mathbf{T}_{train}$) before fixing
their values to 20 and 0.5, respectively. Each model is trained for 100
epochs, each contains $\lfloor|\mathbf{T}_{train}|/|B|\rfloor$ iterations. The
ratio for partitioning $\mathbf{T}_{real}$ is $80:20$ for $\mathbf{T}_{train}$
and $\mathbf{T}_{test}$, respectively, with a total of 100 distinct seeds.

####  4.1.2 Evaluation Metrics and Framework

Figure 2 illustrates our evaluation framework. We evaluate the efficacy of
generative models222For every measure of $\mathbf{T}_{syn}$’s quality while
evaluating some trained generator $G$, we use $G$ to generate 10 samples of
$\mathbf{T}_{syn}$ (each satisfies the two requirements specified in
Definition 2 (Tabular Data Generation).) and evaluate their quality
independently, then average the respective scores. using three criteria:
_conflict_ , _compatibility_ Park et al. (2018), and _utility_.

Figure 2: Workflow for evaluating the efficacy of generative models. The
“Classification models” blocks refer to the same set of classifiers listed in
the definition of Compatibility.. Each of these blocks takes a table as input
for training the classifiers and outputs their predictions on
$\mathbf{T}_{test}$. The blue dashed arrows represent the computation of the
respective scores/metrics.

##### Conflict.

Every row of a synthetic table generated by $\mathcal{G}$ should not contain
too many entries that violate skip logic. Given the $j$th row of
$\mathbf{T}_{syn}$ that is represented as
$\mathbf{\hat{r}}_{j}=\mathbf{\hat{C}}_{j}\oplus\mathbf{\hat{d}}_{1,j}\oplus\ldots\oplus\mathbf{\hat{d}}_{N_{d},j}$
where $\mathbf{\hat{d}}_{i,j}$ is the one-hot vector of the $i$th categorical
feature and $\mathbf{\hat{C}}_{j}$ is the representation of continuous
features in that row333varied across different generative models e.g., CTGAN
uses the proposed _mode-specific normalization_ (and so does our GAN), we
check for each skip constraint $\kappa$ whether the columns of
$\mathbf{\hat{r}}_{j}$ satisfy the condition ($D_{i^{*}}=k^{*}$) i.e., match
$cond$ (left-hand side of Equation 1). If it does, the one-hot vectors of the
features within the chain linked by $\kappa$ must exactly match the mask
vectors associated with said features in the restricted $cond$ vector,
$\zeta(\mathbf{m}_{i})$ (whose construction is defined in Equation 2). Hence,
we quantify the degree of $\kappa$-violation by computing the Hamming distance
between the vectors $(\mathbf{\hat{r}}_{i,j})_{i\in M}=\bigoplus_{i\in
M}\mathbf{\hat{d}}_{i,j}$ and $\bigoplus_{i\in M}\zeta(\mathbf{m}_{i})$, where
$M$ contains the indices for the features within the chain linked by $\kappa$.
The conflict metric of a single row $\mathbf{\hat{r}}_{j}$ is the average
Hamming distance across all applicable skip constraints, and we compute the
conflict of $\mathbf{T}_{syn}$ by taking its average across all rows. Thus, a
synthetic table whose rows adequately conform skip logic yields a low conflict
score in $[0,1]$.

##### Compatibility.

The classification models trained on synthetic data should output prediction
for unseen examples in test data as accurately as those trained on real
(training) data. We train classification models on $\mathbf{T}_{train}$ and on
$\mathbf{T}_{syn}$ (both having the same size as defined in Section 2.1), then
test them using $\mathbf{T}_{test}$ and compare their predictive performance,
which is measured using the standard Area under the Receiver Operating
Characteristic curve (AUROC). We train each of the multiclass classification
models in Problem (B) using the one-vs-all strategy and we compute the
resulting AUROC also using the one-vs-all strategy Fawcett (2006).

We report the average difference in AUROC (across different classification
models Xu et al. (2019) and different partitioning of $\mathbf{T}_{train}$ and
$\mathbf{T}_{test}$) of models trained on $\mathbf{T}_{syn}$ and those trained
on $\mathbf{T}_{train}$. It is expected for the classification models trained
on $\mathbf{T}_{syn}$ to score lower AUROC than those trained on
$\mathbf{T}_{train}$, ideally with a margin as small as possible. Therefore,
the compatibility score of $\mathbf{T}_{syn}$ should be negative and close to
zero. We consider the following classifiers same for both Problems (A) and
(B): elastic-net logistic regression Zou and Hastie (2005), decision tree
i.e., CART Breiman (2001), random forest Ho (1995), XGBoost Chen and Guestrin
(2016), CatBoost Prokhorenkova et al. (2018), 3-layer MLP (100, 100, 10) with
sigmoid/softmax activation in the last layer, and WPFS Margeloiu et al.
(2023).

##### Utility.

Ultimately, when the real training data is augmented by synthetic samples, the
classification models trained on it should excel compared to the models
trained on real data only. We report the average AUROC (across different
classification models and different partitioning of $\mathbf{T}_{train}$ and
$\mathbf{T}_{test}$) of the models trained on the augmented data,
$\mathbf{T}_{train}+\mathbf{T}_{syn}$, and compare it against the average
AUROC of those trained on $\mathbf{T}_{train}$. We consider the same set of
classifiers listed above. Unlike compatibility, we do not compute the
difference in AUROC scores. Hence, the utility of $\mathbf{T}_{syn}$ follows
the same scale as AUROC $\in[0,1]$ and should ideally yield values as high as
possible.

(a) Compatibility (greater i.e., less negative is better)

(b) Utility (greater is better)

Figure 3: Efficacy of considered generative models in Problem (A) (top) and
Problem (B) (bottom). Each column considers one drug. The red dashed lines in
Figure 3(b) mark the average AUROC of classification models trained on
$\mathbf{T}_{train}$ without data augmentation.

###  4.2 Results

We use our evaluation framework to evaluate CTGAN, TVAE Xu et al. (2019),
TableGAN Park et al. (2018), CTABGAN+ Zhao et al. (2022), and TabDDPM
Kotelnikov et al. (2023) in addition to the proposed GAN, which are considered
state-of-the-art for tabular data generation. For baseline, we use Bayesian
networks Pearl (2011).

##### Compatibility of Synthetic Data.

Figure 3(a) shows the performance of all considered generative models in terms
of compatibility. Across all considered drugs in both problems, we see that
the drop in predictive performance of classifications models trained on
synthetic data that are generated from our GAN ($<$0.05 and $<$0.07 lower in
AUROC for Problems A and B, respectively) is considerably smaller with lower
uncertainty compared to the performance drop of those trained on synthetic
data from other generative models. This demonstrates the efficacy of our GAN
in generating synthetic samples that are comparable to real training data.

##### Utility of Synthetic Data.

In terms of utility, the proposed GAN also performs well in each problem
relative to the benchmarks as shown in Figure 3(b). More specifically, the
classification models trained using augmented data from our GAN gain from
8.35% up to 13.4% in AUROC in Problem (A) and 8.66% up to 15.8% in Problem
(B), whereas those trained using augmented data from other models barely show
any improvement. This implies that our GAN is capable of augmenting existing
HDLSS data in order for classification models to effectively improve their
predictive performance.

##### Violation of Skip Logic in Synthetic Data.

The evaluation on degree of skip logic violation is summarized in Table 1.
Since our work is the first to take this criteria into account, our proposed
GAN has a clear advantage over existing models, with $\mathbf{T}_{syn}$
consistently having lower average conflict by a large margin in both
considered problems. The additional runtime for enforcing all 26 skip
constraints is negligible i.e., for approximately 5% longer training time.

Model | Problem (A) | Problem (B)
---|---|---
BN | $0.471\pm 0.113$ | $0.490\pm 0.154$
CTGAN | $0.339\pm 0.074$ | $0.358\pm 0.085$
TVAE | $0.297\pm 0.072$ | $0.325\pm 0.080$
TableGAN | $0.417\pm 0.108$ | $0.452\pm 0.094$
CTABGAN+ | $0.351\pm 0.086$ | $0.379\pm 0.084$
TabDDPM | $0.377\pm 0.107$ | $0.390\pm 0.090$
Our GAN | $\mathbf{0.196\pm 0.048}$ | $\mathbf{0.218\pm 0.039}$
|  |
Table 1: Conflict in $\mathbf{T}_{syn}$ (averaged over all splits of
$\mathbf{T}_{real}$, lower is better). Both problems concern meth usage.

##### Impact of Enforcing Skip Logic.

We also perform an ablation study to understand the practical benefits of
enforcing skip logic during our GAN training. As shown in Table 2, when we
enforce skip logic, the training phase exhibits not only higher stability but
also higher efficiency, as evidenced by the loss of $\mathcal{G}$ at 50 and
100 epochs. Note that the adoption of WGAN (via CTGAN) in our GAN allows us to
interpret the loss in a meaningful way. As a result, the trained GAN is able
to consistently generate high-quality samples, which positively affects the
scores for all three evaluation metrics to some extent. Similar results for
other considered problems can be found in the appendix of our full paper.

Criteria | Enforcing Skip Logic?
---|---
No | Yes
$L^{\mathcal{G}}_{orig}$ @50 ($\downarrow$) | $-2.276\pm 0.849$ | $\mathbf{-3.472\pm 0.653}$
$L^{\mathcal{G}}_{orig}$ @100 ($\downarrow$) | $-4.542\pm 0.933$ | $\mathbf{-4.871\pm 0.795}$
Conflict ($\downarrow$) | $0.358\pm 0.096$ | $\mathbf{0.196\pm 0.048}$
Compatibility ($\uparrow$) | $-0.026\pm 0.017$ | $\mathbf{-0.023\pm 0.009}$
Utility ($\uparrow$) | $0.818\pm 0.038$ | $\mathbf{0.821\pm 0.023}$
|  |
Table 2: Impact of enforcing skip logic during the training of our GAN in
Problem (A) for meth. The (original) loss of $\mathcal{G}$ is recorded at 50
epochs ($L^{\mathcal{G}}_{orig}$ @50) and 100 epochs ($L^{\mathcal{G}}_{orig}$
@100). ($\downarrow$) denotes less is better and ($\uparrow$) denotes greater
is better.

##  5 Conclusion

In this paper, using HDLSS tabular data collected by our team via a survey
that employs skip logic on short-term substance use behavior, we design a
novel GAN for augmenting our limited tabular data in order to help
classification models accurately predict short-term substance use behaviors of
PWUDs: (A) whether they would increase usage of a certain drug and (B) at
which ordinal frequency they would use it within the next 12 months. Our
evaluation results demonstrate the efficacy of the proposed GAN. The resulting
predictions for the two defined problems can ultimately be leveraged by
relevant substance use organizations as a complementary forecasting tool when
determining the most appropriate resource to allocate to PWUDs that need the
most help.

## Acknowledgements

This project is supported by the National Institute of General Medical
Sciences of the National Institutes of Health [P20GM130461], the Rural Drug
Addiction Research Center at the University of Nebraska-Lincoln, and the
National Science Foundation under grant IIS:RI #2302999. The content is solely
the responsibility of the authors and does not necessarily represent the
official views of the funding agencies.

## Ethics Statement

Due to the confidentiality agreement and IRB approval for this study, we do
not have access to sensitive information such as full name and gender
identity. Data will be made available to applicants (with an IRB protocol and
ethical research plan) upon request.

## References

  * Andersson et al. [2023] Helle Wessel Andersson, Mats P Mosti, and Trond Nordfjaern.  Inpatients in substance use treatment with co-occurring psychiatric disorders: a prospective cohort study of characteristics and relapse predictors.  BMC psychiatry, 23(1):1–10, 2023.
  * Arjovsky et al. [2017] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In International conference on machine learning, pages 214–223. PMLR, 2017.
  * Borisov et al. [2022] Vadim Borisov, Tobias Leemann, Kathrin Seßler, Johannes Haug, Martin Pawelczyk, and Gjergji Kasneci.  Deep neural networks and tabular data: A survey.  IEEE Transactions on Neural Networks and Learning Systems, 2022.
  * Breiman [2001] Leo Breiman.  Random forests.  Machine learning, 45:5–32, 2001.
  * Chawla et al. [2002] Nitesh V Chawla, Kevin W Bowyer, Lawrence O Hall, and W Philip Kegelmeyer.  Smote: synthetic minority over-sampling technique.  Journal of artificial intelligence research, 16:321–357, 2002.
  * Chen and Guestrin [2016] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, pages 785–794, 2016.
  * Colledge-Frisby et al. [2023] Samantha Colledge-Frisby, Sophie Ottaviano, Paige Webb, Jason Grebely, Alice Wheeler, Evan B Cunningham, Behzad Hajarizadeh, Janni Leung, Amy Peacock, Peter Vickerman, et al.  Global coverage of interventions to prevent and manage drug-related harms among people who inject drugs: a systematic review.  The Lancet Global Health, 2023.
  * Couper [2008] Mick P Couper.  Designing effective Web surveys. Cambridge University Press, 2008.
  * Dillman et al. [2014] Don A Dillman, Jolene D Smyth, and Leah Melani Christian.  Internet, phone, mail, and mixed-mode surveys: The tailored design method.  John Wiley & Sons, 2014.
  * Fawcett [2006] Tom Fawcett.  An introduction to roc analysis.  Pattern recognition letters, 27(8):861–874, 2006.
  * Fowler [1995] Floyd J Fowler.  Improving survey questions: Design and evaluation.  Sage, 1995.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  Advances in neural information processing systems, 27, 2014.
  * Gulrajani et al. [2017] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron C Courville.  Improved training of wasserstein gans.  Advances in neural information processing systems, 30, 2017.
  * Heckathorn [2014] Douglas D. Heckathorn.  Respondent-Driven Sampling: A New Approach to the Study of Hidden Populations*.  Social Problems, 44(2):174–199, 07 2014.
  * HHS [2023] HHS.  Sidebar: The Many Consequences of Alcohol and Drug Misuse — Surgeon General’s Report on Alcohol, Drugs, and Health [Internet].  https://addiction.surgeongeneral.gov/sidebar-many-consequences-alcohol-and-drug-misuse, 2023.  Cited: 2023-03-18.
  * Ho [1995] Tin Kam Ho.  Random decision forests.  In Proceedings of 3rd international conference on document analysis and recognition, volume 1, pages 278–282. IEEE, 1995.
  * Karamouzian et al. [2022] Mohammad Karamouzian, Andreas Pilarinos, Kanna Hayashi, Jane A Buxton, and Thomas Kerr.  Latent patterns of polysubstance use among people who use opioids: A systematic review.  International Journal of Drug Policy, 102:103584, 2022.
  * Kotelnikov et al. [2023] Akim Kotelnikov, Dmitry Baranchuk, Ivan Rubachev, and Artem Babenko.  Tabddpm: Modelling tabular data with diffusion models.  In International Conference on Machine Learning, pages 17564–17579. PMLR, 2023.
  * Lin et al. [2018] Zinan Lin, Ashish Khetan, Giulia Fanti, and Sewoong Oh.  Pacgan: The power of two samples in generative adversarial networks.  Advances in neural information processing systems, 31, 2018.
  * Linden-Carmichael et al. [2022] Ashley N Linden-Carmichael, Natalia Van Doren, Bethany C Bray, Kristina M Jackson, and Stephanie T Lanza.  Stress and affect as daily risk factors for substance use patterns: An application of latent class analysis for daily diary data.  Prevention science, 23(4):598–607, 2022.
  * Liu et al. [2022] Tennison Liu, Zhaozhi Qian, Jeroen Berrevoets, and Mihaela van der Schaar.  Goggle: Generative modelling for tabular data by learning relational structure.  In The Eleventh International Conference on Learning Representations, 2022.
  * Lorvick et al. [2023] Jennifer Lorvick, Jordana Hemberg, Madeleine J George, Joy Piontak, and Megan L Comfort.  Protocol: Understanding polysubstance use at the daily and event levels: protocol for a mixed-methods qualitative and ecological momentary assessment study in a community-based sample of people who use illicit drugs in oakland, california, usa.  BMJ Open, 13(9), 2023.
  * Margeloiu et al. [2023] Andrei Margeloiu, Nikola Simidjievski, Pietro Lio, and Mateja Jamnik.  Weight predictor network with feature selection for small sample tabular biomedical data.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, pages 9081–9089, 2023.
  * Newcomb and Bentler [1989] Michael D Newcomb and Peter M Bentler.  Substance use and abuse among children and teenagers.  American psychologist, 44(2):242, 1989.
  * NIDA [2020] NIDA.  Introduction [Internet].  https://nida.nih.gov/research-topics/commonly-used-drugs-charts, 2020.  Cited: 2023-03-18.
  * NIDA [2023] NIDA.  Drug Overdose Death Rates [Internet].  https://nida.nih.gov/research-topics/trends-statistics/overdose-death-rates, 7 2023.  Cited: 2023-07-11.
  * Odena et al. [2017] Augustus Odena, Christopher Olah, and Jonathon Shlens.  Conditional image synthesis with auxiliary classifier gans.  In International conference on machine learning, pages 2642–2651. PMLR, 2017.
  * ODPHP [2020] ODPHP.  Drug and alcohol use.  https://health.gov/healthypeople/objectives-and-data/browse-objectives/drug-and-alcohol-use, 2020.
  * Ouimette et al. [1997] Paige Crosby Ouimette, John W Finney, and Rudolf H Moos.  Twelve-step and cognitive-behavioral treatment for substance abuse: A comparison of treatment effectiveness.  Journal of consulting and clinical psychology, 65(2):230, 1997.
  * Park et al. [2018] Noseong Park, Mahmoud Mohammadi, Kshitij Gorde, Sushil Jajodia, Hongkyu Park, and Youngmin Kim.  Data synthesis based on generative adversarial networks.  arXiv preprint arXiv:1806.03384, 2018.
  * Patki et al. [2016] Neha Patki, Roy Wedge, and Kalyan Veeramachaneni.  The synthetic data vault.  In 2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA), pages 399–410. IEEE, 2016.
  * Pearl [2011] Judea Pearl.  Bayesian networks.  UCLA: Department of Statistics, 2011.
  * Prokhorenkova et al. [2018] Liudmila Prokhorenkova, Gleb Gusev, Aleksandr Vorobev, Anna Veronika Dorogush, and Andrey Gulin.  Catboost: unbiased boosting with categorical features.  Advances in neural information processing systems, 31, 2018.
  * Ray et al. [2020] Lara A Ray, Lindsay R Meredith, Brian D Kiluk, Justin Walthers, Kathleen M Carroll, and Molly Magill.  Combined pharmacotherapy and cognitive behavioral therapy for adults with alcohol or substance use disorders: a systematic review and meta-analysis.  JAMA network open, 3(6):e208279–e208279, 2020.
  * Russell et al. [2021] Cayley Russell, Farihah Ali, Frishta Nafeh, Sean LeBlanc, Sameer Imtiaz, Tara Elton-Marshall, and Jürgen Rehm.  A qualitative examination of substance use service needs among people who use drugs (pwud) with treatment and service experience in ontario, canada.  BMC Public Health, 21(1):2021, 2021.
  * SAMHSA [2022] SAMHSA.  Key Substance Use and Mental Health Indicators in the United States: Results from the 2021 National Survey on Drug Use and Health [Internet].  https://www.samhsa.gov/data/sites/default/files/reports/rpt39443/2021NSDUHFFRRev010323.pdf, 2022.  Cited: 2023-07-11.
  * Shorten and Khoshgoftaar [2019] Connor Shorten and Taghi M Khoshgoftaar.  A survey on image data augmentation for deep learning.  Journal of big data, 6(1):1–48, 2019.
  * Street et al. [1993] W Nick Street, William H Wolberg, and Olvi L Mangasarian.  Nuclear feature extraction for breast tumor diagnosis.  In Biomedical image processing and biomedical visualization, volume 1905, pages 861–870. SPIE, 1993.
  * Tanner-Smith et al. [2016] Emily E Tanner-Smith, Katarzyna T Steinka-Fry, Heather H Kettrey, and Mark W Lipsey.  Adolescent substance use treatment effectiveness: A systematic review and meta-analysis.  Peabody Research Institute, Vanderbilt University Nashville, TN, 2016.
  * SAMHSA [2021] SAMHSA.  Treatment episode data set: Admissions (teds-a).  https://www.datafiles.samhsa.gov/dataset/treatment-episode-data-set-admissions-2021-teds-2021-ds0001, 2021.
  * United Nations [2023a] United Nations.  Ensure healthy lives and promote well-being for all at all ages.  https://sdgs.un.org/goals/goal3, 2023.
  * United Nations [2023b] United Nations.  World drug report 2023.  https://www.unodc.org/unodc/en/data-and-analysis/wdr-2023-online-segment.html, 2023.
  * Wilson and Brown [2023] Anna M Wilson and Aaron R Brown.  Barriers to utilizing substance use disorder treatment and harm reduction services in appalachia.  Journal of Rural Mental Health, 2023.
  * Xu et al. [2019] Lei Xu, Maria Skoularidou, Alfredo Cuesta-Infante, and Kalyan Veeramachaneni.  Modeling tabular data using conditional gan.  Advances in neural information processing systems, 32, 2019.
  * Zhao et al. [2022] Zilong Zhao, Aditya Kunar, Robert Birke, and Lydia Y Chen.  Ctab-gan+: Enhancing tabular data synthesis.  arXiv preprint arXiv:2204.00401, 2022.
  * Zou and Hastie [2005] Hui Zou and Trevor Hastie.  Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 67(2):301–320, 2005.


HIN-18-008

HIN-18-008

# Two-particle azimuthal correlations in ${\PGg\Pp}$ interactions
using ${\Pp}\text{Pb}$ collisions at $\sqrtsNN=8.16\TeV$

###### Abstract

The first measurements of the Fourier coefficients ($V_{n\Delta}$) of the
azimuthal distributions of charged hadrons emitted from photon-proton
(${\PGg\Pp}$) interactions are presented. The data are extracted from 68.8of
ultra-peripheral proton-lead (${\Pp}\text{Pb}$) collisions at
$\sqrtsNN=8.16\TeV$ using the CMS detector. The high energy lead ions produce
a flux of photons that can interact with the oncoming proton. This ${\PGg\Pp}$
system provides a set of unique initial conditions with multiplicity lower
than in photon-lead collisions but comparable to recent electron-positron and
electron-proton data. The $V_{n\Delta}$ coefficients are presented in ranges
of event multiplicity and transverse momentum () and are compared to
corresponding hadronic minimum bias ${\Pp}\text{Pb}$ results. For a given
multiplicity range, the mean of charged particles is smaller in ${\PGg\Pp}$
than in ${\Pp}\text{Pb}$ collisions. For both the ${\PGg\Pp}$ and
${\Pp}\text{Pb}$ samples, ${V_{1\Delta}}$ is negative, ${V_{2\Delta}}$ is
positive, and ${V_{3\Delta}}$ consistent with 0. For each multiplicity and
range, ${V_{2\Delta}}$ is larger for ${\PGg\Pp}$ events. The ${\PGg\Pp}$ data
are consistent with model predictions that have no collective effects thus
suggesting the absence of collectivity in the ${\PGg\Pp}$ system over the
multiplicity range explored in this work.

##  0.1 Introduction

A wide variety of measurements suggest the existence of collectivity in the
collisions of small systems such as the proton-proton (${\Pp}{\Pp}$) [1, 2, 3,
4, 5] and proton-nucleus (pA) [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]
collisions. Such collectivity could indicate the formation of a hot, strongly
interacting “quark gluon plasma” (QGP), characterized by nearly ideal
hydrodynamic behavior [18, 19, 20], or could alternatively arise from gluon
saturation in the initial state [21, 22]. Properties of the QGP have been
previously studied in a wide range of high-energy nucleus-nucleus (AA)
collisions at the CERN LHC and BNL RHIC [23, 24, 25, 26, 27, 28, 29, 30, 31].
In these studies, collectivity is observed via the azimuthal correlations of
particles that are far apart in rapidity. This phenomenon is known as the
“ridge” [21], and has been observed in ${\Pp}{\Pp}$ and ${\Pp}\text{Pb}$
collisions that produce at least 110 charged particles in the final state [1,
2, 3, 4, 5, 7, 8, 6, 9, 10, 11, 12, 13, 14, 15, 16, 17]. The two-particle
azimuthal correlations can be characterized by their Fourier components
($V_{n\Delta}$) where $n$ represents the order of the moment. If the two-
particle correlations can be factorized into the product of the corresponding
single particle azimuthal distributions, then the single-particle azimuthal
anisotropy Fourier coefficients $v_{n}$ can be extracted as
$v_{n}=\sqrt{V_{n\Delta}}$ [32]. The second ($v_{2}$) and third ($v_{3}$)
coefficients are known as elliptic and triangular flow, respectively, and are
directly related to the initial collision geometry and its fluctuations, which
influence the medium evolution and and provide information about its
fundamental transport properties [33, 34, 35, 36].

In high-multiplicity events, $v_{2}$ and $v_{3}$ depend upon the hadron
species [15, 37, 38, 39, 40, 41] and scale with the number of valence quarks
in the hadron [15]. Such results suggest a common origin of the collectivity
seen in PbPb, as well as in high-multiplicity ${\Pp}{\Pp}$ and
${\Pp}\text{Pb}$ events, where a hydrodynamic description can be used to
reasonably reproduce the measurements in each case [42, 43, 44, 45]. Probing
systems with even smaller interaction regions is therefore important to
understand the reach of such a hydrodynamic description. The search for
collectivity has been recently extended to electron-positron, electron-proton,
and photon-nucleus interactions [46, 47, 48]. None of these systems exhibit
evidence of the long-range correlations seen in hadronic collisions.

High-energy ${\Pp}\text{Pb}$ ultra-peripheral collisions at the LHC, where the
impact parameter is larger than the nucleus radius provide a new system to
extend the search of long-range correlations to photon-proton (${\PGg\Pp}$)
collisions. At energies, the lead (Pb) nuclei generate a very large quasi-real
photon flux [49]. In the equivalent photon approximation [50, 51, 52], this
flux can be considered as $\PGg$ beams of virtuality $Q^{2}<1/R^{2}$, where
$R$ is the effective radius of the charge distribution. For Pb nuclei at
2.56with radius $R\approx 7\unit{fm}$, the quasi-real photon beams have
virtualities $Q^{2}<10^{-3}\GeV^{2}$, but very large longitudinal energy, up
to $E_{\PGg}=\hbar c/\alpha R\approx 80\GeV$, where $\alpha$ is the reciprocal
Lorentz relativistic factor.

This Letter presents a first search for long-range correlations and the first
measurement of $V_{n\Delta}$ coefficients for the ${\PGg\Pp}$ system. The CMS
detector has been used to collect a large sample of ${\PGg\Pp}$ interactions
that occur in ultra-peripheral ${\Pp}\text{Pb}$ collisions. The beam energies
were 6.50for the protons and 2.56per nucleon for the Pb nuclei, resulting in a
center-of-mass energy per nucleon pair ($\sqrtsNN$) of 8.16. The ${\PGg\Pp}$
results are compared to both hadronic minimum bias (MB) ${\Pp}\text{Pb}$
collisions (previously studied in Ref. [53]) and predictions of the v8.2 [54]
model interfaced with the Delphes v3.4.2 fast simulation package [55]. The
minimum bias data are compared to predictions from the v2.1 generator [56]
coupled to a full simulation of the detector [57].

##  0.2 Experimental apparatus and data sample

The central feature of the CMS apparatus is a superconducting solenoid of 6 m
internal diameter, providing a magnetic field of 3.8 T. Within the solenoid
volume is the silicon tracker, a lead tungstate crystal electromagnetic
calorimeter, and a brass and scintillator hadron calorimeter, each composed of
a barrel and two endcap sections that cover the range $\abs{\eta}<3.0$. The
silicon tracker measures charged particles within the range $\abs{\eta}<2.5$.
It consists of 1440 silicon pixel and 15 148 silicon strip detector modules,
and provides an impact parameter resolution of about $15\mum$ and a transverse
momentum () resolution better than 1.5% at $\pt\approx 100\GeVc$. Event
selection for this analysis makes use of detectors in the forward region:
hadron forward (HF) calorimeters that use quartz fibers embedded in a steel
absorber covering the region $3.0<\abs{\eta}<5.2$ and the two Zero Degree
Calorimeters (ZDCs) which measure neutral particles with $\abs{\eta}>8.3$
[58]. Analysis in the midrapidity region is based upon objects produced by the
CMS particle-flow (PF) algorithm [59], which reconstructs and identifies
final-state particles with an optimized combination of information from the
various elements of the CMS detector. A more detailed description of the CMS
detector, together with a definition of the coordinate system used and the
relevant kinematic variables, can be found in Ref. [60].

The analysis is performed using data recorded by CMS during the LHC
${\Pp}\text{Pb}$ run in 2016 with an integrated luminosity of 68.8 . The
proton-going direction is towards the side of the detector with positive
$\eta$. As a result of the energy difference between the colliding beams, the
nucleon-nucleon (NN) center-of-mass for ${\Pp}\text{Pb}$ collisions is not at
rest with respect to the laboratory frame. Massless particles emitted at
$\eta_{\text{cm}}=0$ in the NN center-of-mass frame will be detected at
$\eta=+0.465$ in the laboratory frame. The event samples were collected by the
CMS experiment with a two-level trigger system [53] consisting in the level-1
(L1), where events are selected by custom hardware processors and the high-
level trigger (HLT), that uses fast versions of the offline software.

##  0.3 Event selection

Figure 1: The $N_{\text{trk}}$ spectra for ${\PGg\Pp}$ and MB samples. The
simulated ${\PGg\Pp}$ distribution has been normalized to the same event yield
of the data ${\PGg\Pp}$ enhanced sample.

Samples of both ${\PGg\Pp}$-enhanced and MB events were collected requiring
energy deposits in at least one of the HF calorimeters above a threshold of
approximately 1at L1. The HLT system requires the presence of at least one
charged particle (track) with $\pt>0.4\GeVc$ in the pixel tracker. Track
reconstruction was performed online as part of the HLT trigger with a
reconstruction algorithm that is identical to the one used offline [61]. More
details of the MB trigger can be found in Ref. [62]. For each event, the
reconstructed vertex with the highest number of associated tracks was selected
as the primary vertex. A zero bias trigger requiring only the presence of
proton and lead bunches in the CMS detector was used to independently study
the trigger efficiency ($\varepsilon_{\text{trig}}$). The beam bunches were
detected by induction counters placed 175m from the interaction point on each
side of the experiment. In addition, a sample of events with neither beam
present was collected for noise studies.

For both ${\PGg\Pp}$ and MB samples, the reconstructed primary vertex was
required to be within 15cm of the nominal interaction point along the beam
axis ($z$) and within 0.15cm in the transverse plane. Only high-purity tracks,
described in Ref. [61], were used in this analysis. The impact parameter
significance of reconstructed tracks with respect to the primary vertex in the
longitudinal and transverse directions was required to be $<$3 standard
deviations. Finally, the relative uncertainty in the of the track was required
to be $<$10%. At least two reconstructed tracks with $\abs{\eta}<2.4$ and
$\pt>0.4\GeVc$ were required to be associated with the primary vertex. Beam-
related background was suppressed by rejecting events for which $<$25% of all
reconstructed tracks pass the standard track selection criteria.

Typical ${\Pp}\text{Pb}$ collisions produce particles at both positive and
negative rapidity [53, 63, 38]. However, ${\PGg\Pp}$ events are expected to be
very asymmetric in the laboratory frame since the photon energy is generally
much smaller than the proton beam energy.

For the ${\PGg\Pp}$-enhanced selection, a rapidity gap is defined as a
continuous region in which there is low detector activity, as done in Ref.
[64]. The detector acceptance $\abs{\eta}<5.0$ is divided into 20 bins.
Threshold values are assigned to each $\eta$ bin, they delimit the energy from
all PF candidates that can be considered significant and which contain at
least 99.7% of detector activity caused by detector noise or by beam-gas
events. These thresholds were obtained by studying the zero-bias events
triggered on non-colliding bunches. For each event, a given $\eta$ bin was
considered to be empty if the energy registered from the PF candidates was
below its assigned threshold value. For the 10 bins in the regions
$\abs{\eta}<2.5$ the energy threshold was 6and no high-purity tracks with
$\pt>200\MeVc$ were allowed. Activity in the forward lead-going region
$-5.0<\eta<-3.0$ was suppressed by implementing a set of threshold values that
range from 13.4 to 16.9. For the bin $-2.5>\eta>-3.0$ only neutral hadrons
were considered and the energy threshold was 13.4. The forward rapidity gap
(${\Delta\eta^{F}}$) variable was then defined as the difference from
$\eta=-5.0$ to the lower edge of the first non-empty $\eta$ bin.

The MB selection requires the coincidence of at least one tower with energy
above 3.0in both HF calorimeters and at least two tracks with
$\abs{\eta}<2.5$. In contrast, a ${\PGg\Pp}$-enhanced selection is designed to
capture events with an intact Pb nucleus, particle production in the positive
$\eta$ region, and a large rapidity gap [65, 66, 67]. The first two
requirements are met by requiring no neutrons in the ZDC on the Pb-going side
and at least 10in the highest energy tower of the HF calorimeter on the
p-going side. To ensure a large rapidity gap, we require
$5.0<{\Delta\eta^{F}}<7.5$. This corresponds to not having a particle within
the negative-$\eta$ region. A total of $8.6\times 10^{6}$ ${\PGg\Pp}$-enhanced
and $1.0\times 10^{9}$ MB candidate events were selected. From Ref. [64] the
purity of the ${\PGg\Pp}$-enhanced sample without the ZDC cut is estimated to
be about 95%. The requirement of no neutron emission used in this analysis
gives an additional suppression of pomeron-Pb events.

The reconstructed track multiplicity ($N_{\text{trk}}$) is defined as the
number of tracks from the primary vertex with $\pt>0.4\GeVc$, and
$\abs{\eta}<2.4$. Figure 1 shows the $N_{\text{trk}}$ spectra for the
${\PGg\Pp}$-enhanced and MB data samples along with simulations from the 8 and
event generators. For the ${\PGg\Pp}$-simulated sample, the events are
restricted to those with no tracks in the $\eta<0$ region and normalized to
the ${\PGg\Pp}$-enhanced yield. In contrast to the MB sample, the
${\PGg\Pp}$-enhanced spectrum drops very rapidly with multiplicity up to a
limiting value of 34. The $\langle N_{\text{trk}}\rangle$ value corresponding
to the $2\leq N_{\text{trk}}<35$ range for the ${\PGg\Pp}$-enhanced sample is
${\approx}2.9$ and about 16.6 for the MB sample. The $N_{\text{trk}}$
distribution from the zero bias data control sample has $\langle
N_{\text{trk}}\rangle\approx 0.84$. The ${\PGg\Pp}$-simulated sample shows a
shape and range that is consistent with the ${\PGg\Pp}$-enhanced data sample.
Three $N_{\text{trk}}$ bins are used to analyze the ${\PGg\Pp}$-enhanced
events: $2\leq N_{\text{trk}}<5$, $5\leq N_{\text{trk}}<10$, $10\leq
N_{\text{trk}}<35$. The first two deliver a comparable number of particle
pairs and the third one aims to probe the higher $N_{\text{trk}}$ domain by
averaging the last part of the distribution. Table 0.3 indicates the $\langle
N_{\text{trk}}\rangle$ values for the data and simulated ${\PGg\Pp}$ and MB
samples. The mean , $\langle p_{\mathrm{T}}\rangle$, values of charged
particles in the ${\PGg\Pp}$ and MB data samples are $0.67\pm 0.01$ and
$0.74\pm 0.01\GeVc$ respectively.

Mean $N_{\text{trk}}$ for the ${\PGg\Pp}$-enhanced and the MB data sets for
five classes of $N_{\text{trk}}$. Statistical uncertainties are negligible.
Sample $2\leq N_{\text{trk}}<5$ $5\leq N_{\text{trk}}<10$ $10\leq
N_{\text{trk}}<35$ $5\leq N_{\text{trk}}<35$ $2\leq N_{\text{trk}}<35$
${\PGg\Pp}$-enhanced 2.6 5.8 11.3 6.0 2.9 ${\PGg\Pp}$-simulated 2.6 5.9 11.4
6.2 2.9 MB 3.0 6.9 21.5 18.5 16.6 MB-simulated 3.1 6.9 20.7 17.2 15.7

##  0.4 Analysis technique

To ensure a high tracking efficiency, only tracks with $0.3<\pt<3.0\GeVc$ are
used in the analysis. The analysis techniques for two-particle correlations
are identical to those used for recent measurements [3, 6, 26] and are
described hereafter. For each multiplicity class, tracks tagged as “trigger
particles” are those whose , labeled as
${p\mathstrut}^{\text{\mathstrut{trig}}}_{\mathrm{T}}$, is within a given
range. The number of trigger particles in the event is denoted by
$N_{\text{trig}}$. Particle pairs are then formed by associating each trigger
particle with the remaining tracks within a specified
$\pt={p\mathstrut}^{\text{assoc}}_{\mathrm{T}}$ interval. For this analysis,
identical ranges are used for
${p\mathstrut}^{\text{\mathstrut{trig}}}_{\mathrm{T}}$ and
${p\mathstrut}^{\text{assoc}}_{\mathrm{T}}$. Two different ranges are studied,
i.e. $[0.3,3.0]$ and $[1.0,3.0]\GeVc$. These are the same as those used in
previous studies of the ridge [6] and observations of correlations between
$v_{n}$ coefficients [53] in ${\Pp}\text{Pb}$ collisions.

The two-dimensional correlation function is defined as

$\frac{1}{N_{\text{trig}}}\frac{\rd^{2}N^{\text{pair}}}{\rd{\Delta\eta}\,\rd{\Delta\phi}}=B(0,0)\frac{S({\Delta\eta},{\Delta\phi})}{B({\Delta\eta},{\Delta\phi})},$
(1)

Figure 2: Two-dimensional (left) and one-dimensional (right) correlation plots
for ${\PGg\Pp}$-enhanced (upper) and MB (lower) events for
$0.3<p_{\mathrm{T}}<3.0\GeVc$ and $2\leq N_{\text{trk}}<35$. For the two-
dimensional distributions, the jet peak centered at
${\Delta\eta}={\Delta\phi}=0$ is truncated to increase visibility. The
rapidity gap requirement for the ${\PGg\Pp}$-enhanced sample limits the
$\abs{\Delta\eta}$ range to $\abs{\Delta\eta}<2.5$. The one-dimensional
${\Delta\phi}$ distributions are symmetrized by construction around
${\Delta\phi}=0$ and $\pi$. The Fourier coefficients, $V_{n\Delta}$ in the
right column are fit over the ${\Delta\phi}$ range $[0,\pi]$. Points outside
this range are shown as open circles and are obtained by symmetrization of
those in $[0,\pi]$. Statistical error bars are shown for both one-dimensional
distributions.

where ${\Delta\eta}$ and ${\Delta\phi}$ are the differences in $\eta$ and
$\phi$ of the pair and $N^{\text{pair}}$ is the number of pairs. The same-
event pair distribution, $S({\Delta\eta},{\Delta\phi})$, represents the yield
of particle pairs from the same event in a given (${\Delta\eta}$,
${\Delta\phi}$) bin. The mixed-event pair distribution
$B({\Delta\phi},{\Delta\eta})$ is constructed by pairing the trigger particles
in each event with the associated charged particles from 100 different
randomly selected events in the same 0.5 wide vertex range and from the same
track multiplicity class. The same-event and mixed-event pair distributions
are first calculated for each event, and then averaged over the events within
the track multiplicity class. The mixed-event distribution is normalized by
the sum of background events. The ratio $B(0,0)/B({\Delta\eta},{\Delta\phi})$
accounts for pair acceptance effects, with $B(0,0)$ representing the mixed-
event associated yield for both particles of the pair going in the same
direction and thus having maximum pair acceptance.

Figure 2 (left) shows the two-particle correlation functions for
${\PGg\Pp}$-enhanced (upper row) and MB (lower row) events within the
multiplicity range $2\leq N_{\text{trk}}<35$ as functions of ${\Delta\eta}$
and ${\Delta\phi}$. For the ${\PGg\Pp}$ distribution, the ${\Delta\eta}$ range
is limited to $\abs{{\Delta\eta}}<2.5$ by the ${\Delta\eta^{F}}$ selection and
the acceptance of the tracker. Both distributions show a large jet peak
centered at ${\Delta\eta}={\Delta\phi}=0$, as well as a broader distribution
from the recoiling jet centered at ${\Delta\eta}=0$ and ${\Delta\phi}=\pi$.
Neither distribution displays a “ridge”-like structure at
$\abs{\Delta\phi}\approx 0$ for $\abs{\Delta\eta}>2$. Figure 2 (right) shows
the projections of the two-dimensional correlation functions onto the
${\Delta\phi}$ axis for $\abs{{\Delta\eta}}>2$, away from the jet
fragmentation peak. These distributions are fitted over the ${\Delta\phi}$
range $[0,\pi]$ to a Fourier decomposition series $\propto
1+\sum_{n}2V_{n\Delta}\cos(n\Delta\phi)$, from where the measured
$V_{n\Delta}$ are extracted. Only the first three terms are included in the
fit, since additional terms have a negligible effect on its quality.

Figure 3: Dependence of $V_{n\Delta}$ on $N_{\text{trk}}$ for ${\PGg\Pp}$ and
MB events for two different ranges. Systematic uncertainties are shown by the
shaded bars in the two panels. The $2\leq N_{\text{trk}}<5$, $5\leq
N_{\text{trk}}<10$, $10\leq N_{\text{trk}}<35$ are used for the lower range
and $2\leq N_{\text{trk}}<5$ and $5\leq N_{\text{trk}}<35$ for the higher
range. The points are placed at the mean value of the corresponding
$N_{\text{trk}}$ range. Lines indicate the prediction for ${\PGg\Pp}$ events
from 8.

##  0.5 Systematic uncertainties

The systematic uncertainties of the experimental procedure are evaluated by
varying the analysis conditions and extracting new $V_{n\Delta}$ coefficients.
The following effects were considered:

  1. 1.

The systematic uncertainties associated with the choice of the
${\Delta\eta^{F}}$ range were estimated by repeating the analysis with
${\Delta\eta^{F}}\in[4.5,5.0)$, just below the range of the nominal analysis.
This alternative selection decreases the purity of the ${\PGg\Pp}$-enhanced
sample by at most 8% [64]. The size of this uncertainty thus estimated to be
7% for ${V_{1\Delta}}$ and 27% for ${V_{2\Delta}}$. For the MB data there is
no rapidity gap requirement and so no systematic uncertainty is assigned for
this effect.

  2. 2.

The effect of tracking inefficiency and misreconstructed track rate was
studied by varying the track quality requirements. The selection thresholds on
the significance of the transverse and longitudinal track impact parameter
were varied from 2 to 5 standard deviations. In addition, the relative
uncertainty is varied from 0.05 to 0.10. This translates into a 3.5%
uncertainty in ${V_{1\Delta}}$ for the $2\leq N_{\text{trk}}<5$ category.

  3. 3.

The sensitivity of the results to the primary vertex position along the beam
axis ($z_{\text{vtx}}$) was quantified by comparing events with different
$z_{\text{vtx}}$ locations from $-15$ to $+15\unit{cm}$. The magnitude of this
systematic effect goes up to 150% for ${V_{3\Delta}}$ with numerical
estimations of ${\pm}0.003$ for $5\leq N_{\text{trk}}<10$ and $10\leq
N_{\text{trk}}<35$ respectively, in the $0.3<p_{\mathrm{T}}<3.0\GeVc$
category, and up to ${\pm}0.013$ for $1.0<p_{\mathrm{T}}<3.0\GeVc$.

  4. 4.

The trigger efficiency depends upon $N_{\text{trk}}$. It decreases
substantially for $N_{\text{trk}}<10$, reaching 70% for $N_{\text{trk}}=2$. To
study this effect, a parallel data sample with weighted events as
(1/$\varepsilon_{\text{trig}}$) was produced. The full difference of the
$V_{n\Delta}$ with and without the correction was taken as the uncertainty.
This uncertainty is 2.3% for ${V_{1\Delta}}$ and 17% for ${V_{2\Delta}}$ for
the sample with $2\leq N_{\text{trk}}<5$.

The systematic uncertainties were added in quadrature. For the
${\PGg\Pp}$-enhanced sample with $N_{\text{trk}}<35$ the final uncertainties
in $V_{n\Delta}$ are 8.4 and 31% for $n=1$ and 2, respectively. For the
minimum bias sample the uncertainties for ${V_{2\Delta}}$ are 11% for $2\leq
N_{\text{trk}}<5$ and smaller than 2.6% for the rest of the $N_{\text{trk}}$
range. Since ${p\mathstrut}^{\text{\mathstrut{trig}}}_{\mathrm{T}}$ and
${p\mathstrut}^{\text{assoc}}_{\mathrm{T}}$ have the same range, the
fractional uncertainties in $v_{n}$ are half those of $V_{n\Delta}$.

The $V_{n\Delta}$ coefficients for ${\PGg\Pp}$-enhanced events, as functions
of and $N_{\text{trk}}$. Statistical and systematic uncertainties are added in
quadrature.  range $2\leq N_{\text{trk}}<5$ $5\leq N_{\text{trk}}<10$ $10\leq
N_{\text{trk}}<35$ ${V_{1\Delta}}$ $-0.086\pm 0.006$ $-0.075\pm 0.005$
$-0.074\pm 0.007$ $0.3<p_{\mathrm{T}}<3.0\GeVc$ ${V_{2\Delta}}$ $0.012\pm
0.004$ $0.015\pm 0.004$ $0.026\pm 0.006$ ${V_{3\Delta}}$ $-0.002\pm 0.001$
$-0.002\pm 0.004$ $-0.010\pm 0.006$ $2\leq N_{\text{trk}}<5$ $5\leq
N_{\text{trk}}<35$ ${V_{1\Delta}}$ $-0.271\pm 0.021$ $-0.221\pm 0.017$
$1.0<p_{\mathrm{T}}<3.0\GeVc$ ${V_{2\Delta}}$ $0.077\pm 0.027$ $0.059\pm
0.017$ ${V_{3\Delta}}$ $-0.015\pm 0.009$ $-0.007\pm 0.013$

##  0.6 Results

In order to reduce the contribution to $v_{n}$ coefficients from back-to-back
jet correlations, one can correct $v_{n}$ by subtracting correlations from
very low-multiplicity events ($v^{\mathrm{sub}}_{n}$), as done in Refs. [4,
68, 53]. Because of the limited $N_{\text{trk}}$ range for the ${\PGg\Pp}$
events (up to ${\sim}35$), no low multiplicity subtraction is applied to these
data. Figure 3 and Table 0.5 show the measured $V_{n\Delta}$ coefficients as a
function of $N_{\text{trk}}$ for the two different $p_{\mathrm{T}}$ ranges for
the ${\PGg\Pp}$ and MB ${\Pp}\text{Pb}$ samples. For the MB sample, the
results are consistent with those in [53] before the subtraction procedure.
Both the ${\PGg\Pp}$ and MB distributions show a negative ${V_{1\Delta}}$, a
positive ${V_{2\Delta}}$ of smaller magnitude than ${V_{1\Delta}}$, and a
${V_{3\Delta}}$ that is consistent with zero. For a given $N_{\text{trk}}$ and
range, both ${V_{1\Delta}}$ and ${V_{2\Delta}}$ are larger in the ${\PGg\Pp}$
samples than in the MB results. For both samples, the magnitude of
${V_{1\Delta}}$ tends to decrease with $N_{\text{trk}}$, while ${V_{2\Delta}}$
has at most a weak $N_{\text{trk}}$ dependence. Their magnitude are both
larger in the higher range.

Figure 4:  Single-particle azimuthal anisotropy $v_{2}$ versus
$N_{\text{trk}}$ for ${\PGg\Pp}$-enhanced and ${\Pp}\text{Pb}$ samples in two
regions. Systematic uncertainties are shown by the shaded bars in the two
panels. Predictions from the 8 and generators are shown for the ${\PGg\Pp}$
and MB ${\Pp}\text{Pb}$ samples respectively. For the ${\PGg\Pp}$ events, same
$N_{\text{trk}}$ bin arrangement as in Figure 3 is kept while for
${\Pp}\text{Pb}$ the bins $[2,5)$, $[5,10)$, $[10,15)$ and $[15,20)$ are used.

Figure 3 also shows predictions from the 8 generator for $V_{n\Delta}$ from
${\PGg\Pp}$ collisions. The predictions of ${V_{2\Delta}}$ and ${V_{3\Delta}}$
from 8 are reasonably consistent with the ${\PGg\Pp}$ data and have a similar
dependence upon $p_{\mathrm{T}}$ and $N_{\text{trk}}$. The ${V_{1\Delta}}$
prediction is smaller in magnitude than the measured values for the low
$p_{\mathrm{T}}$.

Figure 4 shows $v_{2}$ as a function of $N_{\text{trk}}$ and for both
${\PGg\Pp}$ and MB data sets. For $0.3<p_{\mathrm{T}}<3.0\GeVc$, the MB
results are consistent with previously published CMS results [53]. Predictions
from the 8 and generators are also shown for ${\PGg\Pp}$ and MB
${\Pp}\text{Pb}$ interactions respectively, none of the models include
collective effects. For both data and simulations, $v_{2}$ varies slowly with
track multiplicity for the ${\PGg\Pp}$ and ${\Pp}\text{Pb}$ samples. At a
given $N_{\text{trk}}$, $v_{2}$ is larger in the higher range. This is similar
to trends observed in ep collisions [47]. The increase of $v_{2}$ with is also
present in the simulations although both generators slightly overshoot the
data at higher . For ${\Pp}\text{Pb}$ collisions it has been shown that
fluctuations in the proton shape can increase $v_{2}$ [69]. It is noticeable
that at a given and $N_{\text{trk}}$, $v_{2}$ is higher for ${\PGg\Pp}$ than
for ${\Pp}\text{Pb}$ interactions. Tabulated results are provided in the
HEPData record for this analysis [70].

##  0.7 Summary

For the first time, the study of long-range particle correlations has been
extended to photon-proton (${\PGg\Pp}$) interactions. This study used proton-
lead (${\Pp}\text{Pb}$) collisions at $\sqrtsNN=8.16\TeV$ recorded with the
CMS detector. The two-particle $V_{n\Delta}$ Fourier coefficients and
corresponding single-particle $v_{2}$ azimuthal anisotropies are reported as
functions of the multiplicity of charged hadrons ($N_{\text{trk}}$) for two
transverse momentum () ranges. For the ${\PGg\Pp}$ sample, the largest
observed multiplicity was $N_{\text{trk}}\sim 35$. The mean of charged
particles is smaller in the ${\PGg\Pp}$ sample than for ${\Pp}\text{Pb}$
collisions within the same multiplicity range. No evidence for a long-range
near-side ridge-like structure was found for either the ${\PGg\Pp}$ or
hadronic minimum bias ${\Pp}\text{Pb}$ (MB) samples within this multiplicity
range. In all $N_{\text{trk}}$ and ranges, ${V_{1\Delta}}$ is negative,
${V_{2\Delta}}$ is positive with a smaller magnitude than ${V_{1\Delta}}$, and
${V_{3\Delta}}$ is consistent with zero. The magnitudes of both
${V_{1\Delta}}$ and ${V_{2\Delta}}$ increase with . This increase has also
been seen in electron-proton collisions. At a given and track multiplicity,
$v_{2}$ is larger for ${\PGg\Pp}$-enhanced events than for MB ${\Pp}\text{Pb}$
interactions. Predictions from the 8 and models describe well the ${\PGg\Pp}$
and ${\Pp}\text{Pb}$ MB data at low , but slightly overestimate the data at
higher . Since these models do not have collective effects, these data suggest
the absence of collectivity in the ${\PGg\Pp}$ system over the multiplicity
range explored in this work.

## References

  * [1] CMS Collaboration, “Observation of long-range, near-side angular correlations in proton-proton collisions at the LHC”,  JHEP 09 (2010) 091, 10.1007/JHEP09(2010)091, arXiv:1009.4122.
  * [2] ATLAS Collaboration, “Observation of long-range elliptic azimuthal anisotropies in $\sqrt{s}$ $=$ 13 and 2.76 pp collisions with the ATLAS detector”,  Phys. Rev. Lett. 116 (2016) 172301, 10.1103/PhysRevLett.116.172301, arXiv:1509.04776.
  * [3] CMS Collaboration, “Measurement of long-range near-side two-particle angular correlations in pp collisions at $\sqrt{s}=13\TeV$”,  Phys. Rev. Lett. 116 (2015) 172302, 10.1103/PhysRevLett.116.172302, arXiv:1510.03068.
  * [4] CMS Collaboration, “Evidence for collectivity in pp collisions at the LHC”,  Phys. Lett. B 765 (2017) 193, 10.1016/j.physletb.2016.12.009, arXiv:1606.06198.
  * [5] ATLAS Collaboration, “Measurement of azimuthal anisotropy of muons from charm and bottom hadrons in pp collisions at $\sqrt{s}=13\TeV$ with the ATLAS detector”,  Phys. Rev. Lett. 124 (2020) 082301, 10.1103/PhysRevLett.124.082301, arXiv:1909.01650.
  * [6] CMS Collaboration, “Observation of long-range near-side angular correlations in pPb collisions at the LHC”,  Phys. Lett. B 718 (2013) 795, 10.1016/j.physletb.2012.11.025, arXiv:1210.5482.
  * [7] PHENIX Collaboration, “Measurement of long-range angular correlations and azimuthal anisotropies in high-multiplicity p$+$Au collisions at $\sqrtsNN=200$ GeV”,  Phys. Rev. C 95 (2017) 034910, 10.1103/PhysRevC.95.034910, arXiv:1609.02894.
  * [8] PHENIX Collaboration, “Creation of quark-gluon plasma droplets with three distinct geometries”,  Nature Phys. 15 (2019) 214, 10.1038/s41567-018-0360-0, arXiv:1805.02973.
  * [9] ATLAS Collaboration, “Observation of associated near-side and away-side long-range correlations in $\sqrtsNN=5.02\TeV$ proton-lead collisions with the ATLAS detector”,  Phys. Rev. Lett. 110 (2013) 182302, 10.1103/PhysRevLett.110.182302, arXiv:1212.5198.
  * [10] ATLAS Collaboration, “Measurement with the ATLAS detector of multi-particle azimuthal correlations in p$+$Pb collisions at $\sqrtsNN=5.02\TeV$”,  Phys. Lett. B 725 (2013) 60, 10.1016/j.physletb.2013.06.057, arXiv:1303.2084.
  * [11] ALICE Collaboration, “Long-range angular correlations on the near and away side in p-Pb collisions at $\sqrtsNN=5.02\TeV$”,  Phys. Lett. B 719 (2013) 29, 10.1016/j.physletb.2013.01.012, arXiv:1212.2001.
  * [12] LHCb Collaboration, “Measurements of long-range near-side angular correlations in $\sqrtsNN=5\TeV$ proton-lead collisions in the forward region”,  Phys. Lett. B 762 (2016) 473, 10.1016/j.physletb.2016.09.064, arXiv:1512.00439.
  * [13] ALICE Collaboration, “Long-range angular correlations of $\pi$, and in pPb collisions at $\sqrtsNN=5.02\TeV$”,  Phys. Lett. B 726 (2013) 164, 10.1016/j.physletb.2013.08.024, arXiv:1307.3237.
  * [14] CMS Collaboration, “Evidence for collective multiparticle correlations in p-Pb collisions”,  Phys. Rev. Lett. 115 (2015) 012301, 10.1103/PhysRevLett.115.012301, arXiv:1502.05382.
  * [15] CMS Collaboration, “Long-range two-particle correlations of strange hadrons with charged particles in pPb and PbPb collisions at LHC energies”,  Phys. Lett. B 742 (2015) 200, 10.1016/j.physletb.2015.01.034, arXiv:1409.3392.
  * [16] ATLAS Collaboration, “Measurement of multi-particle azimuthal correlations in pp, p$+$Pb and low multiplicity Pb$+$Pb collisions with the ATLAS detector”,  Eur. Phys. J. C. 77 (2017) 428, 10.1140/epjc/s10052-017-4988-1, arXiv:1705.04176.
  * [17] ATLAS Collaboration, “Measurement of long-range multiparticle azimuthal correlations with the subevent cumulant method in pp and p$+$Pb collisions with the ATLAS detector at the CERN large hadron collider”,  Phys. Rev. C 97 (2018) 024904, 10.1103/PhysRevC.97.024904, arXiv:1708.03559.
  * [18] J.-Y. Ollitrault, “Anisotropy as a signature of transverse collective flow”,  Phys. Rev. D 46 (1992) 229, 10.1103/PhysRevD.46.229.
  * [19] U. Heinz and R. Snellings, “Collective flow and viscosity in relativistic heavy-ion collisions”,  Ann. Rev. Nucl. Part. Sci. 63 (2013) 123, 10.1146/annurev-nucl-102212-170540, arXiv:1301.2826.
  * [20] C. Gale, S. Jeon, and B. Schenke, “Hydrodynamic modeling of heavy-ion collisions”,  Int. J. Mod. Phys. A 28 (2013) 1340011, 10.1142/S0217751X13400113, arXiv:1301.5893.
  * [21] K. Dusling, W. Li, and B. Schenke, “Novel collective phenomena in high-energy proton–proton and proton–nucleus collisions”,  Int. J. Mod. Phys. E 25 (2016) 1630002, 10.1142/S0218301316300022, arXiv:1509.07939.
  * [22] J. L. Nagle and J. Orjuela Koop, “A quasiparticle transport explanation for collectivity in the smallest of collision systems ($\Pp+\Pp$ and $\Pe^{+}\Pe^{-}$)”,  Nucl. Phys. A 982 (2019) 455, 10.1016/j.nuclphysa.2018.09.005, arXiv:1807.04619.
  * [23] STAR Collaboration, “Elliptic flow in au$+$au collisions at $\sqrtsNN=130$ GeV”,  Phys. Rev. Lett. 86 (2001) 402, 10.1103/PhysRevLett.86.402, arXiv:nucl-ex/0009011.
  * [24] PHOBOS Collaboration, “System size dependence of cluster properties from two-particle angular correlations in Cu+Cu and Au+Au collisions at $\sqrtsNN=200$ GeV”,  Phys. Rev. C 81 (2010) 024904, 10.1103/PhysRevC.81.024904, arXiv:0812.1172.
  * [25] CMS Collaboration, “Long-range and short-range dihadron angular correlations in central PbPb collisions at a nucleon-nucleon center of mass energy of 2.76 ”,  JHEP 10 (2011) 076, 10.1007/JHEP07(2011)076, arXiv:1105.2438.
  * [26] CMS Collaboration, “Centrality dependence of dihadron correlations and azimuthal anisotropy harmonics in PbPb collisions at $\sqrtsNN=2.76\TeV$”,  Eur. Phys. J. C 72 (2012) 10052, 10.1140/epjc/s10052-012-2012-3, arXiv:1201.3158.
  * [27] STAR Collaboration, “Long range rapidity correlations and jet production in high energy nuclear collisions”,  Phys. Rev. C 80 (2009) 064912, 10.1103/PhysRevC.80.064912, arXiv:0909.0191.
  * [28] STAR Collaboration, “Three-particle coincidence of the long range pseudorapidity correlation in high energy nucleus-nucleus collisions”,  Phys. Rev. Lett. 105 (2010) 022301, 10.1103/PhysRevLett.105.022301, arXiv:0912.3977.
  * [29] PHENIX Collaboration, “Formation of dense partonic matter in relativistic nucleus-nucleus collisions at RHIC: Experimental evaluation by the PHENIX collaboration”,  Nucl. Phys. A 757 (2005) 184, 10.1016/j.nuclphysa.2005.03.086, arXiv:nucl-ex/0410003.
  * [30] STAR Collaboration, “Beam-energy dependence of identified two-particle angular correlations in $\sqrtsNN=7.7-200\GeV$ Au+Au collisions”,  Phys. Rev. C 101 (2020) 014916, 10.1103/PhysRevC.101.014916, arXiv:1906.09204.
  * [31] STAR Collaboration, “Correlation measurements between flow harmonics in Au+Au collisions at RHIC”,  Phys. Lett. B 783 (2018) 459, 10.1016/j.physletb.2018.05.076, arXiv:1803.03876.
  * [32] S. Voloshin and Y. Zhang, “Flow study in relativistic nuclear collisions by fourier expansion of azimuthal particle distributions”,  Phys. Rev. C. 70 (1996) 665, 10.1007/s002880050141, arXiv:hep-ph/9407282.
  * [33] B. H. Alver, C. Gombeaud, M. Luzum, and J.-Y. Ollitrault, “Triangular flow in hydrodynamics and transport theory”,  Phys. Rev. C 82 (2010) 034913, 10.1103/PhysRevC.82.034913, arXiv:1007.5469.
  * [34] C. G. B. Schenke, S. Jeon, “Elliptic and triangular flow in event-by-event $\PD=3+1$ viscous hydrodynamics”,  Phys. Rev. Lett 106 (2011) 042301, 10.1103/PhysRevLett.106.042301, arXiv:1009.3244.
  * [35] C. S. Z. Qiu and U. Heinz, “Hydrodynamic elliptic and triangular flow in Pb-Pb collisions at $\sqrtsNN=2.76\TeV$”,  Phys. Lett. B 707 (2012) 151, 10.1016/j.physletb.2011.12.041, arXiv:1110.3033.
  * [36] B. Alver and G. Roland, “Collision geometry fluctuations and triangular flow in heavy-ion collisions”,  Phys. Rev. C 81 (2010) 054905, 10.1103/PhysRevC.82.039903, arXiv:1003.0194. [Erratum: Phys. Rev. C 82 (2010) 039903].
  * [37] CMS Collaboration, “Elliptic flow of charm and strange hadrons in high-multiplicity pPb collisions at $\sqrtsNN=8.16\TeV$”,  Phys. Rev. Lett. 121 (2018) 082301, 10.1103/PhysRevLett.121.082301, arXiv:1804.09767.
  * [38] CMS Collaboration, “Centrality and pseudorapidity dependence of the transverse energy density in pPb collisions at $\sqrtsNN=5.02\TeV$”,  Phys. Rev. C 100 (2019) 024902, 10.1103/PhysRevC.100.024902, arXiv:1810.05745.
  * [39] ALICE Collaboration, “Multi-strange baryon production in pPb collisions at $\sqrtsNN=5.02\TeV$”,  Phys. Lett. B 758 (2016) 389, 10.1016/j.physletb.2016.05.027, arXiv:1512.07227.
  * [40] PHENIX Collaboration, “Elliptic flow of identified hadrons in Au+Au collisions at $\sqrtsNN=200$ GeV”,  Phys. Rev. Lett. 91 (2003) 182301, 10.1103/PhysRevLett.91.182301, arXiv:nucl-ex/0305013.
  * [41] STAR Collaboration, “Distributions of charged hadrons associated with high transverse momentum particles in pp and Au+Au collisions at $\sqrtsNN=200$ GeV”,  Phys. Rev. Lett. 95 (2005) 152301, 10.1103/PhysRevLett.95.152301, arXiv:nucl-ex/0501016.
  * [42] R. D. Weller and P. Romatschke, “One fluid to rule them all: viscous hydrodynamic description of event-by-event central p+p, p+Pb and Pb+Pb collisions at $\sqrt{s}=5.02\TeV$”,  Phys. Lett. B 774 (2017) 351, 10.1016/j.physletb.2017.09.077, arXiv:1701.07145.
  * [43] P. Bozek, “Collective flow in p-Pb and d-Pb collisions at energies”,  Phys. Rev. C 85 (2012) 014911, 10.1103/PhysRevC.85.014911, arXiv:1112.0915.
  * [44] P. Bozek and W. Broniowski, “Correlations from hydrodynamic flow in p-Pb collisions”,  Phys. Lett. B 718 (2013) 1557, 10.1016/j.physletb.2012.12.051, arXiv:1211.0845.
  * [45] P. Bozek and W. Broniowski, “Collective dynamics in high-energy proton-nucleus collisions”,  Phys. Rev. C 88 (2013), no. 1, 014903, 10.1103/PhysRevC.88.014903, arXiv:1304.3044.
  * [46] A. Badea et al., “Measurements of two-particle correlations in $\Pe^{+}\Pe^{-}$ collisions at 91 with ALEPH archived data”,  Phys. Rev. Lett. 123 (2019) 212002, 10.1103/PhysRevLett.123.212002, arXiv:1906.00489.
  * [47] ZEUS Collaboration, “Two-particle azimuthal correlations as a probe of collective behaviour in deep inelastic ep scattering at HERA”,  JHEP 04 (2020) 070, 10.1007/JHEP04(2020)070, arXiv:1912.07431.
  * [48] ATLAS Collaboration, “Two-particle azimuthal correlations in photonuclear ultraperipheral Pb$+$Pb collisions at 5.02 with ATLAS”,  Phys. Rev. C 104 (2021) 014903, 10.1103/PhysRevC.104.014903, arXiv:2101.10771.
  * [49] A. J. Baltz, “The physics of ultraperipheral collisions at the LHC”,  Phys. Rept. 458 (2008) 1, 10.1016/j.physrep.2007.12.001, arXiv:0706.3356.
  * [50] C. F. von Weizsäcker, “Radiation emitted in collisions of very fast electrons”,  Z. Phys. 88 (1934) 612, 10.1007/BF01333110.
  * [51] E. J. Williams, “Nature of the high-energy particles of penetrating radiation and status of ionization and radiation formulae”,  Phys. Rev. 45 (1934) 729, 10.1103/PhysRev.45.729.
  * [52] E. Fermi, “On the theory of collisions between atoms and electrically charged particles”,  Nuovo Cim. 2 (1925) 143, 10.1007/BF02961914, arXiv:hep-th/0205086.
  * [53] CMS Collaboration, “Observation of correlated azimuthal anisotropy fourier harmonics in pp and p$+$Pb collisions at the LHC”,  Phys. Rev. Lett. 120 (2018) 092301, 10.1103/PhysRevLett.120.092301, arXiv:1709.09189.
  * [54] I. Helenius and C. O. Rasmussen, “Hard diffraction in photoproduction with 8”,  Eur. Phys. J. C 79 (2019), no. 5, 413, 10.1140/epjc/s10052-019-6914-1, arXiv:1901.05261.
  * [55] DELPHES 3 Collaboration, “DELPHES 3, a modular framework for fast simulation of a generic collider experiment”,  JHEP 02 (2014) 057, 10.1007/JHEP02(2014)057, arXiv:1307.6346.
  * [56] X.-N. Wang and M. Gyulassy, “: A monte carlo model for multiple jet production in pp, pA and AA collisions”,  Phys. Rev. D 44 (1991) 3501, 10.1103/PhysRevD.44.3501.
  * [57] GEANT4 Collaboration, “—a simulation toolkit”,  Nucl. Instrum. Meth. A 506 (2003) 250, 10.1016/S0168-9002(03)01368-8.
  * [58] O. Surányi, “Performance of the CMS zero degree calorimeters in the 2016 pPb run”,  J. Phys. Conf. Ser. 1162 (2019) 012005, 10.1088/1742-6596/1162/1/012005.
  * [59] CMS Collaboration, “Particle-flow reconstruction and global event description with the CMS detector”,  JINST 12 (2017) P10003, 10.1088/1748-0221/12/10/P10003, arXiv:1706.04965.
  * [60] CMS Collaboration, “The CMS experiment at the CERN LHC”,  JINST 3 (2008) S08004, 10.1088/1748-0221/3/08/S08004.
  * [61] CMS Collaboration, “Description and performance of track and primary-vertex reconstruction with the CMS tracker”,  JINST 9 (2014) P10009, 10.1088/1748-0221/9/10/P10009, arXiv:1405.6569.
  * [62] CMS Collaboration, “The CMS trigger system”,  JINST 12 (2017) P01020, 10.1088/1748-0221/12/01/P01020, arXiv:1609.02366.
  * [63] CMS Collaboration, “Pseudorapidity distributions of charged hadrons in proton-lead collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}=$ 5.02 and 8.16 ”,  JHEP 01 (2018) 045, 10.1007/JHEP01(2018)045, arXiv:1710.09355.
  * [64] CMS Collaboration, “First measurement of the forward rapidity gap distribution in $\mathrm{p}\mathrm{Pb}$ collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}=8.16\TeV$”, Technical Report CMS-PAS-HIN-18-019, 2020\.
  * [65] ATLAS Collaboration, “Rapidity gap cross sections measured with the ATLAS detector in pp collisions at $\sqrt{s}=7\TeV$”,  Eur. Phys. J. C 72 (2012) 1926, 10.1140/epjc/s10052-012-1926-0, arXiv:1201.2808.
  * [66] CMS Collaboration, “Measurement of diffraction dissociation cross sections in pp collisions at $\sqrt{s}=7\TeV$”,  Phys. Rev. D 92 (2015) 012003, 10.1103/PhysRevD.92.012003, arXiv:1503.08689.
  * [67] J. D. Bjorken, “Rapidity gaps and jets as a new physics signature in very high-energy hadron hadron collisions”,  Phys. Rev. D 47 (1993) 101, 10.1103/PhysRevD.47.101.
  * [68] CMS Collaboration, “Multiplicity and transverse momentum dependence of two- and four-particle correlations in pPb and PbPb collisions”,  Phys. Lett. B 724 (2013) 213, 10.1016/j.physletb.2013.06.028, arXiv:1305.0609.
  * [69] H. M$\ddot{a}$ntysaari, B. Schenke, C. Shen, and P. Tribedy, “Imprints of fluctuating proton shapes on flow in proton-lead collisions at the LHC”,  Phys. Lett. B 772 (2017) 681, 10.1016/j.physletb.2017.07.038, arXiv:1705.03177.
  * [70] “HEPData record for this analysis”, 2022.  10.17182/hepdata.89877.