# Classification of the Cervical Vertebrae Maturation (CVM) stages Using the
Tripod Network

###### Abstract

We present a novel deep learning method for fully automated detection and
classification of the Cervical Vertebrae Maturation (CVM) stages. The deep
convolutional neural network consists of three parallel networks (TriPodNet)
independently trained with different initialization parameters. They also have
a built-in set of novel directional filters that highlight the Cervical Verte
edges in X-ray images. Outputs of the three parallel networks are combined
using a fully connected layer. 1018 cephalometric radiographs were labeled,
divided by gender, and classified according to the CVM stages. Resulting
images, using different training techniques and patches, were used to train
TripodNet together with a set of tunable directional edge enhancers. Data
augmentation is implemented to avoid overfitting. TripodNet achieves the
state-of-the-art accuracy of 81.18% in female patients and 75.32% in male
patients. The proposed TripodNet achieves a higher accuracy in our dataset
than the Swin Transformers and the previous network models that we
investigated for CVM stage estimation.

${}^{\ddagger}$${}^{\ddagger}$footnotetext: These authors contributed equally
to this work.

Index Terms—  Deep learning, cervical vertebrae maturation, tripod network,
vision transformers.

##  1 Introduction

The success of orthodontic/orthopedic treatment depends on optimal treatment
timing, especially in addressing craniofacial skeletal imbalances. The optimal
treatment timing relies on identifying craniofacial skeletal maturity stages.
Bone age assessment using radiographic analyses was reported to be more
accurate than chronological age in determining skeletal maturation, growth
rate, the peak period of growth, and the remaining growth potential [1, 2, 3].
Cervical vertebra maturation (CVM) staging in lateral cephalometric
radiographs is a method to determine skeletal maturation. The validity and
reliability of the CVM staging have been supported by multiple studies [4, 5].
Cervical vertebrae are the first seven bones of the spinal column. Vertebral
growth involves changes in the size of vertebral bodies and the shape of the
upper and lower borders of C2, C3, C4 vertebrae. These changes have been
described into six stages, correlating with morphological modifications of the
vertebral shapes. The major limitation of the CVM method is that it is not
user-friendly and needs experienced practitioners; researchers reported poor
reproducibility among nonexpert examiners [6].

The use of machine learning (ML) techniques in the field of medical imaging is
rapidly evolving, and a fully automated diagnostic approach has gained
attention with its promise of reducing human error as well as the time and
effort needed for the task [7]. The application of Deep Learning (DL) to study
human growth and development from radiographs is a promising idea that needs
to be explored. The present study aims to apply a custom-designed DL method to
develop a fully automated machine-learning system to detect and classify the
CVM stages. There have been recent studies to use pre-trained networks to
create fully automated detection and classification of the CVM stages where
each utilizes a different dataset [8, 9]. In this study, we propose a custom-
designed network model to develop a fully automated system to detect and
classify the CVM stages. Our DL network has a tripod-like structure consisting
of three parallel networks which are independently trained with different
initialization parameters [10]. They also have a built-in set of novel
directional filters that highlight the edges of the cervical vertebrae in
X-ray images.

The TripodNet has some of the features of the transformer networks: (i) Input
images are divided into patches, (ii) input patches are augmented as in
transformers; and (iii) the network has a multiheaded structure. ResNet-20 is
used as the backbone of TripodNet. Output feature maps of the three parallel
ResNets are combined using a fully connected layer to produce the final
decision as shown in Fig. 1. Moreover, the age information of the patients is
also fed to the network to increase the accuracy of classification. The
dataset is divided by gender as male patients can have a different growth rate
than female patients. The proposed model achieves state-of-the-art performance
with 81.18% in female patients and 75.32% in male patients on the dataset
collection, which is superior to the previous results that used DL in a
traditional way including the straightforward implementation of the vision
transformer [11].

Fig. 1: The Model Diagram of the TriPodNet. Overlapping vertebrae images are
fed to each pod of the network.

##  2 Methodology

###  2.1 Overview of the CVM Dataset

The dataset used in developing our algorithm consists of digitized images of
scanned lateral cephalometric films for subjects aged between 4 and 29
obtained from the American Association of Orthodontists Foundation (AAOF)
Craniofacial Growth Legacy Collections, an open data source [12]. The images
were studied and labeled by the third author (MHE) who is an expert
orthodontist scientist with more than ten years of experience in classifying
CVM. Cervical maturation stages were classified into six stages (CS1- CS6).
Six stages are shown in Fig. 2 [2]. It is visible from Fig. 2 that the main
difference among the classes arises in the size and shape of C2, C3, and C4
vertebrae. The change may also happen in one vertebrae only, which may cause
confusion for a traditional CNN model which uses the entire image as its
input. For example, two classes CS1 and CS2 are only separated by the shape of
the lower bound of C2 vertebra whereas C3 and C4 vertebrae have the same
shape.

Fig. 2: Six CVM Stages and the corresponding X-ray images. The top drawing is
adapted from reference [2].

The main dataset was classified into six stages of CVM (CS1-CS6) and we will
examine the model performances on the six-stage classification problem. The
dataset consists of 1012 images classified by the third author (MHE) as the
principal evaluator. Out of 1012 images, 478 images belong to female patients,
and 534 images belong to male patients. In our entire dataset, the number of
lateral cephalograms belonging to cervical stages CS1, CS2, CS3, CS4, CS5, and
CS6 are 153, 182, 174, 159, 167, and 177, respectively. Data augmentation
methods such as random translation, AutoContrast [13], AugMix [14], and
RandAugment [15] are implemented with various combinations to avoid the
overfitting of the DL networks to the training dataset.

As it is stated earlier, the changes in C2, C3, and C4 vertebrae create 6
classes; therefore, their shape and size play a significant role in the
determination of CVM stages. TriPodNet requires three inputs to create an
output. In this study, for the network to distinguish the differences among
the stages, we used three patches per image where each patch containing one of
the C2, C3, or C4 vertebrae. Moreover, every patch is augmented at the
beginning of the TriPodNet to increase the performance of each model. Since
the size of the dataset is much smaller than the benchmark dataset used to
train ResNet-20, data augmentation and patch augmentation methods are
necessary to train the TriPodNet. Note that the patch augmentation is
different from the data augmentation method we used to avoid overfitting.
Rotation and grayscale jittering are implemented to create a different
instance of the same input each time an input is about to be fed into the
network.

###  2.2 Structure of the Multi-Pod Network (MultiPodNet)

A typical MultiPodNet utilizes two or more parallel Convolutional Neural
Networks (CNNs) performing the same computations and may process the input
image as sequential image patches in parallel as in transformer networks. The
original input image and its augmented versions are fed into convolutional
networks forming the MultiPod network and the output feature maps of parallel
convolutional networks are concatenated before the fully connected dense
layer. In this study, we use TriPodNet because it generated the best results
in our dataset and it is the best performer among the other models compared to
our arXiv manuscript [10]. Similar to [10], we also use ResNets as the
baseline network in this article. The TriPodNet model structure is shown in
Fig 1.

We first segment a given head and neck X-ray image and identify the spine
(cervical vertebrae) region using the so-called Aggregate Channel Features
(ACF) object detector [16] as shown in Fig. 1. This avoids the process of
manually cropping the spine region in each image in the database. As a result,
the skull, jaw, and irrelevant background regions are removed before the
images are applied to the deep learning algorithm. The ACF object detector
automatically extracts the Region of Interest (RoI) in the images thereby
reducing the search space of the deep learning structure. Because all the
segmented images have variable sizes, they are resized to a common size of
$77\times 35$.

After this step, we use image patches to create input image sequences for the
TriPodNet model. Shapes of C2, C3 and C4 vertebrae determine the stage of the
CVM. Therefore, we crop the segmented vertebrae image into $35\times 35$
patches that contain only one vertebra as shown in Fig. 1. Three patches are
derived from each image before the patch augmentation. The location of each
vertebra in the image is used to create the patches since C2, C3 and C4 are
always in similar positions in any image. We used overlapped windows to create
patches with the size of $35\times 35$. Before using the patches in model
training, we also benefit from patch augmentation where we rotate the images
randomly by 5 degrees. The rotation is necessary to help the model generalize
as the alignment of each vertebra depends on the posture of the patient.
Moreover, we apply grayscale jittering on the patches. With the grayscale
jittering, patches become brighter or darker, randomly. Patch augmentation is
visualized in Fig. 3. Similar to other image classification methods,
augmentation is only applied to the training images. We do not apply patch
augmentation on testing images.

Fig. 3: Patch augmentation on the training images.

Next, instead of feeding the image patches directly to the ResNets, the edges
of the vertebral body in the patches are emphasized using eight directional
filters described in [17, 18]. The outputs of directional filters are then fed
to the ResNets. The same directional filters in [11] are used before ResNets
to highlight the edges of CVM images. The motivation behind using the
directional filters is to start the deep learning model using our prior domain
knowledge. Since the angles of CVM bones are the distinguishing factor in
estimating the CVM stages highlighting the edges of the bones in multiple
directions will give an advantage to the deep network model.

As pointed out above each pod of the TriPodNet is a ResNet-20. The structure
of ResNet-20 with the directional filters and patch augmentation is summarized
in Table 1. The ResNets are implemented in parallel. The output feature maps
of ResNets are concatenated together with the age of the subject. We also use
the chronological age information as the input to the fully connected layer.
Obviously, the chronological age of a subject is correlated with the maturity
of a patient. The age information is repeated six times in the vector and
Gaussian noise with zero mean and 0.01 variance is added to secure its impact
before the output layer.

Layer | Output Shape | Implementation Details
---|---|---
PatchAug | $35\times 35\times 8$ | -
DirFilts | $35\times 35\times 8$ | $7\times 7,8$
Conv1 | $35\times 35\times 16$ | $3\times 3,16$
Conv2_x | $35\times 35\times 16$ | $\left[\begin{array}[]{c}3\times 3,16\\\ 3\times 3,16\end{array}\right]\times 3$
Conv3_x | $17\times 17\times 32$ | $\left[\begin{array}[]{c}3\times 3,32\\\ 3\times 3,32\end{array}\right]\times 3$
Conv4_x | $8\times 8\times 64$ | $\left[\begin{array}[]{c}3\times 3,64\\\ 3\times 3,64\end{array}\right]\times 3$
AAP | $64$ | Adaptive Average Pooling
Table 1: Structure of each pod (ResNet with the directional filters) of the
TriPodNet. PatchAug stands for patch augmentation. DirFilts stands for
directional filters [11].

We selected the TriPodNet with random translation and AutoContrast
augmentation methods [13] as the best network model in our dataset based on
our experiments. In the next section, we present our experimental results.

##  3 Experimental Results

We studied different networks, and different data augmentation methods to
determine the best possible network structure. As introduced in [10],
MultiPodNet can be constructed from two or more pod networks processing the
input in parallel. Outputs of pod networks can be combined by adding the
feature maps or by concatenation that produced a higher accuracy than adding
the feature maps of individual networks. We studied a single channel ResNet,
DuPodNet which consists of two pod networks, TriPodNet, and QuadPodNet which
consists of four parallel pod networks. The TriPodNet with input overlapping
input image patching strategy as shown in Fig. 1 produced the highest accuracy
in our dataset. We compare the TripodNet constructed from ResNet-20s [19] with
the Swin Transformer [20], Xception [21], MobileNet-V1 [22], MobileNet-V2
[23], and the custom designed CNN with the directional filters that achieved
the previous best result in [11]. In addition to random translation and
AutoContrast data augmentation methods we studied ”randAugment”, ”AugMix”
methods and without any augmentation.

To train the networks, an SGD optimizer with a weight decay of 0.0001 and
momentum of 0.9 is used. These models are trained with a batch size of 32, an
initial learning rate of 0.1 for 100 epochs, and the learning rate is reduced
to 1/10 at epochs 25, 50, and 75. The experiments are implemented using
PyTorch in Python 3.

The accuracy results of the TripodNet with and without directional filters and
with different augmentation methods are summarized in Table 2. Directional
filters improve TripodNet’s accuracy by 4.88% (from 76.29% to 81.17%) in the
dataset containing female subjects. Similarly, they improve the accuracy by
4.17% (from 71.15% to 75.32%) in the dataset containing male subject images.
As pointed out above we augment both the entire image and the patches. We
trained the network using random translations with AutoContrast as the data
augmentation method as shown in Table 3. Furthermore, we augment the input
image patches to make the system robust to changes in posture and exposure as
shown in Table 2 during training.

Directional | Augmentation | Accuracy
---|---|---
Filters | Data | Patch | Female | Male
$\times$ | $\times$ | $\times$ | 67.05% | 65.38%
$\times$ |  | $\times$ | 75.11% | 70.19%
|  | $\times$ | 78.64% | 70.19%
$\times$ |  |  | 76.29% | 71.15%
|  |  |  81.17% |  75.32%
Table 2: Accuracy results of the TriPodNet with and without directional filters, augmentation of entire input images, and augmentation of image patches. Data Augmentation Method | Female | Male
---|---|---
No augmentation | 68.23% | 66.34%
Random augmentation | 74.11% | 70.18%
AugMix | 77.64% | 69.23%
Random translation, AutoContrast | 81.17% | 75.32%
Table 3: Accuracy results of various data augmentation methods. Directional
filters and patch augmentation are applied. AugMix, AutoContrast are PyTorch
commands [13].

The accuracy of different MultiPod networks is presented in Table 4. TriPodNet
achieves the highest accuracy. QuadPodNet has more parameters compared to
TriPodNet but it produced a lower accuracy. This may be due to the small
dataset size and training issues. StackNet in Table 4 has also three parallel
pod networks. In the StackNet, we stack the input patches and feed the
augmented stack patches to the pod networks at the same time.

Model | Parameters | Female | Male
---|---|---|---
ResNet | 0.27M | 75.29% | 73.07%
DuPodNet | 0.54M | 75.29% | 68.27%
StackNet | 0.82M | 80.03% | 70.19%
TriPodNet | 0.81M | 81.17% | 75.32%
QuadPodNet | 1.08M | 76.84% | 71.43%
Table 4: Accuracy results of MultiPod networks. DuPodNet uses two ResNet pods,
TriPodNet uses three, and QuadPodNet uses four. StackNet feeds all three
patches stacked to three pods.

(a)  Female patients. (b) Male patients.

Fig. 4: Test accuracy versus the number of epochs of different networks. The
TripodNet reaches the highest accuracy on both genders.

Comparison of the TripodNet with other state-of-art networks is presented in
Table 5 and Fig. 4. Fig. 4 shows the test accuracy versus the number of epochs
of different networks. In Table 5, we compare the TripodNet with ResNet-20
[19], Swin Transformer [20], Xception [21], MobileNet-V1 [22], MobileNet-V2
[23], and the custom-designed CNN with the directional filters [11]. We adjust
the image sizes to pre-trained networks: we remove some of the downsampling
layers in the MobileNet-V1, MobileNet-V2, and Xception. We implement the Swin
Transformer using two stages, where each stage contains a patch merging layer
with 4 Swin Transformer blocks. Moreover, it is well-known that after pre-
training on a large dataset, a transformer can get a much higher accuracy
(this is known as transfer learning). Therefore, to compare with the Swin
Transformer, we also use the ImageNet-1k-pre-trained Swin-s and Swin-t
transformers. In these cases, the input tensors are interpolated to $224\times
224$ to adjust the CVM image sizes to transformer input sizes. Table 5 and
Fig. 4 show that our proposed TripodNet produces the highest accuracy in our
data set. The TripodNet gets better results than the pre-trained Swin-s and
Swin-t transformers. This may be due to the ”tiny” size of the CVM data set.
It is probably too small for the transformer networks. Another reason is that
transfer learning may not work properly because the CVM images are quite
different from the ImageNet images.

Model | Parameters | Female | Male
---|---|---|---
MobileNet-V1 | 3.21M | 73.20% | 60.40%
MobileNet-V2 | 2.29M | 72.16% | 68.32%
Swin Transformer | 4.12M | 63.58% | 62.50%
Pre-trained Swin-s | 48.8M | 75.26% | 73.27%
Pre-trained Swin-t | 28.5M | 72.16% | 74.26%
Directional CNN [11] | 0.71M | 70.33% | 71.88%
ResNet-20 | 0.27M | 75.29% | 73.07%
Xception | 33.0M | 78.82% | 71.15%
TripodNet | 0.81M | 81.17% | 75.32%
Table 5: Accuracy results of various networks. Swin-s and Swin-t were pre-
trained on ImageNet-1K.

Fig. 5 shows the confusion matrices of the TriPodNet in the male and female
patient datasets.

(a) Female patients. (b) Male patients.

Fig. 5: The confusion matrices of the TriPodNet.

##  4 Conclusion

In this paper, we present a new method for CVM classification. We introduce a
novel neural network which is a combination of three parallel networks for
this purpose. Similar to the transformer networks, TriPodNet performs its
computations in parallel using traditional CNNs and the output feature maps of
parallel CNNs are combined using a fully connected layer to produce the final
result. We also compared the results of two, three, four, or more parallel
networks. The TriPodNet with three parallel ResNet-20s produced the best
accuracy result.

The transformer networks did not produce as good results as the TriPodNet.
This is probably because our tiny data set is too small for the transformer
networks.

## References

  * [1] Sunjay Suri, Chandrakala Prasad, Bryan Tompson, and Wendy Lou,  “Longitudinal comparison of skeletal age determined by the greulich and pyle method and chronologic age in normally growing children, and clinical interpretations for orthodontics,”  American Journal of Orthodontics and Dentofacial Orthopedics, vol. 143, no. 1, pp. 50–60, 2013.
  * [2] Lorenzo Franchi, Tiziano Baccetti, Laura De Toffol, Antonella Polimeni, and Paola Cozza,  “Phases of the dentition for the assessment of skeletal maturity: a diagnostic performance study,”  American Journal of Orthodontics and Dentofacial Orthopedics, vol. 133, no. 3, pp. 395–400, 2008.
  * [3] Zachary J Mellion, Rolf G Behrents, and Lysle E Johnston Jr,  “The pattern of facial skeletal growth and its relationship to various common indexes of maturation,”  American Journal of Orthodontics and Dentofacial Orthopedics, vol. 143, no. 6, pp. 845–854, 2013.
  * [4] Leonard S Fishman,  “Radiographic evaluation of skeletal maturation: a clinically oriented method based on hand-wrist films,”  The Angle Orthodontist, vol. 52, no. 2, pp. 88–112, 1982.
  * [5] Brent Hassel and Allan G Farman,  “Skeletal maturation evaluation using cervical vertebrae,”  American Journal of Orthodontics and Dentofacial Orthopedics, vol. 107, no. 1, pp. 58–66, 1995.
  * [6] Trenton S Nestman, Steven D Marshall, Fang Qian, Nathan Holton, Robert G Franciscus, and Thomas E Southard,  “Cervical vertebrae maturation method morphologic criteria: poor reproducibility,”  American Journal of Orthodontics and Dentofacial Orthopedics, vol. 140, no. 2, pp. 182–188, 2011.
  * [7] June-Goo Lee, Sanghoon Jun, Young-Won Cho, Hyunna Lee, Guk Bae Kim, Joon Beom Seo, and Namkug Kim,  “Deep learning in medical imaging: general overview,”  Korean journal of radiology, vol. 18, no. 4, pp. 570–584, 2017\.
  * [8] Hyejun Seo, JaeJoon Hwang, Taesung Jeong, and Jonghyun Shin,  “Comparison of deep learning models for cervical vertebral maturation stage classification on lateral cephalometric radiographs,”  Journal of Clinical Medicine, vol. 10, no. 16, pp. 3591, 2021.
  * [9] Hossein Mohammad-Rahimi, Saeed Motamadian, Mohadeseh Nadimi, Sahel Hassanzadeh-Samani, Mohammad Minabi, Erfan Mahmoudinia, Victor Lee, and Mohammad Hossein Rohban,  “Deep learning for the classification of cervical maturation degree and pubertal growth spurts: A pilot study,”  Korean journal of orthodontics, vol. 52, pp. 112–122, 03 2022.
  * [10] Hongyi Pan, Salih Atici, and Ahmet Enis Cetin,  “Multipod convolutional network,”  arXiv preprint arXiv:2210.00689, 2022.
  * [11] Salih Furkan Atici, Rashid Ansari, Veerasathpurush Allareddy, Omar Suhaym, Ahmet Enis Cetin, and Mohammed H Elnagar,  “Fully automated determination of the cervical vertebrae maturation stages using deep learning with directional filters,”  Plos one, vol. 17, no. 7, pp. e0269198, 2022.
  * [12] “Aaof craniofacial growth legacy collection,” https://www.aaoflegacycollection.org/aaof_home.html Accessed: 2021-10-03.
  * [13] “Pytorch autocontrast,” https://pytorch.org/vision/stable/generated/torchvision.transforms.functional.autocontrast.html Accessed: 2022-10-03.
  * [14] Dan Hendrycks, Norman Mu, Ekin D Cubuk, Barret Zoph, Justin Gilmer, and Balaji Lakshminarayanan,  “Augmix: A simple data processing method to improve robustness and uncertainty,”  arXiv preprint arXiv:1912.02781, 2019.
  * [15] Ekin D Cubuk, Barret Zoph, Jonathon Shlens, and Quoc V Le,  “Randaugment: Practical automated data augmentation with a reduced search space,”  in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition workshops, 2020, pp. 702–703.
  * [16] Piotr Dollár, Ron Appel, Serge Belongie, and Pietro Perona,  “Fast feature pyramids for object detection,”  IEEE transactions on pattern analysis and machine intelligence, vol. 36, no. 8, pp. 1532–1545, 2014.
  * [17] Ahmet M Bagci, Rashid Ansari, and William D Reynolds,  “Low-complexity implementation of non-subsampled directional filter banks using polyphase representations and generalized separable processing,”  in 2007 IEEE International Conference on Electro/Information Technology. IEEE, 2007, pp. 422–427.
  * [18] Alican Bozkurt, Alexander Suhre, and A Enis Cetin,  “Multi-scale directional-filtering-based method for follicular lymphoma grading,”  Signal, Image and Video Processing, vol. 8, no. 1, pp. 63–70, 2014\.
  * [19] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [20] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo,  “Swin transformer: Hierarchical vision transformer using shifted windows,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2021, pp. 10012–10022.
  * [21] François Chollet,  “Xception: Deep learning with depthwise separable convolutions,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 1251–1258.
  * [22] Andrew G Howard, Menglong Zhu, Bo Chen, Dmitry Kalenichenko, Weijun Wang, Tobias Weyand, Marco Andreetto, and Hartwig Adam,  “Mobilenets: Efficient convolutional neural networks for mobile vision applications,”  arXiv preprint arXiv:1704.04861, 2017.
  * [23] Mark Sandler, Andrew Howard, Menglong Zhu, Andrey Zhmoginov, and Liang-Chieh Chen,  “Mobilenetv2: Inverted residuals and linear bottlenecks,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2018, pp. 4510–4520.


# Medium-Assisted Enhancement of Exotic Hadron Production
from Small to Large Colliding Systems

Yu Guo  Guangdong Provincial Key Laboratory of Nuclear Science, Institute of
Quantum Matter, South China Normal University, Guangzhou 510006, China
Guangdong-Hong Kong Joint Laboratory of Quantum Matter, Southern Nuclear
Science Computing Center, South China Normal University, Guangzhou 510006,
China  Xingyu Guo<EMAIL_ADDRESS>Guangdong Provincial Key Laboratory of
Nuclear Science, Institute of Quantum Matter, South China Normal University,
Guangzhou 510006, China  Guangdong-Hong Kong Joint Laboratory of Quantum
Matter, Southern Nuclear Science Computing Center, South China Normal
University, Guangzhou 510006, China  Jinfeng Liao<EMAIL_ADDRESS>Physics
Department and Center for Exploration of Energy and Matter, Indiana
University, 2401 N Milo B. Sampson Lane, Bloomington, IN 47408, USA  Enke Wang
<EMAIL_ADDRESS>Guangdong Provincial Key Laboratory of Nuclear Science,
Institute of Quantum Matter, South China Normal University, Guangzhou 510006,
China  Guangdong-Hong Kong Joint Laboratory of Quantum Matter, Southern
Nuclear Science Computing Center, South China Normal University, Guangzhou
510006, China  Hongxi Xing<EMAIL_ADDRESS>Guangdong Provincial Key
Laboratory of Nuclear Science, Institute of Quantum Matter, South China Normal
University, Guangzhou 510006, China  Guangdong-Hong Kong Joint Laboratory of
Quantum Matter, Southern Nuclear Science Computing Center, South China Normal
University, Guangzhou 510006, China

###### Abstract

Studies of exotic hadrons such as the $\chi_{c1}(3872)$ state provide crucial
insights into the fundamental force governing the strong interaction dynamics,
with an emerging new frontier to investigate their production in high energy
collisions where a partonic medium is present. Latest experimental
measurements from the Large Hadron Collider show an intriguing evolution
pattern of the $\chi_{c1}(3872)$-to-$\psi(2S)$ yield ratio from proton-proton
collisions with increasing multiplicities toward proton-lead and lead-lead
collisions. Here we propose a novel mechanism of medium-assisted enhancement
for the $\chi_{c1}(3872)$ production, which competes with the more
conventional absorption-induced suppression and results in a non-monotonic
trend from small to large colliding systems. Realistic simulations from this
model offer the first quantitative description of all available data.
Predictions are made for the centrality dependence of this observable in PbPb
collisions as well as for its system size dependence from OO and ArAr to XeXe
and PbPb collisions. In both cases, a non-monotonic behavior emerges as the
imprint of the competition between enhancement and suppression and can be
readily tested by future data.

Introduction. The overwhelming majority of the energy and mass in the visible
component of our universe comes from the strongly interacting elementary
particles, or hadrons, such as protons, neutrons and pions. According to the
fundamental theory of elementary particles known as the Standard Model, these
hadrons are themselves made from quarks and antiquarks whose interactions are
governed by a basic theory of strong interaction — the Quantum Chromodynamics
(QCD). While the QCD equations are known, their full consequences are
difficult to decipher. One of the outstanding challenges is about the so-
called exotic hadrons, whose quark/antiquark configurations do not follow the
established normal patterns of three quarks forming a baryon (such as the
protons and neutrons) and a quark-anti-quark pair forming a meson (like
pions).

A most notable example of such states is the $\chi_{c1}(3872)$ particle, also
commonly known as $X(3872)$, first discovered by the Belle experimentChoi  _et
al._ (2003) in 2003. Its quark content consists of a pair of charm and anti-
charm quarks as well as another pair of light flavor quark and anti-quark,
i.e. $c\bar{c}q\bar{q}$. Subsequently extensive efforts Acosta  _et al._
(2004); Abazov  _et al._ (2004); Aaij  _et al._ (2020a, b); Ablikim  _et al._
(2020, 2014); Sirunyan  _et al._ (2022, 2022); Lees  _et al._ (2020); Kato
_et al._ (2018); Li and Yuan (2019); Ablikim  _et al._ (2022); Adachi  _et
al._ (2022); LHC (2022); Aubert  _et al._ (2009); Bhardwaj  _et al._ (2011);
Aaij  _et al._ (2014) have been made to measure its quantum numbers as well as
to find many other candidates of exotic hadrons. There are many open questions
in this very active research frontier: What kinds of exotic hadrons could
exist? What are the internal structures and properties of them? What are their
production and decay mechanism?

While properties of $\chi_{c1}(3872)$ have been studied at multiple colliders,
physicists remain puzzled by the nature of this particle despite 20 years past
its initial discovery, with multiple interpretations being proposed for its
internal structures, such as a large-size hadronic molecule versus a compact-
size tetraquark state. See recent reviews in e.g. Brambilla  _et al._ (2011);
Esposito  _et al._ (2015); Briceno  _et al._ (2016); Shepherd  _et al._
(2016); Hosaka  _et al._ (2016); Esposito  _et al._ (2017); Lebed  _et al._
(2017); Olsen  _et al._ (2018); Guo  _et al._ (2018); Wu  _et al._ (2022);
Chen  _et al._ (2023).

A new avenue of investigating exotic hadrons has recently emerged and rapidly
developing, namely to study their formation in high energy hadron and nuclear
collisions where a partonic medium is present. In such collisions, a fireball
with many light flavor quarks/antiquarks (on the order of hundreds to
thousands depending on colliding systems) is created together with a
considerable number of charm quarks/antiquarks. This provides an ideal
environment for creating heavy flavor exotic states and probing their
properties, as demonstrated in the latest theoretical works Cho  _et al._
(2011a, b, 2017); Hong  _et al._ (2018); Fontoura  _et al._ (2019); Zhang  _et
al._ (2021); Esposito  _et al._ (2021); Braaten  _et al._ (2021); Wu  _et al._
(2021); Abreu and Llanes-Estrada (2021); Hu  _et al._ (2021); Chen  _et al._
(2022); Albaladejo  _et al._ (2021); Xu  _et al._ (2021); Abreu  _et al._
(2022a); Jin  _et al._ (2021); Abreu  _et al._ (2022b); Braaten  _et al._
(2022); Yun  _et al._ (2023); Montaña  _et al._ (2022); Hu and Zhang (2023).
Most importantly, experimental measurements of $\chi_{c1}(3872)$ production in
these collisions have started to arrive in the last few years, including LHCb
data from high multiplicity proton-proton (pp) collisions Aaij  _et al._
(2021) and proton-lead (pPb) collisions Durham (2022) as well as CMS data from
lead-lead (PbPb) collisions Sirunyan  _et al._ (2022) at the Large Hadron
Collider (LHC). Already this first batch of empirical information shows an
unusual pattern of the partonic medium’s influence on the $\chi_{c1}(3872)$
yield with respect to the yield of another particle called $\psi(2S)$ which is
a normal hadronic state serving as a benchmark for comparison, by virtue of
its similar heavy flavor content and decay channel ($J/\Psi\pi\pi$) as well as
close mass value to the exotic $\chi_{c1}(3872)$. (These data points are shown
in Fig. 2.) The LHCb pp results suggest the yield ratio of $\chi_{c1}(3872)$
relative to $\psi(2S)$ decreases with increasing event multiplicity, which
would hint at a suppression effect due to the medium. On the other hand, the
LHCb pPb results and the CMS PbPb results, for which the generated medium is
expected to be larger in terms of parton density and system size as compared
with pp collisions, show a strong increase in this observable, which would
indicate an opposite trend to the pp results. So far there has lacked a
consistent explanation that reconciles this intriguing behavior of
$\chi_{c1}(3872)$ production from small to large colliding systems.

In this Letter, we present a phenomenological model for the partonic medium
attenuation effects on the production of $\chi_{c1}(3872)$ in high energy
hadron and nuclear collisions. In particular, a novel mechanism of medium-
assisted enhancement effect will be proposed which competes with the more
conventional absorption-induced suppression effect. Based on this important
feature, it will first be demonstrated qualitatively how the competition leads
to a nontrival pattern in the yield ratio of $\chi_{c1}(3872)$ relative to the
$\psi(2S)$ while the partonic medium evolves from the smaller to the larger
systems. We will then utilize realistic simulations to show how such a model
offers the first quantitative description of all available experimental data.
Predictions will also be made for observables that can be verified in the
future.

Method. In the high $p_{T}$ region where recent CMS and LHCb measurements were
made, the production of $\psi(2S)$ and $\chi_{c1}(3872)$ should dominantly
come from virtual $c\bar{c}$ pairs generated in the initial hard scatterings.
Suppose the number of such pairs that would eventually turn into $\psi(2S)$
and $\chi_{c1}(3872)$, in the absence of any medium effect, would be
$N_{\psi(2S)}$ and $N_{X}$ respectively. However, in nucleus-nucleus (AA) or
high-multiplicity pp and proton-nucleus (pA) collisions, these pairs will need
to first travel through the created partonic medium before producing those
final hadrons. The influence of the medium on the evolution of such $c\bar{c}$
pairs is the focus of our analysis.

The first important effect is the medium absorption. Random collisions with
quarks and gluons from the medium result in the dissociation of the correlated
co-moving $c\bar{c}$ pair, which is akin to the well-known $J/\Psi$
suppression as well as jet quenching that have been observed in AA collisions
Rapp  _et al._ (2010); Rapp (2023). We model this effect as the geometric
absorption along the in-medium path of a $c\bar{c}$ pair:

$\displaystyle\frac{\mathrm{d}N_{i}}{\mathrm{d}x}=-\alpha_{i}n(x)N_{i},$ (1)

where $i\to\psi(2S),\chi_{c1}(3872)$. The $n(x)$ is the local parton density
of the medium along the path of a surviving $c\bar{c}$ pair. The coefficient
$\alpha_{i}$ describes the likelihood of a given state to be dissociated, with
the dimension of a cross-section. While the precise internal structure of
$\chi_{c1}$ is unknown, given the higher mass of $\chi_{c1}$ and its near-
threshold nature, its binding energy is plausibly smaller than that of
$\psi(2S)$. Thus one expects a stronger absorption effect for $\chi_{c1}$,
similarly to the sequential suppression studied for normal charmonium and
bottomonium states Karsch  _et al._ (2006); Chatrchyan  _et al._ (2012);
Aboona  _et al._ (2023); Rapp  _et al._ (2010); Rapp (2023); Du and Rapp
(2015). As is typically done in geometric models for jet energy loss or for
charmonium suppression, one can evaluate the overall suppression by first
integrating the above equation along any given path, then averaging over all
possible in-medium paths, and finally averaging over collision events. This
leads to the following expression for the suppression factor of $\psi(2S)$:

$\displaystyle R^{\psi(2S)}=\langle\langle
e^{-\alpha_{\psi(2S)}\int_{\mathrm{path}}n(x)\mathrm{d}x}\rangle\rangle$ (2)

where the notation $\langle\langle...\rangle\rangle$ means
$\langle\langle...\rangle_{\mathrm{path}}\rangle_{\mathrm{event}}$. We note
that such a suppression effect applies similarly to the $\chi_{c1}(3872)$
production.

For $\chi_{c1}(3872)$, however, there is another medium effect that can
actually help enhance its production. In addition to the $c\bar{c}$, the
formation of $\chi_{c1}(3872)$ requires two light quarks/antiquarks.
Scatterings with the partonic medium, which serves as a reservoir of numerous
light quarks/antiquarks, could lead to “picking up” of light quarks/antiquarks
which then co-move with the $c\bar{c}$ pair. This enhances the probability to
form the $\chi_{c1}(3872)$ state in the end. One could consider this as a two-
step process, in which the $c\bar{c}$ pair picks up the first needed light
parton and subsequently a second needed light parton. Therefore one can model
such a _medium-assisted enhancement_ effect as follows:

$\displaystyle\frac{\mathrm{d}N_{\chi_{c1}}}{\mathrm{d}x}=\beta_{X}n(x)\left[\int_{0}^{x}\beta_{X}n(y)\mathrm{d}y\right]N_{X},$
(3)

where $\beta_{X}$ is a parameter characterizing the probability of picking up
a single light parton, which also has the dimension of a cross-section. An
important feature of this effect is that it scales as square power of the
medium parton density. The enhancement mechanism arises from utilizing light
partons from medium for producing the exotic state along in-medium path of
$c\bar{c}$ pairs, which differs from the regeneration mechanism proposed for
the production of normal quarkonium states Grandchamp and Rapp (2002); Du and
Rapp (2015); Du  _et al._ (2017); Zhou  _et al._ (2014); Rapp  _et al._
(2010).

Combining this enhancement together with the previous suppression effect, one
obtains:

$\displaystyle R^{X}=\langle\langle
e^{\int_{\mathrm{path}}[-\alpha_{X}n(x)+\beta_{X}^{2}n(x)\int_{0}^{x}n(y)\mathrm{d}y]\mathrm{d}x}\rangle\rangle\,.$
(4)

The two parameters $\alpha_{X}$ and $\beta_{X}$ quantify the suppression and
enhancement due to dynamical interactions with medium. Their values should be
sensitive to whether the $\chi_{c1}(3872)$ is a hadronic molecule or a
tetraquark state. Thus extracting them from experimental data can help shed
light on the microscopic structure of $\chi_{c1}(3872)$.

Now we can compare the production of $\chi_{c1}(3872)$ relative to $\psi(2S)$.
This is quantified by the ratio of their baseline pp production cross-section,
modulated by their respective suppression/enhancement effects along the in-
medium paths:

$\displaystyle\frac{N^{X}}{N^{\psi(2S)}}$ $\displaystyle=$
$\displaystyle\frac{\sigma^{X}_{pp}}{\sigma^{\psi(2S)}_{pp}}\times\frac{R^{X}}{R^{\psi(2S)}}$
(5) $\displaystyle\approx$
$\displaystyle\frac{\sigma^{X}_{pp}}{\sigma^{\psi(2S)}_{pp}}\times\mathcal{R}_{med.},$
$\displaystyle\mathcal{R}_{med.}$ $\displaystyle\equiv$
$\displaystyle\langle\langle
e^{\int_{\mathrm{path}}[-(\alpha_{X}-\alpha_{\psi(2S)})n(x)+\beta_{X}^{2}n(x)\int_{0}^{x}n(y)\mathrm{d}y]\mathrm{d}x}\rangle\rangle.$

In the above, the pp baseline $\sigma^{X}_{pp}/\sigma^{\psi(2S)}_{pp}$ could
be inferred from experimental data. We will focus on analyzing the medium
attenuation factor $\mathcal{R}_{med.}$. Clearly, $\mathcal{R}_{med.}>1$
implies an overall medium enhancement while $\mathcal{R}_{med.}<1$ means an
overall medium suppression for the $\chi_{c1}(3872)$ production relative to
the $\psi(2S)$. While $\chi_{c1}(3872)$ and $\psi(2S)$ may have different
formation time, we focus on relatively high $p_{T}$ region where both states
presumably form at late time and assume the same in-medium paths for the
initial $c\bar{c}$ pairs.

Let us first examine the qualitative feature of the partonic medium effect.
Since $\alpha_{X}>\alpha_{\psi(2S)}$ and thus
$\alpha_{X}-\alpha_{\psi(2S)}>0$, the first term in the exponential of
$\mathcal{R}_{med.}$ is a suppression term. Its contribution grows linearly
with the medium parton density and the path length. The second term is an
enhancement term and its contribution grows quadratically with the parton
density and path length. As a result, for relatively low medium density and/or
small medium size, the first term will dominate and therefore the overall
medium effect would be a suppression of $\chi_{c1}(3872)$ relative to
$\psi(2S)$. On the other hand, for relatively high medium density and large
medium size, the second term will dominate and therefore the overall medium
effect would be an enhancement, instead. This nonlinear feature, arising from
the competition between suppression and enhancement, points to a non-monotonic
evolution of medium effect that can help explain and provide a unified
interpretation of the recent measurements by both LHCb and CMS from small to
large colliding systems. For a simple illustration, let us assume that the
average medium effect can be approximated by an average parton density
$\bar{n}$ and average path length $\bar{L}$. Further defining a medium
thickness parameter $\bar{W}\equiv\bar{n}\cdot\bar{L}$, the factor
$\mathcal{R}_{med.}$ can be simplified as

$\displaystyle\mathcal{R}_{med.}=e^{\left[-\left(\alpha_{X}-\alpha_{\psi(2S)}\right)\bar{W}+\frac{1}{2}\beta^{2}_{X}\bar{W}^{2}\right]}.$
(7)

The above result clearly suggests that changing from pp with increasing
multiplicities through pA eventually to AA collisions, the medium thickness
$\bar{W}$ will monotonically increase so that the medium attenuation factor
$\mathcal{R}_{med.}$ should first decrease and then increase, for which a
minimum would occur at
$\bar{W}=\bar{n}\bar{L}=\frac{\alpha_{X}-\alpha_{\psi(2S)}}{\beta^{2}_{X}}$.

Figure 1:  Individual contributions from absorption-induced suppression (blue)
and medium-assisted enhancement (red) as well as the overall
$\mathcal{R}_{med.}$ (black) are plotted as functions of the QGP “brick”
length (in fm unit). The QGP is set at a temperature of $360\rm MeV$
corresponding to an entropy density of $105\rm fm^{-3}$. The dashed line
represents a baseline of $\mathcal{R}_{med.}=1$ in the absence of any medium
effect.

To further demonstrate this feature, let us compute the $\mathcal{R}_{med.}$
for such a simple partonic medium, essentially a QGP “brick” Armesto  _et al._
(2012) with constant and homogeneous temperature which extends along a given
direction with a fixed width. In Fig. 1, the individual contributions from
suppression term (blue) and enhancement term (red) as well as the overall
$\mathcal{R}_{med.}$ (black) are plotted as functions of the QGP “brick”
length, showing the nontrivial decrease-then-increase behavior of
$\mathcal{R}_{med.}$ due to the competition between suppression and
enhancement. This behavior already qualitatively agrees with the trends seen
in experimental data. Of course, the partonic medium created in those
collisions is much more complicated than the simple approximation here, due to
nontrivial initial conditions, event-by-event fluctuations, dynamical
expansions, etc. To fully verify the feasibility of this idea, one needs to
perform quantitative and realistic simulations, which we report next.

One important issue not discussed so far is the kinematic factors such as
transverse momentum. The geometric path approach here is suitable for
describing a rapidly traversing probe and best applicable for extracting the
integrated average medium effects on the relatively high $p_{T}$ region. A
more differential treatment that implements initial $c\bar{c}$ $p_{T}$
spectrum as well as $p_{T}$-dependence for both suppression and enhancement
effects, will be a future goal.

Results. To quantitatively evaluate the medium effect on the $\chi_{c1}(3872)$
production relative to the $\psi(2S)$, we utilize event-by-event simulations
based on the iEBE-VISHNU hydrodynamic model Shen  _et al._ (2016a) which has
been well tested by experimental data from small to large colliding systems
Plumberg (2020); Taghavi (2021); Zhao  _et al._ (2020); Aad  _et al._ (2023);
Shen  _et al._ (2016b); Lacey and Magdy (2022); Oliinychenko and Shen (2021);
Wang  _et al._ (2023); Garcia Montero  _et al._ (2021); Dai  _et al._ (2020);
Mitrankov  _et al._ (2022); Abbasi  _et al._ (2018); Zhao  _et al._ (2018a);
Zhu (2016); Lacey  _et al._ (2022); Xu  _et al._ (2023); Heinz and Liu (2016);
Ru  _et al._ (2018); Chang  _et al._ (2016); Li  _et al._ (2018); Mehrabpour
and Taghavi (2019); Li  _et al._ (2022); Zhu  _et al._ (2015); Mordasini  _et
al._ (2020); Mitrankov  _et al._ (2021); Zhu  _et al._ (2017); Aronson  _et
al._ (2018); Song  _et al._ (2017); Fu (2015); Acharya  _et al._ (2018a);
Devetak (2021); Shen  _et al._ (2016c); Bernhard  _et al._ (2016); Acharya
_et al._ (2020, 2018b, 2019); Pang  _et al._ (2018); Song  _et al._ (2017);
Zhao  _et al._ (2017); Adam  _et al._ (2016); Qian  _et al._ (2016); Zhao  _et
al._ (2018b); Acharya  _et al._ (2018a); Zhao  _et al._ (2018a); Xu  _et al._
(2016); Du  _et al._ (2020); Shen  _et al._ (2015). The iEBE-VISHNU provides
event-wise time-dependent evolution information of the bulk medium as well as
initial conditions for pp, pA, and AA collisions. We generate $c\bar{c}$ pairs
at different spots on the event plane according to the initial binary
collision density profiles. The pairs then move along straight paths whose
directions are randomly chosen, and along each path the bulk medium is
evolving in time. At each spacetime point the local entropy density can be
read from the iEBE-VISHNU hydrodynamic evolution, which is directly related to
active degrees of freedom in the thermal medium and thus can be used to
represent the local parton density up to a proportionality constant. In this
work, we simulated a total of 500,000 events for pp collisions at
$\sqrt{S_{NN}}=8$ TeV, 200,000 events for pPb collisions at
$\sqrt{S_{NN}}=8.16$ TeV and 100,000 events for PbPb collisions at
$\sqrt{S_{NN}}=5.02$ TeV. Model parameters of the hydro code were set in the
same way as previous studies based on the same package in the literature, see
e.g. Zhu  _et al._ (2017); Li  _et al._ (2018); Zhao  _et al._ (2020) . The
bulk properties such as multiplicity and elliptic flow from our simulations
compare well with experimental measurements for all colliding systems. For
each event, we further simulate about 100 to 10000 in-medium paths depending
on the medium size. Due to the substantial amount of needed computing time, we
chose to simplify the calculations by first performing average over the path
integrations in each event and then computing the exponential for the ratio
between $\chi_{c1}(3872)$ and $\psi(2S)$. That is:

$\displaystyle\mathcal{R}_{med.}\approx\langle e^{-\alpha^{\prime}\cdot
P_{1}+\beta^{{}^{\prime}2}\cdot P_{2}}\rangle_{\mathrm{event}},$ (8)

where

$\displaystyle P_{1}$ $\displaystyle=$
$\displaystyle\left\langle\int_{\mathrm{path}}s(x)\mathrm{d}x\right\rangle_{\mathrm{path}},$
(9) $\displaystyle P_{2}$ $\displaystyle=$
$\displaystyle\left\langle\int_{\mathrm{path}}s(x)\left(\int_{0}^{x}s(y)\mathrm{d}y\right)\mathrm{d}x\right\rangle_{\mathrm{path}}.$
(10)

where we introduce $s(x)$ as the entropy density, and
$\alpha^{\prime}=(\alpha_{X}-\alpha_{\psi(2S)})\left(\frac{n}{s}\right)$,
$\beta^{\prime}=\beta_{X}\left(\frac{n}{s}\right)$ whose definitions absorb
the proportionality constant between parton density and entropy density.

To determine the two key parameters $\alpha^{\prime}$ and $\beta^{\prime}$, we
take the LHCb pp (at $\sqrt{S_{NN}}=8$ TeV) and preliminary pPb (at
$\sqrt{S_{NN}}=8.16$TeV) as well as the CMS PbPb (at $\sqrt{S_{NN}}=5.02$ TeV)
data for a global fitting analysis, with results shown in Fig. 2. The best
fit, with $\chi^{2}/d.o.f=1.78$, gives the following numbers with $1\sigma$
level uncertainty: $\alpha^{\prime}=(5.7\pm 2.2)\times 10^{-3}fm^{2}$,
$\beta^{\prime}=(5.7\pm 0.7)\times 10^{-3}fm^{2}$ and
$\sigma^{X}_{pp}/\sigma^{\psi(2S)}_{pp}=0.162\pm 0.037$. Note the analysis
does not assume the sign of $\alpha^{\prime}$ which is found to be positive
from the best fit. As one can see, our model with just two parameters
characterizing a competition between suppression and enhancement can well
describe the quantitative trends of all global data from small to large
systems. The newly proposed medium-assisted enhancement is particularly
important for understanding the rapid increase of $\chi_{c1}(3872)$ yield
relative to $\psi(2S)$ in the pPb and PbPb collisions. The phenomenological
extraction of the two parameters from experimental data can serve as important
constraints on the nature of $\chi_{c1}(3872)$. Microscopic calculations of
these parameters based on different assumptions (e.g. hadronic molecule versus
tetraquark state) could then be compared with the empirical values obtained
here to help decipher its internal structures, which will be an interesting
future project.

Figure 2:  A comparison of the $\chi_{c1}(3872)$ yield relative to $\psi(2S)$
between model simulation results (blue curve) and experimental data from LHCb
pp collisions at $\sqrt{S_{NN}}=8$ TeV (red circle), LHCb preliminary pPb
collisions at $\sqrt{S_{NN}}=8.16$ TeV (orange triangle) and CMS PbPb
collisions at $\sqrt{S_{NN}}=5.02$ TeV (green box) Sirunyan  _et al._ (2022);
Aaij  _et al._ (2021); Durham (2022). The model parameters are determined from
the global fitting analysis with the blue band showing the $1\sigma$ level
uncertainty. (See text for details.)

A natural next step would be testing our model predictions for which
experimental data are not yet available and can serve as a future validation.
To do that, we further investigate the centrality dependence of the medium
effect in the PbPb collisions. In Fig 3, we show the yield ratio of
$\chi_{c1}(3872)$ to $\psi(2S)$ in three centrality bins: $60-90\%$,
$30-60\%$, and $0-30\%$. Interestingly, the results again show a non-monotonic
behavior. In peripheral collisions the ratio is about 0.1, while in middle-
centrality class it will drop to around 0.02 which is comparable to that in
the pPb collisions. Finally in the central collisions, it rockets up to be as
high as 3. This finding also suggests that the most central collisions
actually contribute most of the $\chi_{c1}(3872)$ particles observed in the
minimal bias measurements from CMS. We emphasize that the model parameters
were already fixed in the aforementioned fitting analysis, so the highly
nontrivial centrality trend predicted by the model here will be an important
verification by future measurements.

Figure 3: The predicted centrality dependence of the $\chi_{c1}(3872)$ yield
relative to $\psi(2S)$ in PbPb collisions at $\sqrt{S_{NN}}=5.02$TeV
collisions. The blue uncertainty band is from the same source as in Fig. 2.
Figure 4:  The predicted trend of the $\chi_{c1}(3872)$ yield relative to
$\psi(2S)$ in AA collisions from small to large systems, including: OO
collisions at $\sqrt{S_{NN}}=6.5$TeV, ArAr collisions at
$\sqrt{S_{NN}}=5.85$TeV, XeXe collisions at $\sqrt{S_{NN}}=5.44$TeV, as well
as PbPb collisions at $\sqrt{S_{NN}}=5.02$TeV. The blue uncertainty band is
from the same source as in Fig. 2.

Finally, the system size scan for AA collisions could offer yet another
independent validation of our model predictions. For that purpose, we have
computed the $\chi_{c1}(3872)$ to $\psi(2S)$ yield ratio for the following
systems: OO collisions at $\sqrt{S_{NN}}=6.5$ TeV, ArAr collisions at
$\sqrt{S_{NN}}=5.85$ TeV, and XeXe collisions at $\sqrt{S_{NN}}=5.44$ TeV.
These results are shown in Fig. 4 in comparison with PbPb collisions at
$\sqrt{S_{NN}}=5.02$ TeV. Again, one observes a nontrivial trend that first
decreases and then increases when changing from smaller to larger colliding
systems. Such prediction of the system size dependence in AA collisions, which
is the consequence of competing suppression and enhancement effects, can be
readily tested with future measurements.

Conclusion. To conclude, we present a phenomenological model for the partonic
medium attenuation effects on the production of $\chi_{c1}(3872)$ and
$\psi(2S)$ particles in high energy hadron and nuclear collisions. In
particular, a novel mechanism of medium-assisted enhancement effect is
proposed for the $\chi_{c1}(3872)$ production, which leads to a competition
with the more conventional absorption-induced suppression effect and becomes
more dominant for higher parton densities and larger medium size. As a
consequence of this important feature, the yield ratio of $\chi_{c1}(3872)$
relative to the $\psi(2S)$ develops a nontrivial pattern, first decreasing
then increasing, when the partonic medium evolves from small to large
colliding systems. Utilizing realistic simulations, we show that this model
offers the first quantitative description of all available experimental
measurements, including the LHCb pp (at $\sqrt{S_{NN}}=8$ TeV) and preliminary
pPb (at $\sqrt{S_{NN}}=8.16$TeV) as well as the CMS PbPb (at
$\sqrt{S_{NN}}=5.02$ TeV) data. We further make predictions for the centrality
dependence of the $\chi_{c1}(3872)$-to-$\psi(2S)$ yield ratio in PbPb
collisions as well as for its system size dependence from OO and ArAr to XeXe
and PbPb collisions. In both cases, a non-monotonic pattern emerges as the
imprint of the competition between enhancement and suppression. Given the
expected abundance of experimental data from planned runs as well as
anticipated upgrades at the LHC, it would be exciting to test these
predictions with future high precision measurements.

###### Acknowledgements.

Acknowledgments. This research was supported in part by the National Natural
Science Foundation of China (NSFC) under Grants No. 12035007, No. 12022512 and
No. 11905066, by Guangdong Major Project of Basic and Applied Basic Research
No. 2020B0301030008, by the National Science Foundation in US under Grant No.
PHY-2209183 (J.L.), and by the DOE through ExoHad Topical Collaboration.

## References

  * Choi  _et al._ (2003) S. K. Choi _et al._ (Belle), Phys. Rev. Lett. 91, 262001 (2003), arXiv:hep-ex/0309032 .
  * Acosta  _et al._ (2004) D. Acosta _et al._ (CDF), Phys. Rev. Lett. 93, 072001 (2004), arXiv:hep-ex/0312021 .
  * Abazov  _et al._ (2004) V. M. Abazov _et al._ (D0), Phys. Rev. Lett. 93, 162002 (2004), arXiv:hep-ex/0405004 .
  * Aaij  _et al._ (2020a) R. Aaij _et al._ (LHCb), Phys. Rev. D 102, 092005 (2020a), arXiv:2005.13419 [hep-ex] .
  * Aaij  _et al._ (2020b) R. Aaij _et al._ (LHCb), JHEP 08, 123 (2020b), arXiv:2005.13422 [hep-ex] .
  * Ablikim  _et al._ (2020) M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 124, 242001 (2020), arXiv:2001.01156 [hep-ex] .
  * Ablikim  _et al._ (2014) M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 112, 092001 (2014), arXiv:1310.4101 [hep-ex] .
  * Sirunyan  _et al._ (2022) A. M. Sirunyan _et al._ (CMS), Phys. Rev. Lett. 128, 032001 (2022), arXiv:2102.13048 [hep-ex] .
  * Lees  _et al._ (2020) J. P. Lees _et al._ (BaBar), Phys. Rev. Lett. 124, 152001 (2020), arXiv:1911.11740 [hep-ex] .
  * Kato  _et al._ (2018) Y. Kato _et al._ (Belle), Phys. Rev. D 97, 012005 (2018), arXiv:1709.06108 [hep-ex] .
  * Li and Yuan (2019) C. Li and C.-Z. Yuan, Phys. Rev. D 100, 094003 (2019), arXiv:1907.09149 [hep-ex] .
  * Ablikim  _et al._ (2022) M. Ablikim _et al._ (BESIII), Phys. Rev. D 105, 072009 (2022), arXiv:2202.00759 [hep-ex] .
  * Adachi  _et al._ (2022) I. Adachi _et al._ (Belle), (2022), arXiv:2206.08592 [hep-ex] .
  * LHC (2022) (2022), arXiv:2204.12597 [hep-ex] .
  * Aubert  _et al._ (2009) B. Aubert _et al._ (BaBar), Phys. Rev. Lett. 102, 132001 (2009), arXiv:0809.0042 [hep-ex] .
  * Bhardwaj  _et al._ (2011) V. Bhardwaj _et al._ (Belle), Phys. Rev. Lett. 107, 091803 (2011), arXiv:1105.0177 [hep-ex] .
  * Aaij  _et al._ (2014) R. Aaij _et al._ (LHCb), Nucl. Phys. B 886, 665 (2014), arXiv:1404.0275 [hep-ex] .
  * Brambilla  _et al._ (2011) N. Brambilla _et al._ , Eur. Phys. J. C 71, 1534 (2011), arXiv:1010.5827 [hep-ph] .
  * Esposito  _et al._ (2015) A. Esposito, A. L. Guerrieri, F. Piccinini, A. Pilloni, and A. D. Polosa, Int. J. Mod. Phys. A 30, 1530002 (2015), arXiv:1411.5997 [hep-ph] .
  * Briceno  _et al._ (2016) R. A. Briceno _et al._ , Chin. Phys. C 40, 042001 (2016), arXiv:1511.06779 [hep-ph] .
  * Shepherd  _et al._ (2016) M. R. Shepherd, J. J. Dudek, and R. E. Mitchell, Nature 534, 487 (2016), arXiv:1802.08131 [hep-ph] .
  * Hosaka  _et al._ (2016) A. Hosaka, T. Iijima, K. Miyabayashi, Y. Sakai, and S. Yasui, PTEP 2016, 062C01 (2016), arXiv:1603.09229 [hep-ph] .
  * Esposito  _et al._ (2017) A. Esposito, A. Pilloni, and A. D. Polosa, Phys. Rept. 668, 1 (2017), arXiv:1611.07920 [hep-ph] .
  * Lebed  _et al._ (2017) R. F. Lebed, R. E. Mitchell, and E. S. Swanson, Prog. Part. Nucl. Phys. 93, 143 (2017), arXiv:1610.04528 [hep-ph] .
  * Olsen  _et al._ (2018) S. L. Olsen, T. Skwarnicki, and D. Zieminska, Rev. Mod. Phys. 90, 015003 (2018), arXiv:1708.04012 [hep-ph] .
  * Guo  _et al._ (2018) F.-K. Guo, C. Hanhart, U.-G. Meißner, Q. Wang, Q. Zhao, and B.-S. Zou, Rev. Mod. Phys. 90, 015004 (2018), [Erratum: Rev.Mod.Phys. 94, 029901 (2022)], arXiv:1705.00141 [hep-ph] .
  * Wu  _et al._ (2022) T.-W. Wu, Y.-W. Pan, M.-Z. Liu, and L.-S. Geng, Sci. Bull. 67, 1735 (2022), arXiv:2208.00882 [hep-ph] .
  * Chen  _et al._ (2023) H.-X. Chen, W. Chen, X. Liu, Y.-R. Liu, and S.-L. Zhu, Rept. Prog. Phys. 86, 026201 (2023), arXiv:2204.02649 [hep-ph] .
  * Cho  _et al._ (2011a) S. Cho _et al._ (ExHIC), Phys. Rev. Lett. 106, 212001 (2011a), arXiv:1011.0852 [nucl-th] .
  * Cho  _et al._ (2011b) S. Cho _et al._ (ExHIC), Phys. Rev. C 84, 064910 (2011b), arXiv:1107.1302 [nucl-th] .
  * Cho  _et al._ (2017) S. Cho _et al._ (ExHIC), Prog. Part. Nucl. Phys. 95, 279 (2017), arXiv:1702.00486 [nucl-th] .
  * Hong  _et al._ (2018) J. Hong, S. Cho, T. Song, and S. H. Lee, Phys. Rev. C 98, 014913 (2018), arXiv:1804.05336 [nucl-th] .
  * Fontoura  _et al._ (2019) C. E. Fontoura, G. Krein, A. Valcarce, and J. Vijande, Phys. Rev. D 99, 094037 (2019), arXiv:1905.03877 [hep-ph] .
  * Zhang  _et al._ (2021) H. Zhang, J. Liao, E. Wang, Q. Wang, and H. Xing, Phys. Rev. Lett. 126, 012301 (2021), arXiv:2004.00024 [hep-ph] .
  * Esposito  _et al._ (2021) A. Esposito, E. G. Ferreiro, A. Pilloni, A. D. Polosa, and C. A. Salgado, Eur. Phys. J. C 81, 669 (2021), arXiv:2006.15044 [hep-ph] .
  * Braaten  _et al._ (2021) E. Braaten, L.-P. He, K. Ingles, and J. Jiang, Phys. Rev. D 103, L071901 (2021), arXiv:2012.13499 [hep-ph] .
  * Wu  _et al._ (2021) B. Wu, X. Du, M. Sibila, and R. Rapp, Eur. Phys. J. A 57, 122 (2021), [Erratum: Eur.Phys.J.A 57, 314 (2021)], arXiv:2006.09945 [nucl-th] .
  * Abreu and Llanes-Estrada (2021) L. M. Abreu and F. J. Llanes-Estrada, Eur. Phys. J. C 81, 430 (2021), arXiv:2008.12031 [hep-ph] .
  * Hu  _et al._ (2021) Y. Hu, J. Liao, E. Wang, Q. Wang, H. Xing, and H. Zhang, Phys. Rev. D 104, L111502 (2021), arXiv:2109.07733 [hep-ph] .
  * Chen  _et al._ (2022) B. Chen, L. Jiang, X.-H. Liu, Y. Liu, and J. Zhao, Phys. Rev. C 105, 054901 (2022), arXiv:2107.00969 [hep-ph] .
  * Albaladejo  _et al._ (2021) M. Albaladejo, J. M. Nieves, and L. Tolos, Phys. Rev. C 104, 035203 (2021), arXiv:2102.08589 [hep-ph] .
  * Xu  _et al._ (2021) H.-g. Xu, Z.-L. She, D.-M. Zhou, L. Zheng, X.-L. Kang, G. Chen, and B.-H. Sa, Eur. Phys. J. C 81, 784 (2021), arXiv:2105.06261 [hep-ph] .
  * Abreu  _et al._ (2022a) L. M. Abreu, F. S. Navarra, M. Nielsen, and H. P. L. Vieira, Eur. Phys. J. C 82, 296 (2022a), arXiv:2110.11145 [hep-ph] .
  * Jin  _et al._ (2021) Y. Jin, S.-Y. Li, Y.-R. Liu, Q. Qin, Z.-G. Si, and F.-S. Yu, Phys. Rev. D 104, 114009 (2021), arXiv:2109.05678 [hep-ph] .
  * Abreu  _et al._ (2022b) L. M. Abreu, H. P. L. Vieira, and F. S. Navarra, Phys. Rev. D 105, 116029 (2022b), arXiv:2202.10882 [hep-ph] .
  * Braaten  _et al._ (2022) E. Braaten, R. Bruschini, L.-P. He, K. Ingles, and J. Jiang, (2022), arXiv:2209.04972 [hep-ph] .
  * Yun  _et al._ (2023) H. Yun, D. Park, S. Noh, A. Park, W. Park, S. Cho, J. Hong, Y. Kim, S. Lim, and S. H. Lee, Phys. Rev. C 107, 014906 (2023), arXiv:2208.06960 [hep-ph] .
  * Montaña  _et al._ (2022) G. Montaña, A. Ramos, L. Tolos, and J. M. Torres-Rincon, (2022), arXiv:2211.01896 [hep-ph] .
  * Hu and Zhang (2023) Y. Hu and H. Zhang, (2023), arXiv:2301.06441 [hep-ph] .
  * Aaij  _et al._ (2021) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 126, 092001 (2021), arXiv:2009.06619 [hep-ex] .
  * Durham (2022) J. M. Durham (LHCb), LHCb-CONF-2022-001, CERN-LHCb-CONF-2022-001 (2022).
  * Rapp  _et al._ (2010) R. Rapp, D. Blaschke, and P. Crochet, Prog. Part. Nucl. Phys. 65, 209 (2010), arXiv:0807.2470 [hep-ph] .
  * Rapp (2023) R. Rapp, Nucl. Sci. Tech. 34, 63 (2023).
  * Karsch  _et al._ (2006) F. Karsch, D. Kharzeev, and H. Satz, Phys. Lett. B 637, 75 (2006), arXiv:hep-ph/0512239 .
  * Chatrchyan  _et al._ (2012) S. Chatrchyan _et al._ (CMS), Phys. Rev. Lett. 109, 222301 (2012), [Erratum: Phys.Rev.Lett. 120, 199903 (2018)], arXiv:1208.2826 [nucl-ex] .
  * Aboona  _et al._ (2023) B. Aboona _et al._ (STAR), Phys. Rev. Lett. 130, 112301 (2023), arXiv:2207.06568 [nucl-ex] .
  * Du and Rapp (2015) X. Du and R. Rapp, Nucl. Phys. A 943, 147 (2015), arXiv:1504.00670 [hep-ph] .
  * Grandchamp and Rapp (2002) L. Grandchamp and R. Rapp, Nucl. Phys. A 709, 415 (2002), arXiv:hep-ph/0205305 .
  * Du  _et al._ (2017) X. Du, M. He, and R. Rapp, Phys. Rev. C 96, 054901 (2017), arXiv:1706.08670 [hep-ph] .
  * Zhou  _et al._ (2014) K. Zhou, N. Xu, Z. Xu, and P. Zhuang, Phys. Rev. C 89, 054911 (2014), arXiv:1401.5845 [nucl-th] .
  * Armesto  _et al._ (2012) N. Armesto _et al._ , Phys. Rev. C 86, 064904 (2012), arXiv:1106.1106 [hep-ph] .
  * Shen  _et al._ (2016a) C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, and U. Heinz, Comput. Phys. Commun. 199, 61 (2016a), arXiv:1409.8164 [nucl-th] .
  * Plumberg (2020) C. Plumberg, (2020), arXiv:2010.11957 [nucl-th] .
  * Taghavi (2021) S. F. Taghavi, Phys. Rev. C 104, 054906 (2021), arXiv:1907.12140 [nucl-th] .
  * Zhao  _et al._ (2020) W. Zhao, Y. Zhou, K. Murase, and H. Song, Eur. Phys. J. C 80, 846 (2020), arXiv:2001.06742 [nucl-th] .
  * Aad  _et al._ (2023) G. Aad _et al._ (ATLAS), JHEP 07, 074 (2023), arXiv:2211.15257 [hep-ex] .
  * Shen  _et al._ (2016b) C. Shen, J. F. Paquet, G. S. Denicol, S. Jeon, and C. Gale, Phys. Rev. Lett. 116, 072301 (2016b), arXiv:1504.07989 [nucl-th] .
  * Lacey and Magdy (2022) R. A. Lacey and N. Magdy, (2022), arXiv:2206.05773 [nucl-ex] .
  * Oliinychenko and Shen (2021) D. Oliinychenko and C. Shen, (2021), arXiv:2105.07539 [hep-ph] .
  * Wang  _et al._ (2023) S. Wang, W. Dai, B.-W. Zhang, and E. Wang, Chin. Phys. C 47, 054102 (2023), arXiv:2005.07018 [hep-ph] .
  * Garcia Montero  _et al._ (2021) O. Garcia Montero, N. Löher, A. Mazeliauskas, K. Reygers, and J. Berges, PoS HardProbes2020, 033 (2021).
  * Dai  _et al._ (2020) W. Dai, S. Wang, S.-L. Zhang, B.-W. Zhang, and E. Wang, Chin. Phys. C 44, 104105 (2020), arXiv:1806.06332 [nucl-th] .
  * Mitrankov  _et al._ (2022) I. M. Mitrankov, Y. A. Berdnikov, A. Y. Berdnikov, D. O. Kotov, M. M. Mitrankova, and N. Novitzky (PHENIX), Phys. Scripta 97, 054012 (2022).
  * Abbasi  _et al._ (2018) N. Abbasi, D. Allahbakhshi, A. Davody, and S. F. Taghavi, Phys. Rev. C 98, 024906 (2018), arXiv:1704.06295 [nucl-th] .
  * Zhao  _et al._ (2018a) W. Zhao, L. Zhu, H. Zheng, C. M. Ko, and H. Song, Phys. Rev. C 98, 054905 (2018a), arXiv:1807.02813 [nucl-th] .
  * Zhu (2016) X. Zhu, Adv. High Energy Phys. 2016, 4236492 (2016), arXiv:1607.04003 [nucl-th] .
  * Lacey  _et al._ (2022) R. A. Lacey, N. Magdy, P. Parfenov, and A. Taranenko, (2022), arXiv:2203.10029 [nucl-ex] .
  * Xu  _et al._ (2023) H.-j. Xu, W. Zhao, H. Li, Y. Zhou, L.-W. Chen, and F. Wang, Phys. Rev. C 108, L011902 (2023), arXiv:2111.14812 [nucl-th] .
  * Heinz and Liu (2016) U. W. Heinz and J. Liu, Nucl. Phys. A 956, 549 (2016), arXiv:1512.08276 [nucl-th] .
  * Ru  _et al._ (2018) P. Ru, G. Bary, and W.-N. Zhang, Phys. Lett. B 777, 79 (2018), arXiv:1707.05461 [nucl-th] .
  * Chang  _et al._ (2016) N.-b. Chang _et al._ , Sci. China Phys. Mech. Astron. 59, 621001 (2016), arXiv:1510.05754 [nucl-th] .
  * Li  _et al._ (2018) J. Li, H.-j. Xu, and H. Song, Phys. Rev. C 97, 014902 (2018), arXiv:1707.09742 [nucl-th] .
  * Mehrabpour and Taghavi (2019) H. Mehrabpour and S. F. Taghavi, Eur. Phys. J. C 79, 88 (2019), arXiv:1805.04695 [nucl-th] .
  * Li  _et al._ (2022) Y. Li, S. Wang, and B.-W. Zhang, (2022), arXiv:2209.00548 [hep-ph] .
  * Zhu  _et al._ (2015) X. Zhu, F. Meng, H. Song, and Y.-X. Liu, Phys. Rev. C 91, 034904 (2015), arXiv:1501.03286 [nucl-th] .
  * Mordasini  _et al._ (2020) C. Mordasini, A. Bilandzic, D. Karakoç, and S. F. Taghavi, Phys. Rev. C 102, 024907 (2020), arXiv:1901.06968 [nucl-ex] .
  * Mitrankov  _et al._ (2021) I. M. Mitrankov, E. V. Bannikov, A. Y. Berdnikov, Y. A. Berdnikov, and D. O. Kotov, J. Phys. Conf. Ser. 2103, 012133 (2021).
  * Zhu  _et al._ (2017) X. Zhu, Y. Zhou, H. Xu, and H. Song, Phys. Rev. C 95, 044902 (2017), arXiv:1608.05305 [nucl-th] .
  * Aronson  _et al._ (2018) S. Aronson, E. Borras, B. Odegard, R. Sharma, and I. Vitev, Phys. Lett. B 778, 384 (2018), arXiv:1709.02372 [hep-ph] .
  * Song  _et al._ (2017) H. Song, Y. Zhou, and K. Gajdosova, Nucl. Sci. Tech. 28, 99 (2017), arXiv:1703.00670 [nucl-th] .
  * Fu (2015) J. Fu, Phys. Rev. C 92, 024904 (2015).
  * Acharya  _et al._ (2018a) S. Acharya _et al._ (ALICE), Eur. Phys. J. C 78, 624 (2018a), arXiv:1801.07051 [nucl-ex] .
  * Devetak (2021) D. Devetak, Eur. Phys. J. D 75, 66 (2021).
  * Shen  _et al._ (2016c) C. Shen, C. Park, J.-F. Paquet, G. S. Denicol, S. Jeon, and C. Gale, Nucl. Phys. A 956, 741 (2016c), arXiv:1601.03070 [hep-ph] .
  * Bernhard  _et al._ (2016) J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Phys. Rev. C 94, 024907 (2016), arXiv:1605.03954 [nucl-th] .
  * Acharya  _et al._ (2020) S. Acharya _et al._ (ALICE), Phys. Rev. C 101, 044907 (2020), arXiv:1910.07678 [nucl-ex] .
  * Acharya  _et al._ (2018b) S. Acharya _et al._ (ALICE), JHEP 09, 006 (2018b), arXiv:1805.04390 [nucl-ex] .
  * Acharya  _et al._ (2019) S. Acharya _et al._ (ALICE), Phys. Lett. B 790, 35 (2019), arXiv:1805.04432 [nucl-ex] .
  * Pang  _et al._ (2018) L.-G. Pang, K. Zhou, N. Su, H. Petersen, H. Stöcker, and X.-N. Wang, Nature Commun. 9, 210 (2018), arXiv:1612.04262 [hep-ph] .
  * Zhao  _et al._ (2017) W. Zhao, H.-j. Xu, and H. Song, Eur. Phys. J. C 77, 645 (2017), arXiv:1703.10792 [nucl-th] .
  * Adam  _et al._ (2016) J. Adam _et al._ (ALICE), JHEP 09, 164 (2016), arXiv:1606.06057 [nucl-ex] .
  * Qian  _et al._ (2016) J. Qian, U. W. Heinz, and J. Liu, Phys. Rev. C 93, 064901 (2016), arXiv:1602.02813 [nucl-th] .
  * Zhao  _et al._ (2018b) W. Zhao, Y. Zhou, H. Xu, W. Deng, and H. Song, Phys. Lett. B 780, 495 (2018b), arXiv:1801.00271 [nucl-th] .
  * Xu  _et al._ (2016) H.-j. Xu, Z. Li, and H. Song, Phys. Rev. C 93, 064905 (2016), arXiv:1602.02029 [nucl-th] .
  * Du  _et al._ (2020) Y.-L. Du, K. Zhou, J. Steinheimer, L.-G. Pang, A. Motornenko, H.-S. Zong, X.-N. Wang, and H. Stöcker, Eur. Phys. J. C 80, 516 (2020), arXiv:1910.11530 [hep-ph] .
  * Shen  _et al._ (2015) C. Shen, Z. Qiu, and U. Heinz, Phys. Rev. C 92, 014901 (2015), arXiv:1502.04636 [nucl-th] .