# Sparsity exploitation via discovering graphical models in multi-variate
time-series forecasting

Ngoc-Dung Do 1 3, Truong Son Hy 2 3 4, Duy Khuong Nguyen 1 1 FPT Software AI
Center, Hanoi 10000, Vietnam2 University of California San Diego, La Jolla, CA
92093, USA3 Co-first authors4 Correspondence to<EMAIL_ADDRESS>
(April 2023)

###### Abstract

Graph neural networks (GNNs) have been widely applied in multi-variate time-
series forecasting (MTSF) tasks because of their capability in capturing the
correlations among different time-series. These graph-based learning
approaches improve the forecasting performance by discovering and
understanding the underlying graph structures, which represent the data
correlation. When the explicit prior graph structures are not available, most
existing works cannot guarantee the sparsity of the generated graphs that make
the overall model computational expensive and less interpretable. In this
work, we propose a decoupled training method, which includes a graph
generating module and a GNNs forecasting module. First, we use Graphical Lasso
(or GraphLASSO) to directly exploit the sparsity pattern from data to build
graph structures in both static and time-varying cases. Second, we fit these
graph structures and the input data into a Graph Convolutional Recurrent
Network (GCRN) to train a forecasting model. The experimental results on three
real-world datasets show that our novel approach has competitive performance
against existing state-of-the-art forecasting algorithms while providing
sparse, meaningful and explainable graph structures and reducing training time
by approximately $40\%$. Our PyTorch implementation is publicly available at
https://github.com/HySonLab/GraphLASSO.

###### Index Terms:

Graph Neural Networks, sparse graph learning, graphical lasso, multi-variate
time-series forecasting, spatio-temporal modeling.

##  I Introduction

Many real world problems can be formed in Multi-variate Time Series
Forecasting (MTS) including traffic load forecasting [1, 2], epidemic modeling
[3], retail sales [4], finance [5], etc. that involve recorded temporal data
from multiple sensors being processed and aligned into the same time axis.
Thus, MTS data can be seen as a type of spatial-temporal data, which contains
multiple inter-related time series. A forecasting model needs to capture both
spatial and temporal information in order to perform well in predicting for
future values.

Previous work on MTS forecasting can be divided into statistics-based and deep
learning-based approaches. The former often consider the linear correlations
among variables (i.e. time series) [6, 7]. Regarding the latter, early works
have utilized Convolution Neural Networks (CNN) to capture the correlations
among variables, yet ignored their non-Euclidean pairwise dependencies.
Transformer-based approaches consist of spatial attention mechanism to
discover the reciprocal salience of different spatial locations at each layer
[8, 9]. Recently, Spatial-Temporal Graph Neural Networks (STGNNs) [10, 11, 2],
have attracted an increasing attention because these methods not only archive
outstanding performance on forecasting tasks but are also partially
interpretable based on the graph structures, [12, 13, 14]. In general, STGNNs
combine Graph Convolutional Networks (GCN) [15, 16] to explore the
relationship between time series and sequential models to capture the temporal
dependencies. In this case, each time series is one node of the graph and the
relationship between two time series is represented by an edge. In practice,
these graphs need to be learnt from observations rather than to be pre-
defined, and to be sparse in number of connections.

In the deep learning setting, end-to-end training methods often incorporate a
graph learning module with a forecasting module to maximize performance on the
downstream tasks, for example, GTS [17] and MTGNN [14]. This approach is
particularly sensitive with hyper-parameters and cannot guarantee the sparse
property. Beside, end-to-end training may lead to difficulty in convergence
and usually generates dense graphs. In short, most of the aforementioned
graph-learning methods cannot ensure the sparsity, which is an importance
constraint in real world problems.

In this paper, we propose a simple but reliable approach that can efficiently
capture the graph structures as well as guarantee the sparsity of their
adjacency matrices by using GraphLASSO [18, 19]. GraphLASSO is a sparse
penalized maximum likelihood estimator for the precision matrix of a
multivariate Gaussian distribution, which encodes conditional correlations
between pairs of variables given the remaining variables. It not only provides
an effective way to control the sparsity of the relationships but also reduces
the training time by rapidly generating the graphs representing correlations
without end-to-end training. In summary, our work has three main contributions
as following:

  * •

A method of decoupled training that has a separated process for graph
generation via GraphLASSO,

  * •

Constructing graph structures and exploiting the spatial features of MTS for
both static and times-varying cases,

  * •

Extensive experiments in three real-world datasets have shown that our
approach gives competitive results while requiring significantly less training
time in comparison with other state-of-the-art methods.

##  II Related works

In this paper, we limit the scope to using correlation information between
series to improve forecasting performance. Classical statistical methods, such
as AutoRegressive Integrated Moving Average (ARIMA), are commonly used for
univariate time series prediction, and its variant VAR is often used for
multi-variate time series. In VAR, the relationship between series is
represented by a matrix and therefore is linear. Deep learning-based methods
such as LSTM [20], LSTNet [21] and TPA-LSTM [22] can capture non-linear
dependencies. These methods have a module that learns spatial information and
encodes all of it on a single hidden state vector. This leads to the fact that
the spatial information and the pairwise relationships between variables may
not be fully exploited.

Constructing the relational graphs can help us better describe, analyze, and
visualize the data, especially in domains such as transportation networks,
social networks, and brain networks. In general, there are 2 main approaches,
statistics-based such as probabilistic models, and deep learning-based, which
often incorporate a graph generating module with a predicting model for
downstream tasks in an end-to-end training framework. GraphLASSO [18] is
popular method in the first approach, it calculates an estimation of the
precision matrix, which can represent the conditional correlations between
each pair of variables. To better model a evolving graph, Time-Varying
Graphical Lasso (TVGL) [23] uses temporal information to adjust the covariance
matrix over time. On the another hand, several works in the second approach,
such as [14] and [12], learn an unique weighted adjacency matrix by computing
the inner product between all pairs of nodes’ the embedding vectors. To
guarantee the necessary sparsity, [14] selects top-k position in each row of
the matrix, which is not flexible and dependent significantly on
hyperparameter selection (i.e. $k$ is a hyperparameter). In a probabilistical
point of view, one often assumes that the learnt dense matrix describes
parameters of a high dimensional Bernoulli distribution of the correlation
between series. LDS [24] tackles the graph learning problem for non-temporal
data by using a bi-level optimization routine and a straight-through gradient
trick which, nonetheless, requires dense computation. The NRI approach [25]
learns a latent variable model predicting the interaction of physical objects
via learning edge attributes of a fully-connected graph. GTS [17] simplifies
the NRI module by considering binary relationships only, and integrates graph
inference in a spatiotemporal recurrent graph neural network [13]. Both NRI
and GTS exploit the categorical Gumbel trick [26, 27] to sample a discretized
adjacency matrix while maintaining differentiability but still require dense
computation. Most of aforementioned methods do not guarantee to produce a
sparse graph, which is more intepretable and explanable in many real world
datasets. Some of other works require pre-defined graph based on expert
knowledge, such as RGSL [28], to regularize sparsity of the structure.
However, the expert knowledge is not avaiable in many cases and may not
provide an accurate explanation for the complex relationship among data.

To address these above existing limitations in the field, our approach
estimates the precision matrix with GraphLASSO, and uses it as parameters of a
Gumbel distribution to sample the graph structure (i.e. adjacency matrix).
Then, this structure is the input to a graph recurrent neural network module
for time series forecasting as a downstream task. The decoupled training
procedure in our proposal not only produces an effective graph structure but
also has less training time in comparison with end-to-end training framework.
Furthermore, we propose an unique loss function aiming to minimize the
forecasting loss and construct the sparse graph structure of data
simultaneously.

##  III Background

In this study, we target on the Multi-variate Time-Series Forecasting (MTSF)
problem. MTS can be denoted as a tensor $X\in\mathbb{R}^{N\times T}$, where
$T$ is the number of timesteps, $N$ denote the number of series. At a time
$t$, given a sequence $H$ time steps of historical observation,
$X=\\{X_{t-H+1},X_{t-H+2},...,X_{t}\\}$, we want to predict the future values
of the next $F$ time steps, $\\{X_{t+1},X_{t+2},...,X_{t+F}\\}$. Beside the
observation, we can concatenate the auxiliary features $S_{t}$ to the input
for improving forecasting performance.

Our approach models the relationships among entities (i.e. variables) by
evolving graphs / networks. Since then, MTS data can be treated as spatial-
temporal data. To accomplish the forecasting tasks, we will use STGNNs models.
Here, we give some formal definition for the problem and algorithms.

###  III-A GraphLASSO

GraphLASSO consider a sequence of multivariate observations in
$\mathbb{R}^{p}$ sampled from a distribution $X\sim N(0,\Sigma(t))$. In static
learning setting, the standard variation $\Sigma(t)$ is unique for all $t$. We
estimate the precision matrix $\Theta=\Sigma^{-1}$ with the constraint that
$\Theta$ is symmetric positive-semidefinite by solving the following
optimization problem:

$\min_{\Theta\in S^{P}}\quad
L(\Theta)=-\log\text{det}(\Theta)+\text{trace}(S\Theta)+\lambda\|\Theta\|_{1}$
(1)

where $S$ is the empirical covariance defined as $\frac{1}{n}XX^{T}$. When $S$
is invertible, $L(\Theta)$ encourages $\Theta$ to converge to $S^{-1}$. When
$\Sigma(t)$ changes during time axis, we have an assumption that they are
fixed in each given time interval $(t_{i},t_{i-1})$. Each time interval has at
least 1 observation and we use these to estimate the corresponding precision
matrix for each interval. [23] provided a general formula of Time-Varying
GraphLASSO (TVGL) problem:

$\min_{\\{\Theta_{t}\\}_{t=1}^{T}}\sum_{i=1}^{T}(-l_{i}(\Theta_{i})+\lambda\|\Theta_{i}\|_{1})+\beta\sum_{j=2}^{T}\psi(\Theta_{j}-\Theta_{j-1})$
(2)

where,
$l_{k}(\Theta_{k})=\log\text{det}(\Theta_{k})-\text{trace}(S_{k}\Theta_{k})$,
$\beta\leq 0$ and $\psi$ is a convex function, which gives penalty for the
significant change between two $\Theta$s, see [23] for more details.

###  III-B Graph Convolution Recurrent Network

Graph Convolutional Recurrent Network (GCRN) [29] is a deep learning model
that can predict structured sequences of data. It is a generalization of
classical recurrent neural networks (RNN) to data structured by an arbitrary
graph. The main idea of GCRN is to merge different representations of the data
provided by Graph Convolutional Network (GCN) layers and by recurrent layers.
GCN layers compute the spatial features of each node through a convolution
operator that integrates information from each node’s local neighbours.
According to [30, 31, 32], the graph convolution operator can be well-
approximated by 1st order Chebyshev polynomial expansion and can be
generalized to high-dimensional GCN as:

$G_{\star W}(X,A)=(I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}})XW+W_{b}.$ (3)

The recurrent layers in GCRN are used to exploit dynamic patterns and capture
the long-term dependencies in the sequential data. Basically, the hidden state
at time step $t$ can be computed by a function of hidden at step $t-1$ and
input at time $t$. There are several variations of recurrent can be applied
here such as Long Short-Term Memoery (LSTM) [20] or Gated Recurrent Unit (GRU)
[33], etc.

Figure 1: Decoupling architecture includes two phases. Phase 1 generates a
sparse estimation of precision matrix. Depending on the case of a static
structure or a time-varying structure, the output of Phase 1 is a matrix or a
matrix sequence corresponding to the time intervals. Phase 2 samples the
sparse graph structures based on these matrices and fit them into graph
forecasting module.

##  IV Method

We establish two phases sequentially that will be trained independently (see
Figure 1). The first one is to estimate the precision matrices by using
GraphLASSO. The second one is to forecast the future time series using graph
convolution recurrent networks.

###  IV-A Decoupled training architecture

Depending on whether there is a structural change in the network, we consider
one of the two types of GraphLASSO problem presented in the section III-A. To
align with Gaussian distribution assumption of GraphLASSO, we first detrend
and normalize entire the time-series. The MTS dataset now is assumped that
follow multivarite Gaussian distribution $\mathcal{N}(0,\Sigma)$. The
precision matrix $\Theta=\Sigma^{-1}$ represent the conditional dependence
between any two pair of series. It means that $\Theta_{ij}=0$ if and only if
series $i$ is independent with series $j$ given other variables.

In case of without assuming about the time-varying structures, we can use
standard interior-point methods to solve (1). We learn only one precision
matrix $\Theta$ of entire dataset. In case of dynamically inferring graph
structures, we can use Time-Varying Graphical Lasso (TVGL) algorithm to solve
problem (2) (see [23]). This algorithm is based on the Alternating Direction
Method of Multipliers (ADMM) [34], a well-known approach for convex
optimization that separates the problem into pieces that are easier to solve.
The result is a sequence of $T$ precision matrices
$\\{\Theta_{t}\\}_{t=1}^{T}$ corresponding to $T$ consecutive time intervals.
We choose $\psi(\Theta)=\|\Theta\|_{F}^{2}=\sum_{i,j}\Theta_{ij}^{2}$ (i.e.
$\|\cdot\|_{F}$ denotes the Frobenius norm) which can make smooth transitions
of graph structure during time axis [23].

Once the precision matrix $\Theta$ is estimated, we fit it into the graph
forecasting module and sample the graph adjacency matrix $A$. We also consider
$A_{t}$ to be a random variable with a parameterized distribution $P$ with the
parameter $\Theta_{k}$, where $t$ is in time interval $(t_{i},t_{i-1})$. Then
the loss function of forecasting task will be rewritten as follows:

$\min\mathbb{E}_{A_{t}\sim
P(\Theta)}\sum_{t}\textit{L}(f(X_{t+1:t+T},A_{t}),X_{t+T+1:t+T+F})$ (4)

In practice, we estimate the precision matrix on the rescaled dataset, so that
$0\leq\Theta_{ij}\leq 1$. We directly sample edge from node $i$ to $j$ by from
the Bernoulli distribution ($\Theta_{ij}$), that means $A_{ij}=1$ with
probability $\Theta_{ij}$. Using a sampling technique instead of using a fixed
threshold can help the model more flexible. Moreover, iterative sampling of
the graph structures during training process can be seen as a dropout layer,
which can decrease overfitting [28].

###  IV-B Forecasting with graph neural network

In this part, we present the details of our forecasting module which is
inspired by the GCRN architecture [12, 2]. The module contains a Graph
Convolution Module (GCM) which exploits spatial information and a Gate
Recurrent Unit (GRU) network to capture temporal dependencies. The output of
GCM with input $X$ and the corresponding adjacency matrix $A$ is denoted by
$Z(X,A)$ as:

$Z(X,A)=G_{\star W}(X,A)=(I+D^{-\frac{1}{2}}AD^{-\frac{1}{2}})XW+W_{b},$ (5)

where $I$ is identity matrix, $D$ is degree matrix of the graph, $W$ and
$W_{b}$ are the learnable parameters. The updating mechanism of hidden state
$h_{t}$ is given as follows:

$\begin{split}r_{t}&=\sigma(G_{r}(X_{t}\mathbin{\|}h_{t-1},A),\\\
z_{t}&=\sigma(G_{z}(X_{t}\mathbin{\|}h_{t-1},A),\\\
c_{t}&=\text{tanh}(G_{c}(X_{t}\mathbin{\|}(r_{t}\odot h_{t-1}))),\\\
h_{t}&=z_{t}\odot h_{t-1}+(1-z_{t})\odot c_{t}.\end{split}$ (6)

Assume that we generate a sequence of graph adjacency matrices
$\\{A^{k}:k\in\\{1,2,\ldots,T\\}\\}$ corresponding to structures at $T$ time
intervals. So, there are more than one graphs appear in a input-output
sequence pair at some $t$. We tackle this issue by a simple aggregation
strategy. Now, the graph convolution operator can be rewrited in the form:

$G_{TV}(X)=\frac{1}{n}\sum_{A\in A^{n}}G(X,A),$ (7)

where $A^{n}$ is the set of $n$ graph structures that intersect with data
sequence $X_{t-H+1:t+F}$. Finally, we use Mean Absolute Error (MAE) as the
loss function of forecasting module:

$L_{f}(\bm{W})=\frac{1}{F}\|X_{t+1:t+F}-\hat{X}_{t+1:t+F}\|_{1},$ (8)

where $\bm{W}$ denotes the set of all learnable parameters.

##  V Experiments

In this section, extensive experiments will be conducted to prove the efficacy
and effectiveness of our proposed method, and the corresponding analysis will
also be provided. Details about datasets and code can be found at
https://github.com/HySonLab/GraphLASSO.

###  V-A Datasets

We report results on two types of MTSF: (i) Small number of series including
CA1, CA2, TX1, TX2, WI1, and WI2; and (ii) Large number of series including
CA1-Food1, CA2-Food1, Electricity100, and Solar-energy-1h. Detail about these
datasets can be founded in Appendix A-A.

###  V-B Experiments setting

We compare the performance of our method with baseline methods as following:
(i) HA: Historical Average, (ii) VAR: Vector Auto-Regressive model [6], (iii)
LSTM: Long Short-Term Memory (LSTM) [20], (iv) GTS: An end-to-end framework
for graph structure learning and GCRN forecasting [17] and (v)
$\text{GCRN}^{*}$: a GCRN model [29] that has the same architecture with our
model but the input graph structure is fully connected graph (i.e. complete
graph). We denote our model as $\text{Our}^{s}$ and $\text{Our}^{d}$, in which
$\text{Our}^{s}$ includes a non-dynamic (i.e. static) graph learning module
and $\text{Our}^{d}$ has a dynamic graph learning using time-varying
GraphLASSO module. The architecture of the corresponding GCRN forecasting
modules were presented in IV-B.

In the benchmark on M5 datasets, we use 14 historical timesteps to predict 7
future timesteps. In electricity and solar-energy, we use 168 input steps and
predict 1 output step. Details about other parameters and training for all
models can be founded in Appendix A-D.

###  V-C Results

| HA | VAR | LSTM | GTS | GCRN | Ours
---|---|---|---|---|---|---
CA1 | 0.1180 | 0.1063 | 0.1026 | 0.0933 | 0.0956 | 0.0902
CA2 | 0.1311 | 0.1277 | 0.0692 | 0.0698 | 0.0713 | 0.0691
TX1 | 0.1038 | 0.0620 | 0.0659 | 0.0665 | 0.0671 | 0.0616
TX2 | 0.0854 | 0.0625 | 0.0619 | 0.0623 | 0.0684 | 0.0621
WI1 | 0.1608 | 0.0768 | 0.0787 | 0.0723 | 0.0730 | 0.0715
WI2 | 0.1488 | 0.0881 | 0.0732 | 0.0801 | 0.0843 | 0.0762

TABLE I: Mean absolute error of M5 store-level datatset

| HA | VAR | LSTM | GTS | GCRN | $\text{Our}^{s}$ | $\text{Our}^{d}$
---|---|---|---|---|---|---|---
| MAE | MAE | MAE | MAE | t-time | MAE | t-time | MAE | t-time | MAE | t-time
CA1_Food1 | 1.9869 | 2.2169 | 1.4045 | 1.3601 | 70s | 1.5108 | 46s | 1.4019 | 37s | 1.3916 | 40s
TX1_Food1 | 1.7883 | 2.3486 | 1.6067 | 1.3741 | 75s | 1.4061 | 43s | 1.3629 | 34s | 1.3777 | 38s
WI1_Food1 | 1.8058 | 2.3056 | 1.4335 | 1.0818 | 72s | 1.1238 | 45s | 1.0843 | 35s | 1.0964 | 38s
Electricity100 | 231.3055 | 101.2598 | 91.6094 | 70.4458 | 390s | 76.5002 | 310s | 75.8671 | 215s | 72.6113 | 230s
Solar-energy-1h | 8.2004 | 4.4045 | 2.5334 | 1.8263 | 125s | 2.0925 | 95s | 1.8406 | 75s | 1.1802 | 80s

TABLE II: Forecasting error and time training evaluated on 5 MTS datasets
which have more than 100 time series.

We first report the MAE results of the next 7th step prediction on datasets
that have small number of series in Table I. The historical average method has
the most misleading prediction among the models. Deep learning models have
better results than VAR on all 6 datasets, which indicates that representing
inter-relationships using non-linear models will be more effective than the
linear one. LSTM gives the best performance in 2 stores, TX2 and WI2, while
our model is best in the remaining 4 stores. Among the models using GNNs, our
model gives better results than GTS and $\text{GCRN}^{*}$. This shows that the
graph produced by GraphLasso can provide more effective spatial information
for the GNN forecasting model. Moreover, our model can control the sparsity of
the generated graphs and it can give a reasonable explanation for the
relationship between series.

We report the results of comparing our models and baselines on MTS datasets
with a larger number of series as mentioned above including MAE error and
training time denoted t-time in Table II. Our methods have competitive results
against baselines. Model $\text{Our}^{s}$, which uses static GraphLASSO, has
best MAE loss in $\text{TX1}\\_\text{Food1}$ dataset and our dynamic model
archives the best result in Solar-energy-1h dataset. In particular, there was
a shift distribution in Solar-energy dataset because it depends on sunlight
intensity during the year. So it may contain a complex time-varying
correlation and the relationships between nodes change over time. Model
$\text{Our}^{d}$, which uses time-varying GraphLASSO for graph learning, can
capture this change well and makes better predictions than the others.

In the remaining 3 datasets, our models have results that are not too much
worse than the best baseline but significantly reduce the training time. Our
model computes the relationship graph based on the statistical estimation of
the precision matrix, thus reducing both the computational complexity and the
number of parameters to be learned by the graph learning module compared to an
end-to-end approach. Another advantage of our models is the use of sparse
graphs, which can be implemented in the form of sparse matrices, thus reducing
the computation time of graph convolution operators. In summary, our method
can capture accurately both spatial information and temporal information in
the correlated MTS data and significantly reduce the training time in
comparison with other GNN-based methods.

##  VI Conclusion

In this paper, we have proposed a graph forecasting method consisting of two
separate training modules to solve the MTSF problem. The first module learns
the sparse graph structures that describe the correlation between series in a
data-driven manner. We use GraphLASSO to directly estimate the precision
matrix of data to construct the graph structures and control their sparsity.
We consider two types of learning graph structures: the relationship between
series is unchanged / static (i.e. independent of time) or dependent on each
specific time interval (i.e. time-varying). In the static case, only one graph
is generated. In the time-varying case, we construct a sequence of graphs that
describe the pairwise relationships between series. The second module is built
from GCRN that takes both the data and the graph structures generated from the
first module as input to learn a forecasting model. Taking advantage of the
sparsity, we can reduce the training time of both the graph generating module
and the time-series forecasting module. Our experiments on three datasets can
be divided into two cases: small number of series and large number of series.
In both cases, our models not only show a competitive performance but also
reduce training time by $\sim 40\%$ in comparison with other forecasting
baselines. In future works, we aim to incorporate the end-to-end learning and
experts’ knowledge graphs to making more meaningful explanation and improve
forecasting performance.

## References

  * [1] Z. Cui, K. Henrickson, R. Ke, and Y. Wang, “Traffic graph convolutional recurrent neural network: A deep learning framework for network-scale traffic learning and forecasting,” IEEE Transactions on Intelligent Transportation Systems, vol. 21, no. 11, pp. 4883–4894, 2019.
  * [2] Y. Li, R. Yu, C. Shahabi, and Y. Liu, “Diffusion convolutional recurrent neural network: Data-driven traffic forecasting,” in 6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings, OpenReview.net, 2018.
  * [3] T. S. Hy, V. B. Nguyen, L. Tran-Thanh, and R. Kondor, “Temporal multiresolution graph neural networks for epidemic prediction,” in Proceedings of the 1st Workshop on Healthcare AI and COVID-19, ICML 2022 (P. Xu, T. Zhu, P. Zhu, D. A. Clifton, D. Belgrave, and Y. Zhang, eds.), vol. 184 of Proceedings of Machine Learning Research, pp. 21–32, PMLR, 22 Jul 2022.
  * [4] S. Liao, J. Yin, and W. Rao, “Towards accurate retail demand forecasting using deep neural networks,” in Database Systems for Advanced Applications (Y. Nah, B. Cui, S.-W. Lee, J. X. Yu, Y.-S. Moon, and S. E. Whang, eds.), (Cham), pp. 711–723, Springer International Publishing, 2020.
  * [5] A. Patton, “Chapter 16 - copula methods for forecasting multivariate time series,” in Handbook of Economic Forecasting (G. Elliott and A. Timmermann, eds.), vol. 2 of Handbook of Economic Forecasting, pp. 899–960, Elsevier, 2013.
  * [6] J. Stock and M. Watson, “Vector autoregressions,” Journal of Economic Perspectives, vol. 15, no. 4, p. 101 – 116, 2001.
  * [7] H. Spliid, “A fast estimation method for the vector autoregressive moving average model with exogenous variables,” Journal of the American Statistical Association, vol. 78, no. 384, pp. 843–849, 1983.
  * [8] S. Li, X. Jin, Y. Xuan, X. Zhou, W. Chen, Y.-X. Wang, and X. Yan, “Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting,” Advances in neural information processing systems, vol. 32, 2019.
  * [9] Y. Zhang and J. Yan, “Crossformer: Transformer utilizing cross-dimension dependency for multivariate time series forecasting,” in The Eleventh International Conference on Learning Representations, 2022.
  * [10] B. Yu, H. Yin, and Z. Zhu, “Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting,” in Proceedings of the 27th International Joint Conference on Artificial Intelligence, IJCAI’18, p. 3634–3640, AAAI Press, 2018.
  * [11] D. Cao, Y. Wang, J. Duan, C. Zhang, X. Zhu, C. Huang, Y. Tong, B. Xu, J. Bai, J. Tong, and Q. Zhang, “Spectral temporal graph neural network for multivariate time-series forecasting,” in Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS’20, (Red Hook, NY, USA), Curran Associates Inc., 2020.
  * [12] L. Bai, L. Yao, C. Li, X. Wang, and C. Wang, “Adaptive graph convolutional recurrent network for traffic forecasting,” in Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS’20, (Red Hook, NY, USA), Curran Associates Inc., 2020.
  * [13] Y. Li, R. Yu, C. Shahabi, and Y. Liu, “Diffusion convolutional recurrent neural network: Data-driven traffic forecasting,” arXiv preprint arXiv:1707.01926, 2017.
  * [14] Z. Wu, S. Pan, G. Long, J. Jiang, X. Chang, and C. Zhang, “Connecting the dots: Multivariate time series forecasting with graph neural networks,” in Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD ’20, (New York, NY, USA), p. 753–763, Association for Computing Machinery, 2020.
  * [15] F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini, “The graph neural network model,” IEEE Transactions on Neural Networks, vol. 20, no. 1, pp. 61–80, 2009.
  * [16] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in International Conference on Learning Representations, 2017.
  * [17] C. Shang and J. Chen, “Discrete graph structure learning for forecasting multiple time series,” in Proceedings of International Conference on Learning Representations, 2021.
  * [18] J. Friedman, T. Hastie, and R. Tibshirani, “Sparse inverse covariance estimation with the graphical lasso,” Biostatistics, vol. 9, pp. 432–441, 12 2007.
  * [19] M. Yuan and Y. Lin, “Model selection and estimation in the Gaussian graphical model,” Biometrika, vol. 94, pp. 19–35, 03 2007.
  * [20] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural Comput., vol. 9, p. 1735–1780, nov 1997.
  * [21] G. Lai, W.-C. Chang, Y. Yang, and H. Liu, “Modeling long-and short-term temporal patterns with deep neural networks,” in The 41st international ACM SIGIR conference on research & development in information retrieval, pp. 95–104, 2018.
  * [22] S.-Y. Shih, F.-K. Sun, and H.-y. Lee, “Temporal pattern attention for multivariate time series forecasting,” Machine Learning, vol. 108, pp. 1421–1441, 2019.
  * [23] D. Hallac, Y. Park, S. Boyd, and J. Leskovec, “Network inference via the time-varying graphical lasso,” in Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 205–213, 2017.
  * [24] L. Franceschi, M. Niepert, M. Pontil, and X. He, “Learning discrete structures for graph neural networks,” in Proceedings of the 36th International Conference on Machine Learning (K. Chaudhuri and R. Salakhutdinov, eds.), vol. 97 of Proceedings of Machine Learning Research, pp. 1972–1982, PMLR, 09–15 Jun 2019.
  * [25] T. Kipf, E. Fetaya, K.-C. Wang, M. Welling, and R. Zemel, “Neural relational inference for interacting systems,” in Proceedings of the 35th International Conference on Machine Learning (J. Dy and A. Krause, eds.), vol. 80 of Proceedings of Machine Learning Research, pp. 2688–2697, PMLR, 10–15 Jul 2018.
  * [26] E. Jang, S. Gu, and B. Poole, “Categorical reparameterization with gumbel-softmax,” in 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings, OpenReview.net, 2017.
  * [27] C. J. Maddison, A. Mnih, and Y. W. Teh, “The concrete distribution: A continuous relaxation of discrete random variables,” in 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings, OpenReview.net, 2017\.
  * [28] H. Yu, T. Li, W. Yu, J. Li, Y. Huang, L. Wang, and A. Liu, “Regularized graph structure learning with semantic knowledge for multi-variates time-series forecasting,” in Proceedings of the Thirty-First International Joint Conference on Artificial Intelligence, IJCAI-22 (L. D. Raedt, ed.), pp. 2362–2368, International Joint Conferences on Artificial Intelligence Organization, 7 2022.  Main Track.
  * [29] Y. Seo, M. Defferrard, P. Vandergheynst, and X. Bresson, “Structured sequence modeling with graph convolutional recurrent networks,” in Neural Information Processing: 25th International Conference, ICONIP 2018, Siem Reap, Cambodia, December 13-16, 2018, Proceedings, Part I 25, pp. 362–373, Springer, 2018.
  * [30] J. Bruna, W. Zaremba, A. D. Szlam, and Y. LeCun, “Spectral networks and locally connected networks on graphs,” CoRR, vol. abs/1312.6203, 2013.
  * [31] M. Henaff, J. Bruna, and Y. LeCun, “Deep convolutional networks on graph-structured data,” arXiv preprint arXiv:1506.05163, 2015.
  * [32] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in Proceedings of the 30th International Conference on Neural Information Processing Systems, NIPS’16, (Red Hook, NY, USA), p. 3844–3852, Curran Associates Inc., 2016.
  * [33] K. Cho, B. van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using RNN encoder–decoder for statistical machine translation,” in Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), (Doha, Qatar), pp. 1724–1734, Association for Computational Linguistics, Oct. 2014.
  * [34] S. P. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach. Learn., vol. 3, no. 1, pp. 1–122, 2011\.
  * [35] S. Makridakis, E. Spiliotis, and V. Assimakopoulos, “M5 accuracy competition: Results, findings, and conclusions,” International Journal of Forecasting, vol. 38, no. 4, pp. 1346–1364, 2022.
  * [36] F. Schaipp, O. Vlasovets, and C. L. Müller, “Gglasso - a python package for general graphical lasso computation,” Journal of Open Source Software, vol. 6, no. 68, p. 3865, 2021.
  * [37] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.

##  Appendix A Dataset and model parameter

###  A-A Dataset

We used M5 retail [35] and its sub-dataset in our experiments. The M5 is a
retail sales forecasting dataset that uses real-life, hierarchically
structured sales data with intermittent characteristics. It collects retail
data at 10 stores in 3 states in the United States including California (CA),
Texas (TX) and Wisconsin (WI). Each store contains 3,049 items and they are
divided into 7 departments (e.g., Food1, Food2, Household1, etc.). All the
series were daily recorded in a period of 1,943 days. The information about
items and departments is encoded, hence we cannot provide an expert knowledge
graph for the model. In this paper, we used 6 store-level series as subsets of
original M5, which are CA1, CA2, TX1, TX2, WI1, and WI2. This store-level data
is then be aggregated by its department, so the dataset has 7 series and 1,943
timesteps. We created item-level subset grouped by the department in each
store because the items in same department may have important dependencies.
Data in some departments contains many zero values which are related to
disruption in the selling process of items. To reduce data processing and give
balance comparision with baselines, we only report result on Food1 department
which has 216 items.

We also evaluated our model on other commonly used MTS datasets which are:
Electricity, Solar-energy [14]. The origin Electricity dataset contains 321
series represent the electricity consumption of clients from 2012 to 2014 and
the length of series is 26,304. We only used first 100 series in this dataset
in our experiments. Solar-energy dataset includes solar power output collected
from 137 sensors in Alabama state in 2007. The timestep in Solar-energy is
aggregated to 1 hour, and the number of samples is 8,760. We denoted the two
new datasets are Electricity100 and Solar-energy-1h.

We present some statstic detail about these datasets in III.

Dataset | Mean | Std | Data range
---|---|---|---
Electricity100 | 650.14 | 283.07 | 0 $-$ 15443
Solar-energy-1h | 6.35 | 9.10 | 0 $-$ 86.85
CA1_Food1 | 1.52 | 2.86 | 0 $-$ 92
TX1_Food1 | 1.06 | 2.42 | 0 $-$ 108
WI1_Food1 | 1.44 | 2.83 | 0 $-$ 166
CA1 | 634.28 | 656.17 | 0 $-$ 4042
TX1 | 453.41 | 452.41 | 0 $-$ 2919
WI1 | 511.36 | 493.51 | 0 $-$ 3203

TABLE III: Statistics of MTS datasets used in experiments.

###  A-B Data preprocessing

Because the M5 retail data at item-level is intermittent and erratic, we first
used $\log(1+x)$ transformation on input sale data. For the department-level
M5 data, we used min-max scaling transformation in each department to keep the
data range in (0,1). We applied normalizing transformation for all datasets
before estimating the precision matrix:

$\hat{X}=\dfrac{X-\text{mean}(X)}{\text{std}(X)}$ (9)

We also used day-in-week and hour-in-day as auxiliary temporal data to improve
forecasting performance. We generated sequence-to-sequence type data for deep
learning models. The train/validate/test ratio was separated as ratio
$70\%/10\%/20\%$.

###  A-C Evaluate metrics

We used Mean Absolute Error (MAE) to measure forecasting performance. We
evaluated forecasting results at $F-th$ timestep in the future, we have the
MAE error between ground truth and the prediction of all series as:

$\dfrac{1}{D}\sum_{j=1}^{D}|X_{t+F,j}-\hat{X}_{t+F,j}|$ (10)

where $D$ is the number of series. We also reported the training time on
larger number of series datasets. Total training time including Phase 1 and
Phase 2 of our models were divided by the number of epochs. We compare them
with the average training time of each epoch of other baselines. The unit
measure is second. All the experiments were performed on a same computer with
GPU NVIDIA GTX 1650Ti.

###  A-D Implementation details

The baselines have the detail of implementation as follow:

  * •

HA: We use the historical timestamp value of the past 24 steps in
Electricity100 and Solar-energy-1h datasets to predict 6 steps in the future
and the past 14 steps in the past and predict 7 future steps in the remaining
datasets. All the predictions in the horizon are invariant.

  * •

VAR: We implement the VAR model based on the python package statsmodel. The
lag parameter was selected as 7 for all M5 sub-datasets and 12 for the
remaining.

  * •

LSTM: We implement an encoder-decoder model which two layers of LSTM with
Tensorflow. Each layer has 32 and 64 hidden units for M5 datasets and the
remaining, respectively. A fully-connected layer was added to predict the
output.

  * •

$\text{GCRN}^{*}$: We set the maximum step of randoms walks to 2, the initial
learning rate to 0.01, drop out rate to 0.1 for all datasets, the GRU layer
has 32 and 64 hidden units for M5 datasets and the remaining, respectively.

  * •

GTS: We used the same architecture for the graph learning module with 2 layers
of CNN and 2 layers of MLP as in [17] in the M5 sub-datasets. In the other
datasets, we used 3 layers of CNN and 2 MLP layers. The batch size was
selected by 32 for all the datasets. Learning rate and hidden units were
chosen the same with the $\text{GCRN}^{*}$ model. Max diffusion step was
selected by 2 and the number of RNN layers was 1.

In our two models $\text{Our}^{s}$ and $\text{Our}^{d}$, the precision
matrices were estimated by using the GGLASSO package [36]. The $\lambda$ and
$\beta$ parameters in Equation 2 were 0.1 and 0.05, respectively. For all the
deep learning models including ours, we optimize with the Adam optimizer [37]
for 100 epochs and use an early stopping strategy with the patience of 15
epochs by monitoring the loss in the validation set. After hyperparameter
searching, we chose the initial learning rate as 0.005 for all datasets. Other
parameters were the same with the $\text{GCRN}^{*}$ model [29].


#  $X(3872)$ and $T_{cc}$: structures and productions in heavy ion collisions

Hyeongock Yun<EMAIL_ADDRESS>Department of Physics and Institute of
Physics and Applied Physics, Yonsei University, Seoul 03722, Korea  Daeho Park
<EMAIL_ADDRESS>Department of Physics and Institute of Physics and
Applied Physics, Yonsei University, Seoul 03722, Korea  Sungsik Noh
<EMAIL_ADDRESS>Department of Physics and Institute of Physics and
Applied Physics, Yonsei University, Seoul 03722, Korea  Aaron Park
<EMAIL_ADDRESS>Department of Physics and Institute of Physics and
Applied Physics, Yonsei University, Seoul 03722, Korea  Woosung Park
Department of Physics and Institute of Physics and Applied Physics, Yonsei
University, Seoul 03722, Korea  Sungtae Cho<EMAIL_ADDRESS>Division of Science Education, Kangwon National University, Chuncheon 24341,
Korea  Juhee Hong  Department of Physics and Institute of Physics and Applied
Physics, Yonsei University, Seoul 03722, Korea  Yongsun Kim<EMAIL_ADDRESS>Department of Physics, Sejong University, Seoul, Korea  Sanghoon Lim
<EMAIL_ADDRESS>Department of Physics, Pusan National University, Pusan,
Korea  Su Houng Lee<EMAIL_ADDRESS>Department of Physics and Institute
of Physics and Applied Physics, Yonsei University, Seoul 03722, Korea

###### Abstract

We argue why the recently observed $T_{cc}$ could either be a compact
multiquark configuration or a loosely bound molecular configuration composed
of charmed mesons, whereas the $X(3872)$ is most likely a molecular
configuration. The argument is based on different short range interactions for
these tetraquark states coming from the color-color and color-spin interaction
in a quark model, and the presence of a common strong D-wave mixing at larger
distance similar to the deuteron case, which for the molecular configurations
lead to large sizes. Such an analogy at large distance allows us to calculate
the transverse momentum dependence of the loosely bound molecular
configuration of tetraquarks produced in heavy ion collisions using the
coalescence model that successfully reproduces the deutron data using the
proton spectra. The ratio of the integrated $X(3872)$ yield obtained from our
method to the $\psi(2S)$ yield obtained from statistical hadronization model
method is calculated to be $0.806\pm 0.234$, which is a factor of 2.47 larger
than that obtained by using statistical model predictions for both particles
and in line with the data from the CMS experiment. As the previously
calculated transverse momentum distribution of the $T_{cc}$ assuming the
structure to be a compact multiquark configuration is markedly different,
experimental measurements of the transverse distribution of the tetraquark
states will discriminate between their two possible structures.

##  I Introduction

Besides the historical discovery of Higgs boson, one of the greatest victories
of Large Hadron Collider (LHC) is the discovery of numerous new hadrons, more
than 50 over the past decade CERN-article2021 . In particular, observation of
several resonances composed of four or five quarks (and anti-quarks)
unambiguously confirmed the existence of tetraqurk and pentaquark states, so
called exotic particles. These results provide novel inputs for understanding
the nature of strong force in the hadronic scale, thus providing crucial clue
to comprehend the property of confinementBelle:2003nnu ; LHCb:2015yax ;
LHCb:2019kea ; D0:2016mwd ; LHCb:2021vvq . The most important remaining
question would be the internal structure; are they in molecular system
composed of loosely bound hadrons, or are they tight multiquark states where
the quarks and anti-quarks form compact configurations with the size of a
typical hadron? No clear answer has been reported so far.

In this paper, we discuss $X(3872)$ ($c\bar{c}u\bar{u}$, a hidden-charm
tetraquark) and $T_{cc}$ ($cc\bar{u}\bar{d}$, open-charm tetraquark), which
are of extraordinary interest for their conjugate relation in terms of charm
quark flavor. Furthermore, they are the representative candidates that could
either be compact multiquark configurations or loosely bound molecular
configurations. We provide strong arguments that $X(3872)$ must be in
molecular state of $D\bar{D}^{*}$ based on two evidences; (a) one is the
absence of strong short range attraction required for a compact configuration,
and (b) the other is the presence of a strong D-wave mixing coming from pion
exchange, a mechanism crucial for the binding of the deuteron. On the other
hand, we report that the $T_{cc}$ has possibilities to be either a compact
multiquark state or a loosely bound molecular state of $DD^{*}$ due to two
competing effects: (a) the presence of a strong short range attraction, and
(b) a strong mixing of the D-wave configuration of $DD^{*}$ by the
$\pi$-exchange that is similar but weaker than that present in $X(3872)$.

We also show that the transverse momentum distributions calculated for both
the compact mutiquark and loosely bound molecular states are completely
different, thus allowing the discrimination of their structures through high
statistics measurements at the LHC heavy ion experiments in near future. Such
measurement will need data at lower transverse momentum than that reported in
the recent anomalous enhancement of $X(3872)$ production in heavy ion
collision experimentsCMS:2021znk .

The paper is organized as follows. In Sec. II, we discuss the short range
interaction of tetraquark states based on a quark model. In Sec. III, we
discuss the D-wave mixing. In Sec. IV, we discuss the transverse momentum
distribution for the deuteron and 3He measured in heavy ion collision and the
fit from a coalescence model. In Sec. V, we use the same coalescence model to
predict the distribution for the molecular configuration of tetraquarks and
compare the result assuming a compact configuration. Finally, we give the
summary in Sec. VI. The appendix provides the fitting formula and details of
the coalescence model used in this work.

##  II Short distance attractions

In a previous publicationPark:2019bsz , we have shown that the short distance
parts of the baryon-baryon interactions for various quantum numbers from the
recent lattice calculation can be well reproduced using a constituent quark
model with color-spin interaction. In fact, the full constituent quark model
calculations are found to be comparable to the simple estimates based on the
following color spin factor for a multiquark configurationPark:2019bsz .

$\displaystyle K$
$\displaystyle=-\sum_{i<j}^{n}\frac{1}{m_{i}m_{j}}\lambda^{c}_{i}\lambda^{c}_{j}\sigma_{i}\cdot\sigma_{j},$
(1)

where $\lambda_{i}^{c}$ and $\sigma_{i}$ are, respectively, the color and spin
operators and $m_{i}$ the mass of the $i$’th quark among $n$ quarks. Although
a full quark model calculation will include the spatial wave function, if the
interquark distancees among quarks are all the same, the short distance
interaction between hadrons can be well described using only the color-spin
factor times an overall factor that determines the interaction strength at the
hadron size.

The color-spin wave function for $X(3872)$, which is a
$I^{G}(J^{PC})=0^{+}(1^{++})$ state, is composed of two states. Using the
quark-antiquark basis $(c\bar{c})\otimes(q\bar{q})$ with $c,q$ being the charm
and light quarks, the two states are $(1\otimes 1)_{C=1}(V\otimes V)_{J=1}$
and $(8\otimes 8)_{C=1}(V\otimes V)_{J=1}$, where the first brackets denote
two color singlets $(1)$ or octets $(8)$ combined into color singlet and the
second two spin 1 $(V)$ into spin 1. After subtracting the $K$ factors for the
lowest $D\bar{D}^{*}$ threshold, the color-spin interaction factor is as
follows.

$K_{X(3872)}-K_{D}-K_{\bar{D}^{*}}=$

$\displaystyle\left(\begin{array}[]{cc}\frac{16}{3m_{c}^{2}}+\frac{16}{3m_{q}^{2}}+\frac{32}{3m_{q}m_{c}}&0\\\
0&-\frac{2}{3m_{c}^{2}}-\frac{2}{3m_{q}^{2}}-\frac{4}{3m_{q}m_{c}}\end{array}\right),$
(4)

where the factor $\frac{32}{3m_{q}m_{c}}$ are due to the subtracted threshold.
One notes that the second diagonal component coming from the color octet
element gives an attractive contribution. However, noting that the $K$ factor
coming from the delta nucleon mass difference is $\frac{16}{m_{q}^{2}}$, which
gives a mass difference of around 290 MeV, we find that using a contitutent
quark mass of $m_{c}=1500$ MeV and $m_{q}=300$ MeV, the attraction can only be
around 17.4 MeV assuming that all quarks occupy the size of a nucleon. The
additional kinetic energy needed to bring the two mesons to 0.3 fm separation,
which is the expected typical size of a compact heavy tetraquark
configuration, is larger than 200 MeV. Hence, the attraction is not strong
enough to confine the $X(3872)$ into a compact multiquark configuration.

A detailed quark model calculation introduces some mixing from the singlet
contribution due to the different spatial distributions of heavy and light
quarks. This could introduce some contribution from the color-color type of
interaction and mixing from the color singlet contribution. However, it should
be noted that at very short distance the color-color interaction is dominated
by the Coulomb potential so that it is repulsive for the color octet heavy
quark-antiquark pair and do not introduce any additional attraction. One can
also convert the previous basis into basis involving $D\bar{D}^{*}$ color
configurations. Then the color octet basis is composed dominantly of color
singlet $D$ and color singlet $\bar{D}^{*}$ basis: the singlet-singlet to
octet-octet fraction is 8 to 1. This result furthermore implies that the
attraction between $D\bar{D}^{*}$ is small at short distance.

On the other hand, assuming that $J^{P}=1^{+}$ for the $T_{cc}$, a better way
to describe the color-spin wave function is using the diquark-antidiquark
$(ud)\otimes(\bar{c}\bar{c})$ basis, for which the two states are
$(\bar{3}\otimes 3)_{C=1}(S\otimes V)_{J=1}$ and
$(6\otimes\bar{6})_{C=1}(V\otimes S)_{J=1}$, where the first and second
brackets show the color combination and spin with $S$ denoting the spin zero
scalar state, respectively. Using these basis, the $K$ factor for the $T_{cc}$
after subtracting the $K$ factor for the threshold is given as follows.

$K_{T_{cc}}-K_{D}-K_{D^{*}}=$

$\displaystyle\left(\begin{array}[]{cc}-\frac{8}{m_{q}^{2}}+\frac{8}{3m_{c}^{2}}+\frac{32}{3m_{q}m_{c}}&\frac{8\sqrt{2}}{m_{c}m_{q}}\\\
\frac{8\sqrt{2}}{m_{c}m_{q}}&-\frac{4}{3m_{q}^{2}}+\frac{4}{m_{c}^{2}}+\frac{32}{3m_{q}m_{c}}\end{array}\right).$
(7)

|
---|---
Figure 1: Left - Tensor potential due to pion exchange contributing to the
$X(3872)$ and $T_{cc}$ with additional short range potential for $T_{cc}$ as a
function of distance between two constituents. Right - Wave functions of
$X(3872)$ and $T_{cc}$. The wave functions of the deuteron are similar to
those of the $T_{cc}$.

Using the same parameters as above, we find the attraction from the first
diagonal component to be around 104.4 MeV, much larger compared to 17.4 MeV
for the $X(3872)$. Diagonalizing the matrix, the attraction is expected to be
larger. Futhermore, the color-color interaction is expected to produce further
attraction in the color triplet or antitriplet quark pair. This is so because
$\lambda^{c}_{i}\lambda^{c}_{j}=-\frac{16}{3},-\frac{8}{3}$ and $\frac{2}{3}$
for the $i,j$ quark or antiquark pair in color singlet, triplet(or
antitriplet) and octet quark pairs, respectively. This attraction, coming from
the most attractive diquark, was the original motivation for the possible
compact $T_{cc}$ configurationZouzou:1986qh . Detailed constituent quark model
calculations agree on the deeply bound compact configurations for the
$T_{bb}$, but do not for the $T_{cc}$ Park:2013fda ; Park:2018wjk ;
Noh:2021lqs .

At the same time, by analyzing the energy of a compact $T_{cc}$ configuration,
one can estimate the attraction between $D-D^{*}$ at relative distance
estimated within that configuration. Using a simple Gaussian wave function,
one finds that at $r_{D-D^{*}}=0.3$ fm, the attraction is around 200 MeV. This
is similar to what has been estimated in the lattice gauge
calculationIkeda:2013vwa . Therefore, we expect that in the hadronic picture,
while there is no short range attraction in the $X(3872)$ channel, there will
be an additional short distance attraction in the $T_{cc}$ channel.

Molecule | Constituents | $V_{\text{short}}$ |  $V_{\pi}$ ($S$ wave) |  $V_{\pi}$ ($D$ wave-mixing) |  $\langle r\rangle$ (fm) |  $E_{b}$ (MeV)
---|---|---|---|---|---|---
$d$ | $pn$ | repulsive | attractive | attractive | $1.9$ | 2.2
$X(3872)$ | $D\bar{D^{*}}$ | negligible | negligible | attractive | 3.0 | 0.33
$T_{cc}$ | $DD^{*}$ | attractive | negligible | attractive | 2.2 | 0.65
Table 1: Comparison of molecular structures.

##  III D-wave mixing

To estimate the D-wave mixing, We will use the pion exchange type of potential
between $D-\bar{D^{*}}$ and $D-D^{*}$ for the $X(3872)$ and $T_{cc}$ quantum
numbers, respectively, that were used in Ref.Tornqvist:1993ng ;
Tornqvist:2004qy .

$\displaystyle V(r)=-\gamma V_{0}\Bigg{[}\begin{pmatrix}0&-\sqrt{2}\\\
-\sqrt{2}&1\end{pmatrix}T_{\pi}(r)\Bigg{]},$ (8)

where $V_{0}=1.3$ MeV as determined by the $\pi N$ coupling constant and the
pion decay constant, and $\gamma=3$ and $-3$ for $X(3872)$ and and $T_{cc}$
channels, respectivelyTornqvist:1993ng . For the $X(3872)$, the attractive
tensor potential and the D-wave mixing coming from the off diagonal component
in the tensor potential lead to a loosely bound molecular configuration
Tornqvist:1993ng . While there is a single pion exchange in Tornqvist:1993ng ,
this contribution can be neglected. On the other hand, because the tensor part
is repulsive for the $T_{cc}$, it is not bound despite of the common D-wave
mixing coming from Eq. (8) which lowers the ground state
energyTornqvist:1993ng ; Tornqvist:2004qy .

However, as discussed before, one should note that there is a strong
attraction at short distance in the $T_{cc}$ channel which will be modeled as
follows: $V_{\text{short}}(r)=V_{\omega}e^{-\mu_{\omega}^{2}r^{2}}$, where
$\mu_{\omega}^{2}=m_{\omega}^{2}-(m_{V}-m_{P})^{2}$ is taken to reproduce the
$r$ dependence in the lattice result Ikeda:2013vwa , which can be viewed as
the effect from omega meson exchange. Furthermore, the overall constant
$V_{\omega}=-$448 MeV is taken to reproduce the attraction extracted from the
constituent quark model discussed previously.

This short range potential in addition to the tensor potential (cf. Fig. 1)
leads to a loosely bound state of $T_{cc}$ in the molecular configuration with
binding energy $E_{b}$ and radius $\langle r\rangle$ given in Table 1. As
discussed before, short range attraction obtained in the quark model is
negligible in the $X(3872)$ channel, so that we do not add the additional
attraction in that channel. In fact, even if we add the relatively small
attraction there, the result does not change much. Their corresponding wave
functions are shown in Fig. 1. While the short range attraction is crucial in
providing the additional attraction, because the main attraction is coming
from the D-wave mixing, the wave function is extended to large distances due
to angular momentum barrier so that changing the form or attraction of the
potential below 0.3 fm do not give rise to any significant changes in the
molecular configuration.

When the formalism is applied to the deuteron case, in which case there is an
additional S-wave central potential, one finds $\gamma=25/3$, and the off
diagonal and the lower diagonal components are replaced by $\sqrt{8}$ and $-2$
in the second term of Eq. (8), respectively. Again, the strong attraction
coming from the D-wave mixing leads to the well known deuteron wave function
with a small binding and a large radius (cf. Table 1) Tornqvist:1993ng .

Therefore, we find that the molecular structures of the $X(3872)$ and $T_{cc}$
composed of heavy mesons are similar to that of the deuteron composed of the
nucleons.

##  IV Transverse momentum distribution of $d$ and 3He:

Since both the molecular tetraquarks and the deuteron are weakly bound broad
molecular states, their transverse momentum dependence in heavy ion collisions
should be determined by the distributions of their respective constituents at
the kinetic freeze-out point. This is so because for widely separated
molecules, their cross section with other hadrons should be just the sum of
that for the individual constituents. Such descriptions can be well encoded in
the coalescence model, where one could use the observed distribution of the
constituents.

On the other hand, the total yield could be determined earlier, where the
number of constituents are smaller than that at the kinetic freeze-out point
where the contributions of excited states are added through feed-down.

To check our idea and to determine the number of constituents at the formation
point, which involves subtracting the contributions of the feed-downs at the
formation time, we will try to fit the transverse momentum distributions as
well as the magnitudes of both the deuteron and 3He simultaneously.

For that purpose, we will use a two dimensional coalescence model with a
Wigner function of the Gaussian type scaled by $\sigma$ and the coalescence
volume $A$. We calculate the deuteron distribution in Pb–Pb collisions at
$\sqrt{s_{NN}}=2.76$ TeV ALICE:2015wav . Along the line, it is necessary to
know the transverse momentum distribution of constituents, the $\sigma$, $A$
and $R_{b}$. $\sigma$ is a parameter that can be determined by the size of the
produced hadron, and the $p_{T}$ distribution of constituents were fitted
using the $p_{T}$ distribution measured by ALICE collaboration and
parametrized as in the Appendix. $A$ is the kinetic freeze-out area. Finally,
$R_{b}$ is the ratio between the bare protons and the final protons at the
formation point. $(A,R_{b})$ can be determined by fitting the coalescence
model results to both the deuteron and 3He data. This is possible because in
the coalescence model, the deuteron production is proportional to
$R_{b}^{2}/A$ while the 3He production is proportional to $R_{b}^{3}/A^{2}$.

For the deuteron, we compare the coalescence model result with the $p_{T}$
distribution measured in Pb–Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV and
0–10% centrality by ALICE collaboration. For 3He, we compare the data at 0-20%
by scaling $A$ according to the charged particle multiplicity. The proton data
were taken at their respective centralities.

In the upper graph of Fig. 2, we show the transverse momentum distribution for
deuteron at 0–10% centrality using both the original two-dimensional
coalescence model (10) with different sizes and $\sigma=\infty$ model (14). In
this calculation we use $A=608~{}\text{fm}^{2}$, $R_{b}=0.368$, proton mass
938 MeV and neutron mass 939 MeV. Also, since $\sigma=\sqrt{8/3}r$ Ref.
Sun:2018mqq  and the radius of deuteron is taken to be $r_{d}=1.9$, we used
$\sigma=3.1$ fm. The 0–10% proton distribution was obtained using the average
of the 0–5% distribution and the 5–10% distribution. As can be seen in the
upper graph of Fig. 2, the coalescence formula with the physical radius
$r_{d}=1.9$ fm and that obtained in the $\sigma\rightarrow\infty$ limit lies
almost on top of each other and well reproduces the experimental results. The
effect of small radius becomes visible only when $\sigma<<1/\sqrt{mT}$ ($m$ is
the reduced mass of its constituent) such as $r=0.1$ fm shown in the figure.
In the lower graph of Fig. 2, the 3He data are well reproduced with the same
$R_{b}$, although at higher $p_{T}$ the model lies at the lower boundary of
the error bar.

Therefore, since the expected radius in the molecular configurations of both
the $X(3872)$ and $T_{cc}$ are larger than that of the deuteron, we can
calculate the transverse momentum dependence for these molecular tetraquark
configuration using $\sigma\rightarrow\infty$.

It is interesting to note that the factor $R_{b}$ obtained to simultaneously
fit the deuteron and 3He are close to what is expected of the fraction of bare
proton number using statistical hadronization model(SHM)Hwa:2004yg  at the
chemical freeze-out point. Therefore, we will also use the SHM with charm
quarks to estimate the corresponding $R_{b}$ factor for the $D$ and $D^{*}$
mesons that will be used as the input in the coalescence model for $X(3872)$
and $T_{cc}$.

Figure 2: The transverse momentum distributions of deuteron ALICE:2015wav ;
ALICE:2017nuf  (top figure) and 3He ALICE:2015wav  (bottom figure) in Pb-Pb
collisions at 2.76 TeV and the corresponding fits.

##  V Transverse momentum distribution of $X(3872)$ and $T_{cc}$:

Fig. 3 shows coalescence model result for the transverse momentum distribution
of $T_{cc}$ and $X(3872)$ at centrality 0–10% for Pb–Pb collisions at
$\sqrt{s_{NN}}=5.02$ TeV. Assuming a molecular configuration of $X(3872)$, its
production will be determined by the distribution of $D^{0}$ and
$\bar{D}^{*0}$ at the kinetic freeze-out point. We fitted the observed $D^{0}$
and $\bar{D}^{*+}$ data in Pb–Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV. As
the the $D^{*0}$ $p_{T}$ distribution is not currently measured, we used the
measured $\bar{D}^{*+}$ distribution here. Taking into account feed downs, we
expect that 31$\%$ of the measured $D^{0}$ will participate in the coalescence
(See Appendix for details). Yield of $D^{0}$ in Pb–Pb collisions at 5.02 TeV
is $dN_{D^{0}}/dy$ = 6.819 ALICE:2021rxa , so that the expected bare number is
2.12.

Since the chemical potential can be almost neglected, a molecular
configuration of $T_{cc}$ will have almost identical transverse momentum
distribution of $X(3872)$ as given in Fig. 3. The band in Fig. 3 shows the
statistical uncertainty of the input $D,D^{*}$ data transferred into the
coalescence result. One also notes that for the radii given in Table 1, the
results are almost identical to the result given in the
$\sigma\rightarrow\infty$ limit.

We also plot, the expected transverse momentum distribution of $T_{cc}$ as
given in Ref. Cho:2019syk , obtained assuming that it has a compact multiquark
configuration. The result in Ref. Cho:2019syk  given for collision energy at
2.76 TeV, has been multiplied by 1.63 to take into account the difference in
the charm quark number squared at 5.02 TeV as given in ExHIC:2017smd  to
compare with the result for the molecular configuration calculated here. As
can be seen in the figure, the expected transverse momentum distributions for
two possible configurations of the $T_{cc}$ are markedly different in addition
to the difference in the total yields. Therefore, once the transverse momentum
distribution of $T_{cc}$ is measure, we will be able to discriminate its
structure.

Figure 3: The transverse momentum distribution of $X(3872)$ and $T_{cc}$
(0–10% event). The bands for molecular configurations are due to the
uncertainties in the input data.

##  VI Total yields and summary

One of us has previously noted that for an exotic hadron to be compact, it has
to include multiple heavy quarks and that it could therefore be abundantly
produced in a relativistic heavy ion collision, where many heavy quark-
antiquark pairs are produced Lee:2007tn ; Proceedings:2007ctk . Later, two of
us with the ExHIC collaboration have shown that one could discriminate whether
an exotic is a molecular configuration or a compact multiquark state, from its
production rate in a heavy ion collision ExHIC:2017smd ; ExHIC:2010gcb ;
ExHIC:2011say . Recently, the CMS collaboration at CERN measured the $X(3872)$
for the first time in a heavy ion collision, and found its production ratio to
$\psi(2S)$ to be anomalously larger than that in pp collision CMS:2021znk .
According to the ExHIC papers, enhanced production of $X(3872)$ in heavy ion
collisions implies that its structure is a loosely bound molecular
configuration. Therefore, we estimate the total yields of hadrons per unit
rapidity at central rapidity in a collision event by integrating the
calculated transverse momentum distributions in the $\sigma\rightarrow\infty$
limit.

Table 2 shows the ratio between the total yield of $d$ and 3He calculated
using the coalescence model and the statistical model with $T_{H}=156$ MeV and
$\mu_{B}=0.7$ MeV. The coalescence model reproduces the SHM result for the
deuteron yield but slightly underestimates the corresponding value for 3He.

The total yields for the tetraquark states are given in Table 3. Assuming the
tetraquark states are of molecular configurations, the coalescence model
results, which are the same for the $X(3872)$ and the $T_{cc}$, are
appreciably larger than that from the statistical model result with charm
(SHMc) Andronic:2021erx : the ratio is 2.47. It should be noted that unlike
the deuteron case where one has to consider the feeddown contribution for both
constituents, one only eliminates those for $D$ and not for $D^{*}$ or
$\bar{D}^{*}$. However, assuming that the $T_{cc}$ is composed of a compact
multiquark configuration, the coalescence model gives a smaller value than
that expected from the SHMc. The last colume in Table 3 shows the yield ratio
of tetraquark states to $\psi(2S)$ in the SHMc. As can be seen in the table,
the ratio is $0.806\pm 0.234$ if the tetraquark is assumed to be of molecular
configuration. Although the experimental measurement is at high transverse
momentum and that the data takes into account the different branching ratios
into $J/\psi\pi\pi$, considering further suppression of $\psi(2S)$ in heavy
ion collision, the observed ratio seems consistent with a molecular picture
for $X(3872)$.

Nucleus | $N^{Nucleus}_{SHM}/N^{p}_{SHM}$ | $N^{Nucleus}_{coal}/N^{p}_{SHM}$
---|---|---
$d$ | $9.07\times 10^{-3}$ | $8.84\times 10^{-3}$
3He | $2.68\times 10^{-5}$ | $2.03\times 10^{-5}$
Table 2: The yield ratio of light nucleus with proton in Pb–Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV. For deuteron and 3He the centralities are 0–10 % and 0–20 %, respectively. Tetraquark | $dN_{coal}/dy$ | $N_{coal}/N^{X(3872)}_{SHMc}$ | $N_{coal}/N^{\psi(2S)}_{SHMc}$
---|---|---|---
$DD^{*}$ molecule | $(2.45\pm 0.71)\times 10^{-3}$ | $2.47\pm 0.716$ | $0.806\pm 0.234$
Compact 4$q$ | 6.2 $\times 10^{-4}$ | $6.25\times 10^{-1}$ | 0.204
Table 3: The first column shows the total yield of the tetraquark depending on
its structure calculated by the coalescence model in Pb–Pb collisions at
$\sqrt{s_{NN}}=5.02$ TeV at 0–10% centrality. The remaining columns show their
ratios to the statistical hadronization model with charm (SHMc)
Andronic:2021erx . Here we used $dN_{\psi(2S)}/dy=3.04\times 10^{-3}$ and
$N_{X(3872)}/N_{\psi(2S)}=0.326$ obtained in SHMc.

We have argued why $X(3872)$ should be a molecular configuration while
$T_{cc}$ could be either a molecular configuration or a compact mutlqiaruk
state. Our analysis suggest that the anomalously large yield ratio
$X(3872)/\psi(2S)$ in heavy ion collision indeed seems consistent with the
moelcular picture of $X(3872)$. However, as the present experimental result is
at higher transverse momentum than what we calculated here, further
measurements of these as well as the newly observed exotics Gershon:2022xnn
in heavy ion collision will be able to clearly discriminate a compact
multiquark configuration and a loosely bound molecular configuration.

##  Appendix A $\pi$-exchange potential

We use the following tensor potential as given in Ref. Tornqvist:1993ng ;
Tornqvist:2004qy .

$\displaystyle T_{\pi}(r)$ $\displaystyle=\frac{\mu^{2}r^{2}+3\mu
r+3}{m_{\pi}^{3}r^{3}}e^{-\mu r}-\frac{\Lambda^{2}r^{2}+3\Lambda
r+3}{m_{\pi}^{3}r^{3}}e^{-\Lambda r}$ (9)
$\displaystyle\quad-\frac{(\Lambda^{2}-\mu^{2})(\Lambda
r+1)}{2m_{\pi}^{3}r}e^{-\Lambda r},$

where $\Lambda=$1300 MeV and $\mu^{2}=m_{\pi}^{2}-(m_{V}-m_{P})^{2}$ stemming
from the recoil effect between pseudoscalar and vector mesons.

##  Appendix B Coalescence model

We use the following two particle coalescence model in two-transverse
dimension at central rapidity, which for the deuteron combines two nucleons.

$\displaystyle\frac{d^{2}N_{d}}{d^{2}P_{T}}=g_{d}$ $\displaystyle\int
d^{2}x_{1}d^{2}x_{2}d^{2}p_{1T}d^{2}p_{2T}\frac{d^{2}N_{p}}{dA_{L}d^{2}p_{1T}}\frac{d^{2}N_{n}}{dA_{L}d^{2}p_{2T}}$
$\displaystyle\times
W_{d}(\vec{r},\vec{k})\delta^{(2)}\left(\vec{P}_{T}-\vec{p}_{1T}-\vec{p}_{2T}\right),$
(10)

where $g_{d}$, $A_{L}$ and $W_{d}(\vec{r},\vec{k})$ are the deuteron
statistical factor divided by those of the constituents, the kinetic freeze-
out area in the Lab frame and the Winger function of deuteron in the CM frame,
respectively. Here, we use a two-dimensional Gaussian-type Wigner function
given as follows:

$\displaystyle W_{d}(\vec{r},\vec{k})=4\exp{\left(-\frac{r_{1}^{\prime
2}}{\sigma^{2}}-\sigma^{2}k_{1}^{\prime 2}\right)}.$ (11)

As emphasised in Cho:2019lxb , to properly take into account deuterons with
large transverse momenta, the Wigner function should be defined in the
deuteron rest frame with the following relative distance and momentum:

$\displaystyle\vec{R}^{\,\prime}=\frac{m_{1}\vec{x}^{\,\prime}_{1T}+m_{2}\vec{x}^{\,\prime}_{2T}}{m_{1}+m_{2}},\quad\vec{r}^{\,\prime}_{1}=\vec{x}^{\,\prime}_{1T}-\vec{x}^{\,\prime}_{2T},$
$\displaystyle\vec{k}^{\,\prime}=\vec{p}^{\,\prime}_{1T}+\vec{p}^{\,\prime}_{2T},\quad\vec{k}^{\,\prime}_{1}=\frac{m_{2}\vec{p}^{\,\prime}_{1T}-m_{1}\vec{p}^{\,\prime}_{2T}}{m_{1}+m_{2}}.$
(12)

The parameter $\sigma$ in Eq. (11) is related to the deuteron radius by
$\sigma=\sqrt{8/3}r_{d}$. Assuming a spatially homogeneous
$\frac{d^{2}N_{p}}{dA_{L}d^{2}p_{1T}}$, the transverse momentum distribution
reduces as follows:

$\displaystyle\frac{d^{2}N_{d}}{d^{2}P_{T}}=\frac{g_{d}}{A}$
$\displaystyle(2\sqrt{\pi})^{2}\sigma^{2}\int
d^{2}p_{1T}d^{2}p_{2T}\frac{d^{2}N_{p}}{d^{2}p_{1T}}\frac{d^{2}N_{n}}{d^{2}p_{2T}}$
$\displaystyle\times\exp{(-\sigma^{2}k_{1}^{\prime
2})}\delta^{(2)}\left(\vec{P}_{T}-\vec{p}_{1T}-\vec{p}_{2T}\right).$ (13)

In the $\sigma\rightarrow\infty$ limit,

$\displaystyle\frac{d^{2}N_{d}}{d^{2}P_{T}}$
$\displaystyle=g_{d}(2\pi)^{2}\left(\frac{\gamma}{A}\right)\frac{d^{2}N_{p}}{d^{2}p_{pT}}|_{\vec{p}_{pT}=\frac{m_{p}}{M}\vec{P}_{T}}\frac{d^{2}N_{n}}{d^{2}p_{nT}}|_{\vec{p}_{nT}=\frac{m_{n}}{M}\vec{P}_{T}},$
(14)

For the deuteron distribution in Pb–Pb collisions at in $\sqrt{s_{NN}}=2.76$
TeV with 0–10% centrality, we use $A=608~{}\text{fm}^{2}$. For 3He, we compare
our calculation to the data by ALICE collaboration measured for 0–20%
centrality ALICE:2015wav . 0–20% proton distribution was obtained by averaging
the 0–5%, 5–10% and 10–20% distributions. Because we have a different
centrality, the kinetic freeze-out area $A$ that we used in the deuteron fit
will be rescaled by the ratio of the charged particle multiplicities between
different centralities. This gives $A=507$ fm2.

To estimate the feed down contribution, we use the SHM Hwa:2004yg

$\displaystyle N^{stat}_{h}$
$\displaystyle=V_{H}\frac{g_{h}}{2\pi^{2}}\int^{\infty}_{0}\frac{p^{2}dp}{\gamma^{-1}_{h}e^{E_{h}/T_{H}}\pm
1},$
$\displaystyle\approx\frac{\gamma_{h}g_{h}V_{H}}{2\pi^{2}}m_{h}^{2}T_{H}K_{2}\left(\frac{m_{h}}{T_{H}}\right),$
(15)

where $g_{h}$ is the degeneracy factor of hadron, $\gamma_{h}$ is fugacity,
$K_{2}$ is the modified Bessel function of the second kind, and $V_{H}$ and
$T_{H}=156$ MeV are the hadronization volume and temperature at chemical
freeze-out.

For the calculation of molecular configuration of $X(3872)$ and $T_{cc}$, we
present the result for Pb–Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV with 0–10%
centrality. The coalescence area $A$ in Pb–Pb collisions at
$\sqrt{s_{NN}}=5.02$ TeV was rescaled according to the charged particle
multiplicity: using $dN/d\eta^{\text{Pb-Pb},2.76\
\text{TeV}}_{0-10\%}=1447.5\pm 39$ Alice:PbPb_2.76  and $dN/d\eta^{\text{Pb-
Pb},5.02\ \text{TeV}}_{0-10\%}=1764\pm 25$ Alice:PbPb_5.02 , $A$ was taken to
be $A_{0-10\%}^{5.02\ \text{TeV}}=$ 740.5 fm2. We fitted the observed $D^{0}$
and $\bar{D}^{*+}$ data in Pb–Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV. As
the $D^{*0}$ $p_{T}$ distribution is not currently measured, we used the
measured $\bar{D}^{*+}$ instead. We again use the statistical model to
estimate the fraction of the $D^{0}$ participating in coalescence. The decay
channels $D^{*}(2007)^{0}\rightarrow D^{0}\pi^{0}$ (BR = 64.7$\%$),
$D^{*}(2007)^{0}\rightarrow D^{0}\gamma$ (BR = 35.3$\%$), and
$D^{*}(2010)^{+}\rightarrow D^{0}\pi^{+}$ (BR = 67.7$\%$) are considered. We
then find that 31$\%$ of the measured $D^{0}$ participates in the coalescence.

##  Appendix C Parameterized transverse momentum distributions

Here, we provide the distributions of the constituents needed in the
coalescence calculations, which were obtained by fitting the transverse
momentum distribution measured in experiments. An exponential function and a
power-law type function were used for low and high $p_{T}$, respectively:

$\displaystyle\frac{dN}{d^{2}p_{T}dy}$
$\displaystyle=ae^{-b(p_{T}/p_{1})^{c}}\quad\quad\qquad(\ p_{T}<p_{c}),$
$\displaystyle=\frac{d}{[1+(p_{T}/p_{0})^{2}]^{e}}\quad(\ p_{T}>p_{c}),$ (16)

where $a,\ b,\ c,\ d,\ e$ and $p_{0}$ are fitting parameters. Here,
$p_{1}=1.0$ GeV and $|y|<0.5$. The values of $p_{c}$ and the fitting
parameters are shown in Table 4. Our form in Eq. (16) is close to the Blast
wave formAlice:PbPb_2.76  but better reproduces the proton data for $p_{T}>3$
GeV.

|  $\sqrt{s_{NN}}$ (GeV) | $\%$ | a (GeV)-2 | b | c | d (GeV)-2 | e |  $p_{0}$ (GeV) |  $p_{c}$ (GeV)
---|---|---|---|---|---|---|---|---|---
$p$ | 2.76 | 0–5 | 4.551 | 0.376 | 2.282 | 6.670 | 7.995 | 3.588 | 2.0
$p$ | 2.76 | 5–10 | 3.727 | 0.3477 | 2.42 | 7.384 | 5.55 | 2.644 | 2.0
$p$ | 2.76 | 10–20 | 2.937 | 0.3812 | 2.333 | 4.415 | 6.719 | 3.149 | 2.0
$D^{0}_{u}$ | 5.02 | 0–10 | 0.7537 | 0.3224 | 1.788 | 3.028$\times 10^{8}$ | 3.053 | 0.08034 | 4.5
$D^{0}_{l}$ | 5.02 | 0–10 | 0.5885 | 0.3463 | 1.836 | 9.959$\times 10^{12}$ | 3.013 | 0.01264 | 4.0
$D^{*+}_{u}$ | 5.02 | 0–10 | 0.2673 | 0.41 | 1.568 | 1.62$\times 10^{10}$ | 2.824 | 0.02545 | 4.5
$D^{*+}_{l}$ | 5.02 | 0–10 | 0.204 | 0.428 | 1.55 | 4.723$\times 10^{8}$ | 3.057 | 0.06464 | 4.5
Table 4:  The fitting parameters. The subscripbs in the first column
corresponds to the fits for the upper (u) and lower (l) bounds in Fig. 3,
respectively.

In Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV Alice:PbPb_2.76 , the
transverse momentum distribution of protons was used for deuteron and 3He
calculations.

In Pb-Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV ALICE:2021rxa , the transverse
momentum distribution of $D$ and $D^{*}$ meson were used for $X(3872)$
calculations.

## Acknowledgements

The work was supported by Samsung Science and Technology Foundation under
Project Number SSTF-BA1901-04, and by the National Research Foundation of
Korea (NRF) grant funded by the Korea government (MSIT) under Project
No.2018R1A5A1025563, No. 2019R1A2C1087107, No. 2021R1A2C1009486, and No.
2022R1A2C1011549.

## References

  * (1) https://home.cern/news/news/physics/59-new-hadrons-and-counting
  * (2) S. K. Choi et al. [Belle], Phys. Rev. Lett. 91, 262001 (2003) doi:10.1103/PhysRevLett.91.262001 [arXiv:hep-ex/0309032 [hep-ex]].
  * (3) R. Aaij et al. [LHCb], Phys. Rev. Lett. 115, 072001 (2015) doi:10.1103/PhysRevLett.115.072001 [arXiv:1507.03414 [hep-ex]].
  * (4) R. Aaij et al. [LHCb], Phys. Rev. Lett. 122, no.22, 222001 (2019) doi:10.1103/PhysRevLett.122.222001 [arXiv:1904.03947 [hep-ex]].
  * (5) V. M. Abazov et al. [D0], Phys. Rev. Lett. 117, no.2, 022003 (2016) doi:10.1103/PhysRevLett.117.022003 [arXiv:1602.07588 [hep-ex]].
  * (6) R. Aaij et al. [LHCb], Nature Phys. 18, no.7, 751-754 (2022) doi:10.1038/s41567-022-01614-y [arXiv:2109.01038 [hep-ex]].
  * (7) A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 128, no.3, 032001 (2022) doi:10.1103/PhysRevLett.128.032001 [arXiv:2102.13048 [hep-ex]].
  * (8) A. Park, S. H. Lee, T. Inoue and T. Hatsuda, Eur. Phys. J. A 56, no.3, 93 (2020) doi:10.1140/epja/s10050-020-00078-z [arXiv:1907.06351 [hep-ph]].
  * (9) S. Zouzou, B. Silvestre-Brac, C. Gignoux and J. M. Richard, Z. Phys. C 30, 457 (1986) doi:10.1007/BF01557611
  * (10) W. Park and S. H. Lee, Nucl. Phys. A 925, 161-184 (2014) doi:10.1016/j.nuclphysa.2014.02.008 [arXiv:1311.5330 [nucl-th]].
  * (11) W. Park, S. Noh and S. H. Lee, Nucl. Phys. A 983, 1-19 (2019) doi:10.1016/j.nuclphysa.2018.12.019 [arXiv:1809.05257 [nucl-th]].
  * (12) S. Noh, W. Park and S. H. Lee, Phys. Rev. D 103, 114009 (2021) doi:10.1103/PhysRevD.103.114009 [arXiv:2102.09614 [hep-ph]].
  * (13) Y. Ikeda, B. Charron, S. Aoki, T. Doi, T. Hatsuda, T. Inoue, N. Ishii, K. Murano, H. Nemura and K. Sasaki, Phys. Lett. B 729, 85-90 (2014) doi:10.1016/j.physletb.2014.01.002 [arXiv:1311.6214 [hep-lat]].
  * (14) N. A. Tornqvist, Z. Phys. C 61, 525-537 (1994) doi:10.1007/BF01413192 [arXiv:hep-ph/9310247 [hep-ph]].
  * (15) N. A. Tornqvist, Phys. Lett. B 590, 209-215 (2004) doi:10.1016/j.physletb.2004.03.077 [arXiv:hep-ph/0402237 [hep-ph]].
  * (16) J. Adam et al. [ALICE], Phys. Rev. C 93 (2016) no.2, 024917 doi:10.1103/PhysRevC.93.024917 [arXiv:1506.08951 [nucl-ex]].
  * (17) R. C. Hwa and X. N. Wang, World Scientific, 2004, ISBN 978-981-238-077-7, 978-981-4488-08-2 doi:10.1142/5029
  * (18) K. J. Sun, C. M. Ko and B. Dönigus, Phys. Lett. B 792, 132-137 (2019) doi:10.1016/j.physletb.2019.03.033 [arXiv:1812.05175 [nucl-th]].
  * (19) S. Acharya et al. [ALICE], Eur. Phys. J. C 77 (2017) no.10, 658 doi:10.1140/epjc/s10052-017-5222-x [arXiv:1707.07304 [nucl-ex]].
  * (20) S. Acharya et al. [ALICE], JHEP 01, 174 (2022) doi:10.1007/JHEP01(2022)174 [arXiv:2110.09420 [nucl-ex]].
  * (21) S. Cho and S. H. Lee, Phys. Rev. C 101, no.2, 024902 (2020) doi:10.1103/PhysRevC.101.024902 [arXiv:1907.12786 [nucl-th]].
  * (22) S. Cho et al. [ExHIC], Prog. Part. Nucl. Phys. 95, 279-322 (2017) doi:10.1016/j.ppnp.2017.02.002 [arXiv:1702.00486 [nucl-th]].
  * (23) S. H. Lee, S. Yasui, W. Liu and C. M. Ko, Eur. Phys. J. C 54, 259-265 (2008) doi:10.1140/epjc/s10052-007-0516-z [arXiv:0707.1747 [hep-ph]].
  * (24) S. H. Lee in “Proceedings, Workshop on Heavy Ion Collisions at the LHC: Last Call for Predictions: Geneva, Switzerland, May 14 - June 8, 2007,” J. Phys. G 35, 054001 (2008) doi:10.1088/0954-3899/35/5/054001 [arXiv:0711.0974 [hep-ph]].
  * (25) S. Cho et al. [ExHIC], Phys. Rev. Lett. 106, 212001 (2011) doi:10.1103/PhysRevLett.106.212001 [arXiv:1011.0852 [nucl-th]].
  * (26) S. Cho et al. [ExHIC], Phys. Rev. C 84, 064910 (2011) doi:10.1103/PhysRevC.84.064910 [arXiv:1107.1302 [nucl-th]].
  * (27) A. Andronic, P. Braun-Munzinger, M. K. Köhler, A. Mazeliauskas, K. Redlich, J. Stachel and V. Vislavicius, JHEP 07, 035 (2021) doi:10.1007/JHEP07(2021)035 [arXiv:2104.12754 [hep-ph]].
  * (28) T. Gershon [LHCb], [arXiv:2206.15233 [hep-ex]].
  * (29) S. Cho, K. J. Sun, C. M. Ko, S. H. Lee and Y. Oh, Phys. Rev. C 101, no.2, 024909 (2020) doi:10.1103/PhysRevC.101.024909 [arXiv:1905.09774 [nucl-th]].
  * (30) B. Abelev et al.  (ALICE Collaboration), Phys. Rev. C 88, 044910
  * (31) J. Adam et al.  (ALICE Collaboration), Phys. Rev. Lett. 166, 222302