autotabularcenter/.style= file=#1, after head=  , table head= , late after
line= , table foot= , late after last line= , command=,
autobooktabularcenter/.style= file=#1, after head=  , table head=  , late
after line= , table foot= , late after last line= , command=,

Research Note

# On the evaluation of (meta-)solver approaches

Roberto Amadini<EMAIL_ADDRESS>Maurizio Gabbrielli
<EMAIL_ADDRESS>Department of Computer Science and Engineering,
University of Bologna, Italy Tong Liu<EMAIL_ADDRESS>Meituan, Beijing, China
Jacopo Mauro<EMAIL_ADDRESS>Department of Mathematics and Computer
Science, University of Southern Denmark, Denmark

###### Abstract

Meta-solver approaches exploits a number of individual solvers to potentially
build a better solver. To assess the performance of meta-solvers, one can
simply adopt the metrics typically used for individual solvers (e.g., runtime
or solution quality), or employ more specific evaluation metrics (e.g., by
measuring how close the meta-solver gets to its virtual best performance). In
this paper, based on some recently published works, we provide an overview of
different performance metrics for evaluating (meta-)solvers, by underlying
their strengths and weaknesses.

##  1 Introduction

A famous quote attributed to Aristotle says that “the whole is greater than
the sum of its parts”. This principle has been applied in several contexts,
including the field of constraint solving and optimization. Combinatiorial
problems arising from application domains such as scheduling, manufacturing,
routing or logistics can be tackled by combining and leveraging the
complementary strengths of different solvers to create a better global _meta-
solver_.111Meta-solvers are sometimes referred in the literature as _portfolio
solvers_ , because they take advantage of a “portfolio” of different solvers.

Several approaches for combining solvers and hence creating effective meta-
solvers have been developed. Over the last decades we witnessed the creation
of new Algorithm Selection Kotthoff (2016) and Configuration Hoos (2012)
approaches222A fine-tuned solver can be seen as a meta-solver where we
consider different configurations of the same solver as different solvers.
that reached peak results in various solving competitions SAT competition
(2021); Stuckey . (2010); ICAPS (2021). To compare different meta-solvers, new
competitions were created, e.g., the 2015 ICON challenge Kotthoff . (2017) and
2017 OASC competition on algorithm selection Lindauer . (2019). However, the
discussion of why a particular evaluation metric has been chosen to rank the
solvers is lacking.

We believe that further study on this issue is necessary because often meta-
solvers are evaluated on heterogeneous scenarios, characterized by a different
number of problems, different timeouts, and different individual solvers from
which the meta-solvers approaches are built. In this paper, starting from some
surprising results presented by sunnyas2 showing dramatic ranking changes with
different, but reasonable, metrics we would like to draw more attention on the
evaluation of meta-solvers approaches by shedding some light on the strengths
and weaknesses of different metrics.

##  2 Evaluation metrics

Before talking about the evaluation metrics, we should spend some words on
what we need to evaluate: the solvers. In our context, a solver is a program
that takes as input the description of a computational problem in a given
language, and returns an observable outcome providing zero or more solutions
for the given problem. For example, for decision problems the outcome may be
simply “yes” or “no” while for optimization problems we might be interested in
the sub-optimal solutions found along the search. An _evaluation metric_ , or
performance metric, is a function mapping the outcome of a solver on a given
instance to a number representing “how good” the solver on this instance is.

An evaluation metric is often not just defined by the output of the
(meta-)solver but can also be influenced by other actors such as the
computational resources available, the problems on which we evaluate the
solver, and the other solvers involved in the evaluation. For example, it is
often unavoidable to set a _timeout_ $\tau$ on the solver’s execution when
there is no guarantee of termination in a reasonable amount of time (e.g., NP-
hard problems). Timeouts makes the evaluation feasible, but inevitably couple
the evaluation metric to the execution context. For this reason, the
evaluation of a meta-solver should also take into account the _scenario_ that
encompasses the solvers to evaluate, the instances used for the validation,
and the timeout. Formally, at least for the purposes of this paper, we can
define a scenario as a triple $(\mathcal{I},\mathcal{S},\tau)$ where:
$\mathcal{I}$ is a set of problem instances, $\mathcal{S}$ is a set of
individual solvers, $\tau\in(0,+\infty)$ is a timeout such that the outcome of
solver $s\in S$ over instance $i\in\mathcal{I}$ is always measured in the time
interval $[0,\tau)$.

Evaluating meta-solvers over heterogeneous scenarios is complicated by the
fact that the set of instances, solvers and the timeout can have high
variability. As we shall see in Sect. 2.3, things are even trickier in
scenarios including optimization problems.

###  2.1 Absolute vs relative metrics

A sharp distinction between evaluation metrics can be drawn depending on
whether their value depend on the outcome of other solvers or not. We say that
an evaluation metric is _relative_ in the former case, _absolute_ otherwise.
For example, a well-known absolute metric is the _penalized average runtime_
with penalty $\lambda\geq 1$ (PARλ) that compares the solvers by using the
average solving runtime and penalizes the timeouts with $\lambda$ times the
timeout.

Formally, let $\mathtt{time}(i,s,\tau)$ be the function returning the runtime
of solver $s$ on instance $i$ with timeout $\tau$, assuming
$\mathtt{time}(i,s,\tau)=\tau$ if $s$ cannot solve $i$ before the timeout
$\tau$. For optimization problems, we consider the runtime as the time taken
by $s$ to solve $i$ to optimality333If $s$ cannot solve $i$ to optimality
before $\tau$, then $\mathtt{time}(i,s)=\tau$ even if sub-optimal solutions
are found. assuming w.l.o.g. that an optimization problem is always a
minimization problem. We can define PARλ as follows.

######  Definition 1 (Penalized Average Runtime).

Let $(\mathcal{I},\mathcal{S},\tau)$ be a scenario, the PARλ score of solver
$s\in\mathcal{S}$ over $\mathcal{I}$ is given by
$\dfrac{1}{|\mathcal{I}|}\sum\limits_{i\in\mathcal{I}}\mathtt{par}_{\lambda}(i,s,\tau)$
where:

$\mathtt{par}_{\lambda}(i,s,\tau)=\begin{cases}\mathtt{time}(i,s,\tau)&\text{if
$\mathtt{time}(i,s,\tau)<\tau$}\\\
\lambda\cdot\tau&\text{otherwise.}\end{cases}$

Well-known PAR measures are, e.g., the $\mathsf{PAR}_{2}$ adopted in the SAT
competitions SAT competition (2021) or the $\mathsf{PAR}_{10}$ used by
lindauer2019algorithm. Evaluating (meta-)solvers with scenarios having
different timeouts should imply a _normalization_ of PARλ in a fixed range to
avoid misleading comparisons.

Another absolute metric for decision problems is the number (or percentage) of
instances solved where ties are broken by favoring the solver minimizing the
average running time, i.e., minimizing the $\mathsf{PAR}_{1}$ score. This
metric has been used in various tracks of the planning competition ICAPS
(2021), the XCSP competition XCSP Competition (2019), and the QBF evaluations
QBFEVAL (2021).

A well-established relative metric is instead the _Borda count_ , adopted for
example by the MiniZinc Challenge Stuckey . (2014) for both single solvers and
meta-solvers. The Borda count is a family of voting rules that can be applied
to the evaluation of a solver by considering the comparison as an election
where the solvers are the candidates and the problem instances are the voters.
The MiniZinc challenge uses a variant of Borda444In the original definition,
the lowest-ranked candidate gets 0 points, the next-lowest 1 point, and so on.
where each solver scores points proportionally to the number of solvers it
beats. Assuming that $\mathtt{obj}(i,s,t)$ is the best objective value found
by solver $s$ on optimization problem $i$ at time $t$, with
$\mathtt{obj}(i,s,t)=\infty$ when no solution is found at time $t$, the
MiniZinc challenge score is defined as follows.

######  Definition 2 (MiniZinc challenge score).

Let $(\mathcal{S},\mathcal{I},\tau)$ be a scenario where
$\mathcal{I}=\mathcal{I}_{dec}\cup\mathcal{I}_{opt}$ with $\mathcal{I}_{dec}$
decision problems and $\mathcal{I}_{opt}$ optimization problems. The MiniZinc
challenge (MZNC) score of $s\in\mathcal{S}$ over $\mathcal{I}$ is
$\sum\limits_{i\in\mathcal{I},s^{\prime}\in
S\setminus\\{s\\}}m_{s}(i,s^{\prime},\tau)$ where:

$m_{s}(i,s^{\prime},\tau)=\begin{cases}0&\text{if
$\mathtt{unknown}(i,s)\vee\mathtt{better}(i,s^{\prime},s)$}\\\ 1&\text{if
$\mathtt{better}(i,s,s^{\prime})$}\\\ 0.5&\text{if
$\mathtt{time}(i,s,\tau)=\mathtt{time}(i,s^{\prime},\tau)$}\\\ ~{}&\text{and
$\mathtt{obj}(i,s,\tau)=\mathtt{obj}(i,s^{\prime},\tau)$}\\\
\dfrac{\mathtt{time}(i,s^{\prime},\tau)}{\mathtt{time}(i,s,\tau)+\mathtt{time}(i,s^{\prime},\tau)}&\text{otherwise}\end{cases}$

where the predicate $\mathtt{unknown}(i,s)$ holds if $s$ does not produce a
solution within the timeout: unknown(i,s) = (i ∈I_dec ∧time(i,s,τ)=τ) ∨(i
∈I_opt ∧obj(i,s,τ)=∞) and $\mathtt{better}(i,s,s^{\prime})$ holds if $s$
finishes earlier than $s^{\prime}$ or it produces a better solution:
better(i,s,s’) = time(i,s,τ) < time(i,s’,τ) = τ∨obj(i,s,τ) < obj(i,s’,τ)

This is clearly a relative metric because changing the set of available
solvers can affect the MiniZinc scores.

To handle the disparate nature of the scenarios when comparing meta-solvers
approaches, the evaluation function adopted in the ICON and OASC challenges
was relative: the _closed gap score_. This metric assigns to a meta-solver a
value in $[0,1]$ proportional to how much it closes the gap between the best
individual solver available, or single best solver ($\mathit{SBS}$), and the
virtual best solver ($\mathit{VBS}$), i.e., an oracle-like meta-solver always
selecting the best individual solver. The closed gap is actually a “meta-
metric”, defined in terms of another evaluation metric. Formally, if
$(\mathcal{I},\mathcal{S},\tau)$ is a scenario and $m$ an evaluation metric to
minimize, we have $m(i,\mathit{VBS},\tau)=\min\\{m(i,s,\tau)\mid
s\in\mathcal{S}\\}$ for each $i\in\mathcal{I}$ and
$\mathit{SBS}=\underset{s\in\mathcal{S}}{\operatorname{argmin}}\sum_{i\in\mathcal{I}}\,m(i,s,\tau)$.
We can define the closed gap as follows.

######  Definition 3 (Closed gap).

Let $m$ be an evaluation metric to minimize for a scenario
$(\mathcal{I},\mathcal{S},\tau)$, and let
$m_{\mathit{VBS}}=\sum_{i\in\mathcal{I}}m(i,\mathit{VBS},\tau)$ and
$m_{\mathit{SBS}}=\sum_{i\in\mathcal{I}}m(i,\mathit{SBS},\tau)$. The _closed
gap_ of a (meta-)solver $s$ w.r.t. $m$ on that scenario is:

$\dfrac{m_{\mathit{SBS}}-\sum_{i\in\mathcal{I}}m(i,s,\tau)}{m_{\mathit{SBS}}-m_{\mathit{VBS}}}$

If not specified, we will assume the closed gap computed w.r.t. the
$\mathsf{\mathsf{PAR}_{10}}$ score as done in the AS challenges 2015 and
2017.555In the 2015 edition, the closed gap was computed as
$1-\frac{m_{\mathit{SBS}}-m(i,s,\tau)}{m_{\mathit{SBS}}-m_{\mathit{VBS}}}=\frac{m_{\mathit{VBS}}-m(i,s,\tau)}{m_{\mathit{VBS}}-m_{\mathit{SBS}}}$.

Table 1: Comparison ASAP vs RF. The MZNC column reports the average MZNC score per scenario. Negative scores are in bold font. | Closed gap | MZNC | Better than other
---|---|---|---
Scenario | ASAP | RF | ASAP | RF | ASAP | RF
ASP-POTASSCO | 0.7444 | 0.5314 | 2.2235 | 2.6163 | 275 | 671
BNSL-2016 | 0.8463 | 0.7451 | 1.2830 | 3.0250 | 98 | 993
CPMP-2015 | 0.6323 | 0.1732 | 2.0501 | 2.3660 | 137 | 334
CSP-MiniZinc-Time-2016 | 0.6251 | 0.2723 | 2.1552 | 2.7214 | 17 | 53
GLUHACK-2018 | 0.4663 | 0.4057 | 1.9040 | 2.4528 | 62 | 147
GRAPHS-2015 | 0.758 | -0.6412 | 2.3045 | 3.3731 | 489 | 3663
MAXSAT-PMS-2016 | 0.5734 | 0.3263 | 1.4747 | 2.8616 | 66 | 439
MAXSAT-WPMS-2016 | 0.7736 | -1.1826 | 1.5168 | 2.4043 | 126 | 386
MAXSAT19-UCMS | 0.6583 | -0.2413 | 2.0893 | 2.5189 | 145 | 269
MIP-2016 | 0.35 | -0.3626 | 2.4035 | 2.4239 | 81 | 105
QBF-2016 | 0.7568 | -0.1366 | 1.8642 | 2.7154 | 193 | 467
SAT03-16_INDU | 0.3997 | 0.1503 | 2.1508 | 2.5812 | 491 | 1116
SAT12-ALL | 0.7617 | 0.6528 | 1.6785 | 2.8250 | 262 | 1227
SAT18-EXP | 0.5576 | 0.3202 | 1.9239 | 2.4998 | 61 | 164
TSP-LION2015 | 0.4042 | -19.1569 | 2.4352 | 2.6979 | 1115 | 1949
Tot. | 9.3077 | -18.1439 | 29.4573 | 40.0826 | 3618 | 11983
Tot.$-$TSP-LION2015 | 8.9035 | 1.013 | 27.0221 | 37.3846 | 2503 | 10034

###  2.2 A surprising outcome

An interesting outcome reported in sunnyas2 was the profound difference
between the closed gap and the MiniZinc challenge scores. Liu . compared the
performance of six meta-solvers approaches across 15 decision-problems
scenarios taken from ASlib Bischl . (2016) and coming from heterogeneous
domains such as Answer-Set Programming, Constraint Programming, Quantified
Boolean Formula, Boolean Satisfiability.

Tab. 1 reports the performance of ASAP and RF, respectively the best approach
according to the closed gap score and the MZNC score. The scores in the four
leftmost columns clearly show a remarkable difference rank if we swap the
evaluation metric. With the closed gap, ASAP is the best approach and RF the
worst among all the meta-solvers considered, while with the MZNC score RF
climbs to the first position while ASAP drops to the last position.

Another thing that catches the eye in Tab. 1 is the presence of _negative
scores_. This happens because, by definition, the closed gap has upper bound 1
(no meta-solver can improve the $\mathit{VBS}$) but not a fixed lower bound.
Hence, when the performance of the meta-solver is worse than the performance
of the single best solver, the closed gap drops below zero. While on a first
glance this seems reasonable—meta-solvers should perform no worse than the
individual solvers—it is worth noting that the penalty for performing worse
than the SBS also depends on the denominator
$m_{\mathit{SBS}}-m_{\mathit{VBS}}$. This means that in scenarios where the
performance of the $\mathit{SBS}$ is close to the perfect performance of the
$\mathit{VBS}$ this penalty can be significantly magnified. The TSP-LION2015
scenario is a clear example: the RF approach gets a penalization of more than
19 points, meaning that RF should perform flawlessly in about 20 other
scenarios to expiate this punishment. In fact, in TSP-LION2015 the
$\mathsf{\mathsf{PAR}_{10}}$ distributions of $\mathit{SBS}$ and
$\mathit{VBS}$ are very close: the $\mathit{SBS}$ is able to solve 99.65% of
the instances solved by the $\mathit{VBS}$, leaving little room for
improvement. RF scores -19.1569 while still solving more than 90% of the
instances of the scenario and having a difference with ASAP of slightly more
than 5% instances solved.

Why are the closed gap and the MZNC rankings so different? Looking at the
rightmost two columns in Tab. 1 showing, for each scenario, the number of
instances where an approach is faster than the other, one may conclude that RF
is far better than ASAP. In all the scenarios the number of instances where
its runtime is lower than ASAP runtime is greater than the number of instances
where ASAP is faster. Overall, it is quite impressive to see that RF beats
ASAP on 11983 instances while ASAP beats RF on 3618 times only.

An initial clue of why this happens is revealed in sunnyas2, where a
parametric version of MZNC score is used. In practice, Def. 2 is generalized
by assuming the performance of two solvers equivalent if their runtime
difference is below a given time threshold $\delta$.666Formally, if
$|\mathtt{time}(i,s,\tau)-\mathtt{time}(i,s^{\prime},\tau)|\leq\delta$ then
both $s$ and $s^{\prime}$ scores 0.5 points—note that if $\delta=0$ we get the
original MZNC score as in Def. 2. This variant was considered because a time
difference of few seconds could be considered irrelevant if solving a problem
can take minutes or hours.

Figure 1: Cumulative Borda count by varying the $\delta$ threshold.

The parametric MZNC score is depicted in Fig. 1, where different thresholds
$\delta$ are considered on the x-axis. It is easy to see how the performance
of ASAP and RF reverses when $\delta$ increases: ASAP bubbles from the bottom
to the top, while RF gradually sinks to the bottom.

Let us further investigate this anomaly.

Figure 2: Solve instances difference ASAP vs RF.

Fig. 2 shows the runtime distributions of the instances solved by ASAP and RF,
sorted by ascending runtime. We can see that ASAP solves more instances, but
for around 15k instances RF is never slower than ASAP. Summarizing, ASAP
solves more instances but RF is in general quicker when it solves an (often
easy) instance. This entails the significant difference between closed gap and
Borda metrics.

In our opinion, on the one hand, it is fair to think that ASAP performs better
than RF on these scenarios. The MZNC score seems to over-penalize ASAP w.r.t.
RF. Moreover, from Fig. 1 we can also note that for $\delta\leq 103$ the
parametric MZNC score of RF is still better, but 103 seconds looks quite a
high threshold to consider two performances as equivalent. On the other hand,
the closed gap score can also be over-penalizing due to negative outliers.

We also would like to point out that the definitions of $\mathit{SBS}$ found
in the literature do not clarify how it is computed in scenarios where the set
of instances $\mathcal{I}$ is split into test and training sets. Should the
$\mathit{SBS}$ be computed on the instances of the training set, the test set,
or the whole dataset $\mathcal{I}$? One would be inclined to use the test set
to select the $\mathit{SBS}$, but this choice might be problematic because the
test set is usually quite small w.r.t. the training set when using, e.g.,
cross-validation methods. In this case issues with negative outliers might be
amplified. If not clarified, this could lead to confusion. For example, in the
2015 ICON challenge the $\mathit{SBS}$ was computed by considering the entire
dataset (training and testing instances together). In the 2017 OASC instead
the $\mathit{SBS}$ was originally computed on the test set of the scenarios,
but then the results were amended by computing the $\mathit{SBS}$ on the
training set.

An alternative to the above metrics is the _speedup_ of a single solver w.r.t.
the $\mathit{SBS}$ or the $\mathit{VBS}$, i.e., how much a meta-solver can
improve a baseline solver. Tab. 2 reports, using the data of Liu . (2021), for
each meta-solver $s$ in a scenario $(\mathcal{I},\mathcal{S},\tau)$ the
average speedup computed as
$\dfrac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}\frac{\mathtt{time}(i,\mathit{VBS},\tau)}{\mathtt{time}(i,s,\tau)}$.
Unlike the closed gap, that has no lower bound, the speedup always falls in
$[0,1]$ with bigger values meaning better performance. We compared this with
the average normalized runtime score, computed as
$1-\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}\frac{\mathtt{time}(i,s,\tau)}{\tau}$
and the average closed gap score w.r.t $\mathsf{PAR}_{1}$. We use
$\mathsf{PAR}_{1}$ instead of $\mathsf{PAR}_{10}$ to be consistent with
speedup and normalized runtime, which do not get any penalization.

Table 2: Average closed gap, speedup, and normalized runtime. Peak performance in bold font. Meta-solver | Closed gap | Speedup | Norm. runtime
---|---|---|---
ASAP | 0.4866 | 0.4026 | 0.8829
sunny-as2 | 0.4717 | 0.4122 | 0.8879
autofolio | 0.4713 | 0.4110 | 0.8855
SUNNY-original | 0.4412 | 0.3905 | 0.8790
Zilla | 0.3416 | 0.3742 | 0.8753
Random Forest | -0.1921 | 0.3038 | 0.8507

The rank with speedup and normalised runtime is the same. The podium changes
if we use the closed gap: in this case sunny-as2 and autofolio lose one
position while ASAP rises from third to first position. However, as we shall
see in the next section, the generalization of these metrics to optimization
problems is not trivial.

###  2.3 Optimization problems

So far we have mainly talked about evaluating meta-solvers on decision
problems. While the MZNC score takes into account also optimization problems,
for the closed gap, (normalized) runtime, and speedup the generalization is
not as obvious as it might seem. Here using the runtime might not be the right
choice: often a solver cannot prove the optimality of a solution, even when it
actually finds it. Hence, the obvious alternative is to consider just the
objective value of a solution. But this value needs to be _normalized_ , and
to do so what bounds should we choose? Furthermore: how to reward a solver
that actually proves the optimality of a solution? And how to penalize solvers
that cannot find any solution?

If solving to optimality is not rewarded, metrics such as the ratio score of
the satisfiability track of the planning competition can be used.777This track
includes optimization problems where the goal is to minimize the length of a
plan. This score is computed as the ratio between the best known solution and
the best objective value found by the solver, giving 0 points in case no
solution is found.

A different metric that focuses on quickly reaching good solution is the
_area_ score, introduced in the MZNC starting from 2017. This metric computes
the integral of a step function of the solution value over the runtime
horizon. Intuitively, a solver that finds good solutions earlier can
outperform a solver that finds better solutions much later in the solving
stage.

Other attempts have been proposed to take into account the objective value and
the running times. For example, the ASP competition Calimeri . (2016) adopted
an elaborate scoring system that combines together the percentage of instances
solved within the timeout, the evaluation time, and the quality of a solution.
Similarly, paper_amai proposed a relative metric where each solver $s$ gets a
reward in $\\{0\\}\cup[\alpha,\beta]\cup\\{1\\}$ according to the objective
value $\mathtt{obj}(s,i,\tau)$ of the best solution it finds, with
$0\leq\alpha\leq\beta\leq 1$. If no solution is found then $s$ scores 0, if it
solves $i$ to optimality it scores 1, otherwise the score is computed by
linearly scaling $\mathtt{obj}(s,i,\tau)$ in $[\alpha,\beta]$ according to the
best and worst objective values find by any other available solver on problem
$i$.

###  2.4 Randomness and aggregation

We conclude the section with some remarks about randomness and data
aggregation.

When evaluating a meta-solver $s$ on scenario
$(\mathcal{I},\mathcal{S},\tau)$, it is common practice to partition
$\mathcal{I}$ into a training set $\mathcal{I}_{tr}$, on which $s$ “learns”
how to leverage its individual solvers, and a test set $\mathcal{I}_{ts}$
where the performance of $s$ on unforeseen problems is measured. In
particular, to prevent overfitting, it is possible to use a $k$-fold cross
validation by first splitting $\mathcal{I}$ into $k$ disjoint folds, and then
using in turn one fold as test set and the union of the other folds as
training set. In the AS challenge 2015 Lindauer . (2019) the submissions were
indeed evaluated with a 10-fold cross validation, while in the OASC in 2017
the dataset of the scenarios was divided only in one test set and one training
set. As also underlined by the organizers of the competition, this is risky
because it may reward a lucky meta-solver performing well on that split but
poorly on other splits.

Note that so far we have always assumed deterministic solvers, i.e., solvers
always providing the same outcome if executed on the same instance in the same
execution environment. Unfortunately, the scenario may contain randomized
solvers potentially producing different results with a high variability. In
this case, solvers should be evaluated over a number of runs and particular
care must be taken because the assumption that a solver can never outperform
the VBS would be no longer true.

A cautious choice to decrease the variance of model predictions would be to
repeat the $k$-fold cross validation $n>1$ times with different random splits.
However, this might imply a tremendous computational effort—the training phase
of a meta-solver might take hours or days—and therefore a significant energy
consumption. This issue is becoming an increasing concern. For example, in
their recent work DBLP:conf/cp/MatriconAFSH21 propose an approach to early
stop running an individual solver that it is likely to perform worse than
another solver on a subset of the instances of the scenario. In this way, less
resources are wasted for solvers that most likely will not bring any
improvement.

Finally, we spend a few words on the aggregation of the results. It is quite
common to use the arithmetic mean, or just the sum, when it comes to aggregate
the outcomes of a meta-solver over different problems of the same scenario
(e.g., when evaluating the results on the $n\cdot k$ test sets of a $k$-fold
cross validation repeated $n$ times). The same applies when evaluating
different scenarios. The choice of how to aggregate the metric values into a
unique value should however be motivated since the arithmetic mean can lead to
misleading conclusions when summarizing normalized benchmark Fleming Wallace
(1986). For example, to amortize the effect of outliers, one may use the
median or use the geometric mean to average over normalized numbers.

##  3 Conclusions

As it happens in many other fields, the choice of reasonable metrics can have
divergent effects on the assessment of (meta-)solvers. While these issues are
mitigated when comparing individual solvers in competitions having uniform
scenarios in term of size, difficulty, and nature, the comparison of meta-
solver approaches poses new challenges due to diversity of the scenarios on
which they are evaluated. Although it is impossible to define a fits-all
metric, we believe that we should aim at more robust metrics avoiding as much
as possible the under- and over-penalization of meta-solvers.

Particular care should be taken when using _relative_ measurements, because
the risk is to amplify small performance variations into large differences of
the metric’s value. Presenting the results in terms of orthogonal evaluation
metrics allows a better understanding of the (meta-)solvers performance, and
these insights may help the researchers to build meta-solvers that better fit
their needs, as well as to prefer an evaluation metric over another. Moreover,
well-established metrics may be combined into hybrid “meta-metrics” folding
together different performance aspects and handling the possible presence of
randomness.

## References

  * Amadini . (2016) paper_amaiAmadini, R., Gabbrielli, M. Mauro, J.  2016\.  Portfolio approaches for constraint optimization problems Portfolio approaches for constraint optimization problems.  Annals of Mathematics and Artificial Intelligence761-2229–246.
  * Bischl . (2016) bischl2016aslibBischl, B., Kerschke, P., Kotthoff, L., Lindauer, M., Malitsky, Y., Fréchette, A.Vanschoren, J.  2016\.  Aslib: A benchmark library for algorithm selection Aslib: A benchmark library for algorithm selection.  Artificial Intelligence23741–58.
  * Calimeri . (2016) ASPcompCalimeri, F., Gebser, M., Maratea, M. Ricca, F.  2016\.  Design and results of the Fifth Answer Set Programming Competition Design and results of the fifth answer set programming competition.  Artif. Intell.231151–181.  https://doi.org/10.1016/j.artint.2015.09.008 10.1016/j.artint.2015.09.008
  * Fleming Wallace (1986) DBLP:journals/cacm/FlemingW86Fleming, PJ. Wallace, JJ.  1986\.  How Not To Lie With Statistics: The Correct Way To Summarize Benchmark Results How not to lie with statistics: The correct way to summarize benchmark results.  Commun. ACM293218–221.  https://doi.org/10.1145/5666.5673 10.1145/5666.5673
  * Hoos (2012) DBLP:books/daglib/p/Hoos12Hoos, HH.  2012\.  Automated Algorithm Configuration and Parameter Tuning Automated algorithm configuration and parameter tuning.  Y. Hamadi, É. Monfroy F. Saubion (), Autonomous Search Autonomous search ( 37–71).  Springer.
  * ICAPS (2021) plancompICAPS.  2021\.  The International Planning Competition Web Page. The international planning competition web page.  https://www.icaps-conference.org/competitions/.  Accessed: 2021-12-10
  * Kotthoff (2016) kotthoff2016algorithmKotthoff, L.  2016\.  Algorithm selection for combinatorial search problems: A survey Algorithm selection for combinatorial search problems: A survey.  Data Mining and Constraint Programming Data mining and constraint programming ( 149–190).  Springer.
  * Kotthoff . (2017) kotthoff2015iconKotthoff, L., Hurley, B. O’Sullivan, B.  2017\.  The ICON Challenge on Algorithm Selection The ICON challenge on algorithm selection.  AI Magazine38291–93.
  * Lindauer . (2019) lindauer2019algorithmLindauer, M., van Rijn, JN. Kotthoff, L.  2019\.  The algorithm selection competitions 2015 and 2017 The algorithm selection competitions 2015 and 2017.  Artificial Intelligence27286–100.
  * Liu . (2021) sunnyas2Liu, T., Amadini, R., Mauro, J. Gabbrielli, M.  2021\.  sunny-as2: Enhancing SUNNY for Algorithm Selection sunny-as2: Enhancing SUNNY for algorithm selection.  J. Artif. Intell. Res.72329–376.
  * Matricon . (2021) DBLP:conf/cp/MatriconAFSH21Matricon, T., Anastacio, M., Fijalkow, N., Simon, L. Hoos, HH.  2021\.  Statistical Comparison of Algorithm Performance Through Instance Selection Statistical comparison of algorithm performance through instance selection.  LD. Michel (), 27th International Conference on Principles and Practice of Constraint Programming, CP 2021, Montpellier, France (Virtual Conference), October 25-29, 2021 27th international conference on principles and practice of constraint programming, CP 2021, montpellier, france (virtual conference), october 25-29, 2021 ( 210, 43:1–43:21).  Schloss Dagstuhl - Leibniz-Zentrum für Informatik.
  * QBFEVAL (2021) QBFEVALQBFEVAL.  2021\.  QBF Evaluations Web Page. Qbf evaluations web page.  http://www.qbflib.org/index_eval.php.  Accessed: 2021-12-10
  * SAT competition (2021) satcompSAT competition.  2021\.  The International SAT Competition Web Page. The international sat competition web page.  http://www.satcompetition.org/.  Accessed: 2021-12-10
  * Stuckey . (2010) DBLP:journals/constraints/StuckeyBF10Stuckey, PJ., Becket, R. Fischer, J.  2010\.  Philosophy of the MiniZinc challenge Philosophy of the minizinc challenge.  Constraints An Int. J.153307–316.  https://doi.org/10.1007/s10601-010-9093-0 10.1007/s10601-010-9093-0
  * Stuckey . (2014) mzn-challengeStuckey, PJ., Feydy, T., Schutt, A., Tack, G. Fischer, J.  2014\.  The MiniZinc Challenge 2008-2013 The MiniZinc Challenge 2008-2013.  AI Magazine35255–60.  http://www.aaai.org/ojs/index.php/aimagazine/article/view/2539
  * XCSP Competition (2019) XCSPcompXCSP Competition.  2019\.  XCSP Competition. Xcsp competition.  http://www.cril.univ-artois.fr/XCSP19/.  Accessed: 2021-12-10


# Probing the internal structures of $p\Omega$ and $\Omega\Omega$ with their
production at the LHC

Jie Pu<EMAIL_ADDRESS>Kai-Jia Sun  Corresponding author<EMAIL_ADDRESS>Chun-Wang Ma  Lie-Wen Chen  Corresponding author<EMAIL_ADDRESS>College
of Physics, Henan Normal University, Xinxiang 453007, China  Shanghai Research
Center for Theoretical Nuclear Physics, NSFC and Fudan University, Shanghai,
200438, China  Key Laboratory of Nuclear Physics and Ion-beam Application
(MOE), Institute of Modern Physics, Fudan University, Shanghai, 200433, China
Institute of Nuclear Science and Technology, Henan Academy of Science,
Zhengzhou, 450046, China  School of Physics and Astronomy, Shanghai Key
Laboratory for Particle Physics and Cosmology, and Key Laboratory for Particle
Astrophysics and Cosmology (MOE), Shanghai Jiao Tong University, Shanghai
200240, China

###### Abstract

The strange dibaryons $p\Omega$ (${}^{5}\rm{S}_{2}$) and $\Omega\Omega$
(${}^{1}\rm{S}_{0}$) are likely bound, existing either in molecular states
like the deuteron or as more exotic compact six-quark states. Here, we
investigate the production of these two dibaryons in Pb+Pb collisions at
$\sqrt{s_{NN}}$=2.76 TeV at the CERN Large Hadron Collider (LHC) within a
covariant coalescence model, which employs a blast-wave-like parametrization
for the phase-space configurations of constituent particles at freeze-out. For
the molecular states, the $p\Omega$ and $\Omega\Omega$ are produced via
$p$-$\Omega$ and $\Omega$-$\Omega$ coalescence, respectively, while for the
six-quark states, they are formed through $uudsss$ and $ssssss$ coalescence.
We find that the yield ratio $N_{p\Omega}/N_{\Omega}$ and
$N_{\Omega\Omega}/N_{\Omega}$ have a distinct centrality dependence between
the molecular and multi-quark states, thus offering a promising way for
distinguishing the two states. Our results suggest that the measurements of
$p\Omega$ and $\Omega\Omega$ production in relativistic heavy-ion collisions
can shed light on their internal structures.

††journal: Physics Letters B

##  1 Introduction

Exploring dibaryons is a long-standing theoretical and experimental challenge
in hadron and nuclear physics [1, 2]. A dibaryon is defined in quantum
chromodynamics (QCD) as a six-quark system with baryon number $B=2$ [3, 4],
and the only stable dibaryon observed so far is the deuteron [5]. The $H$
dibaryon, a hypothetical bound state with strangeness $S=-2$ and comprising
$uuddss$ quarks with $J^{\pi}=0^{+}$ and $I=0$, was first predicted by Jaffe
[6] within a bag-model approach. Other interesting candidates for strange
dibaryons include the $N\Omega$ ($uudsss$ or $uddsss$) with $J^{\pi}=2^{+}$
and $I=1/2$, and the $\Omega\Omega$ ($ssssss$ ) with $J^{\pi}=0^{+}$ and $I=0$
[7, 8, 9]. Within various quark models, the $N\Omega$ dibaryon has been
suggested to be a bound state [10, 11] or quasibound state [12]. In addition,
the $\Omega\Omega$ dibaryon was predicted to be a bound state by a chiral
quark model [13, 14, 15, 16, 17], while a weak repulsive $\Omega$-$\Omega$
interaction was suggested by other models [10, 18]. Recent lattice QCD studies
using the time-dependent HAL QCD method have shown the $N\Omega(^{5}S_{2})$
potential to be attractive at all distances, producing a quasi-bound state
near unitarity [19]. The di-Omega ($\Omega\Omega$)(${}^{1}S_{0}$) was also
thoroughly examined using lattice QCD data analyzed with the HAL QCD method
and the results show that the $\Omega\Omega$ system has an overall attraction
and is located near the unitary regime [20]. The $\Omega\Omega$, being the
most strange dibaryon with $S=-6$, is of great interest in QCD [20].

Experimentally, high-energy nuclear collisions provide a unique tool to
extract information on nucleon-hyperon and hyperon-hyperon interactions via
measuring the hadron-pair correlation functions [21, 22, 23]. The STAR
Collaboration has measured the $p$-$\Omega$ correlation function [24] in Au+Au
collisions at $\sqrt{s_{NN}}=200$ GeV and the results suggest a bound state of
$p$-$\Omega$ with a binding energy of $27$ MeV. The ALICE Collaboration’s
measurements in p+p collisions $\sqrt{s_{NN}}=13$ TeV indicated that the
correlation functions predicted by lattice QCD underestimate the data,
warranting further investigations to conclusively ascertain the existence of
this bound state [25].

Exploration of the production of $p\Omega$ and $\Omega\Omega$ dibaryons in
high-energy nuclear collisions may offer insights into their nature [26, 27,
28, 29]. The production of $p\Omega$ and $\Omega\Omega$ dibaryons in high
energy heavy-ion collisions has been investigated within the coalescence model
by assuming the $p\Omega$ and $\Omega\Omega$ dibaryons are bound molecular
states of $p$-$\Omega$ (${}^{5}\rm{S}_{2}$) and $\Omega$-$\Omega$
(${}^{1}\rm{S}_{0}$), respectively. For the $p\Omega$ and $\Omega\Omega$
dibaryons, the Pauli exclusion principle does not operate among their valence
quarks, and these dibaryons can in principle be compact six-quark states [3,
4, 30]. In high-energy heavy-ion collisions, the coalescence production of
clusters generally depends on the internal structures of the hadrons [27, 28].
Therefore, it is interesting to probe the internal structure effects on the
production of $p\Omega$ (${}^{5}S_{2}$) and $\Omega\Omega$ (${}^{1}S_{0}$) in
these collisions by treating them as either molecular or six-quark states.
This is our main motivation of this work.

In the present study, we compute the production of the $p\Omega$ and
$\Omega\Omega$ exotic states in central Pb+Pb collisions at
$\sqrt{s_{NN}}$=2.76 TeV using a covariant coalescence model with a blast-
wave-like parametrization for the phase-space configurations of constituent
particles at freeze-out. Specifically, we calculate the yields of $p\Omega$
and $\Omega\Omega$ through hadron coalescence in the case of molecular states
and through quark coalescence for six-quark states. We find a distinct
difference in the centrality dependence of the yield ratios (scaled by the
corresponding yield of $\Omega$) between the molecular and multi-quark states,
providing a conceivable tool to differentiate between these two structures.

##  2 Covariant Coalescence Model

The main feature of the coalescence model [31, 32, 33] is that the coalescence
probability depends on the details of the phase space configurations of the
constituent particles at freeze-out as well as the statistical weight and wave
function of the coalesced cluster. At this point, we would like to mention
that the details of cluster wavefunction are of no relevance in the thermal
model [34, 35, 36, 37, 38, 39] for cluster production. For particle production
at mid-rapidity in central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV at LHC
considered here, we assume a longitudinal boost-invariant expansion for the
constituent particles which are emitted from a freezeout hypersurface
$\Sigma^{\mu}$, and the Lorentz invariant one-particle momentum distribution
is then given by

$\displaystyle
E\frac{\text{d}^{3}N}{\text{d}^{3}p}=\frac{\text{d}^{3}N}{p_{T}\text{d}p_{T}\text{d}\phi_{p}\text{d}y}=\int\text{d}^{4}xS(x,p),$
(1)

where $p^{\mu}$ is the four-momentum of the emitted particle,$\phi_{p}$ is
azimuthal direction of the emitted particle,
$p_{T}=\sqrt{p_{x}^{2}+p_{y}^{2}}$ is the transverse momentum and $S(x,p)$ is
the emission function which reads

$\displaystyle
S(x,p)\text{d}^{4}x=m_{T}\text{cosh}(\eta-y)f(x,p)J(\tau)\text{d}\tau\text{d}\eta
r\text{d}r\text{d}\phi_{s}.$ (2)

In above expressions, we use longitudinal proper time
$\tau=\sqrt{t^{2}-z^{2}}$, space-time rapidity
$\eta=\frac{1}{2}\text{ln}\frac{t-z}{t+z}$, cylindrical coordinates ($r$,
$\phi_{s}$), rapidity $y=\frac{1}{2}\ln(\frac{E+p_{z}}{E-p_{z}})$, transverse
momentum ($p_{T},\phi_{p}$), and transverse mass
$m_{T}=\sqrt{m^{2}+p_{T}^{2}}$. In the relativistic heavy ion collisions, the
statistical distribution function $f(x,p)$ is given by

$\displaystyle f(x,p)=g(2\pi)^{-3}[\exp(p^{\mu}u_{\mu}/kT)/\xi\pm 1]^{-1}$ (3)

with $g$ being spin degeneracy factor, $\xi$ the fugacity, $u_{\mu}$ the four-
velocity of a fluid element in the fireball, and $T$ the local temperature.
The detailed information can be found in Ref. [40, 41]. The freeze-out process
generally occurs gradually over a period of time, and here assume the freeze-
out time follows a Gaussian distribution [42]

$\displaystyle
J(\tau)=\frac{1}{\Delta\tau\sqrt{2\pi}}\exp[-\frac{(\tau-\tau_{0})^{2}}{2(\Delta\tau)^{2}}]$
(4)

with a mean value $\tau_{0}$ and a dispersion $\Delta\tau$. The transverse
rapidity distribution of the fluid element in the fireball is parametrized as
$\rho=\rho_{0}r/R_{0}$ with $\rho_{0}$ being the maximum transverse rapidity
and $R_{0}$ the transverse radius of the fireball. The phase-space freeze-out
configuration of the constituent particles are thus determined by six
parameters, i.e., $T$, $\rho_{0}$, $R_{0}$, $\tau_{0}$, $\Delta\tau$ and
$\xi$.

The cluster production probability can be calculated through the overlap of
the cluster Wigner function with the constituent particle phase-space
distribution at freeze-out. If $M$ particles are coalesced into a cluster, the
invariant momentum distribution of the cluster can be obtained as

$\displaystyle E\frac{\text{d}^{3}N_{c}}{\text{d}^{3}P}$ $\displaystyle=$
$\displaystyle
Eg_{c}\int\bigg{(}\prod_{i=1}^{M}\frac{\text{d}^{3}p_{i}}{E_{i}}\text{d}^{4}x_{i}S(x_{i},p_{i})\bigg{)}\times$
$\displaystyle\rho_{c}^{W}(x_{1},...,x_{M};p_{1},...,p_{M})\delta^{3}(\mathbf{P}-\sum_{i=1}^{M}\mathbf{p_{i}}),$

where $N_{c}$ is the cluster multiplicity, $E$ ($\mathbf{P}$) is its energy
(momentum), $g_{c}$ is the coalescence factor [32] which can be expressed as
$g_{c}=\frac{2S+1}{2^{M}3^{M}}$ including the spin and color degrees of
freedom in the quark coalescence, but
$g_{c}=\frac{2S+1}{(2s_{1}+1)(2s_{2}+1)...(2s_{M}+1)}$ in the hadron
coalescence including the spin degrees of freedom, where $s_{M}$ is spin of
particle, and $S$ is spin of cluster. $\rho_{c}^{W}$ is the Wigner function
and $\delta$-function is adopted to ensure momentum conservation. In this
study, the harmonic oscillator wave functions are assumed and the cluster
Wigner function is expressed as

$\displaystyle\rho_{c}^{W}(x_{1},...,x_{M};p_{1},...,p_{M})$
$\displaystyle\quad=\rho^{W}(q_{1},\cdot\cdot\cdot,q_{M-1},k_{1},\cdot\cdot\cdot,k_{M-1})$
$\displaystyle\quad=8^{M-1}\exp[-\sum_{i=1}^{M-1}(q_{i}^{2}/\sigma_{i}^{2}+\sigma_{i}^{2}k_{i}^{2})],$
(6)

where
$\mu_{i-1}=\frac{i}{i-1}\frac{m_{i}\sum_{k=1}^{i-1}m_{k}}{\sum_{k=1}^{k=i}m_{k}},(i\geq
2)$ is the reduced mass in the center-of-mass frame,
$\sigma_{i}^{2}=(\mu_{i}w)^{-1}(1\leq i\leq M-1)$, and $w$ is the harmonic
oscillator frequency. The detailed information about the coordinate
transformation from $(x_{1},...,x_{M}),(p_{1},...,p_{M})$ to relative
coordinates $(q_{1},...,q_{M-1}),(k_{1},...,k_{M-1})$ can be found in Ref.
[41]. If the masses of constituent particles are equal, we then have
$\mu_{i}=m$. The mean-square radius is given by

$\displaystyle\langle
r^{2}_{M}\rangle=\frac{3}{2Mw}[\sum_{i=1}^{M}\frac{1}{m_{i}}-\frac{M}{\sum_{i=1}^{M}m_{i}}].$
(7)

The integral (LABEL:Eq:Coal) can be directly calculated through multi-
dimensional numerical integration by Mote-Carlo method [41, 43]. It should be
emphasized that since the constituent particles may have different freeze-out
time, in the numerical calculation, the particles that freeze out earlier are
allowed to propagate freely until the time when the last particle in the
cluster freezes out in order to make the coalescence at equal time [41, 44,
45].

##  3 Results and Discussion

###  3.1 Molecular states from hadron coalescence

We first consider the $p\Omega$ and $\Omega\Omega$ dibaryons as molecular
states, and their production in central Pb+Pb collisions at
$\sqrt{s_{NN}}=2.76$ TeV can be described by $p$-$\Omega$ and
$\Omega$-$\Omega$ coalescence, respectively, in the covariant coalescence
model. The freeze-out configurations of nucleons (denoted FOPb-$p$) and
$\Omega$ (denoted FOPb-$\Omega$) for central Pb+Pb collisions at
$\sqrt{s_{NN}}=2.76$ TeV have been obtained by fitting the experimental
spectra [46, 47, 48, 49] of protons, $\Omega$, deuterons, and 3He [40]. The
parameter values of FOPb-$p$ and FOPb-$\Omega$ are summarized in Tab. 1, and
the details can be found in Refs.[41, 50, 51]. The freeze-out temperature $T$
of nucleons [51] can be extracted from measured charged particle spectra
expressed as:

$\displaystyle
T=T_{0}+T_{1}\left[1+(q-1)\times\frac{\text{d}N_{ch}/\text{d}\eta}{M}\right]^{-\frac{1}{q-1}},$
(8)

by using the four parameters $T_{0}=80.6\pm 31.0$ MeV, $T_{1}=83.0\pm 46.9$
MeV, $q=3.33\pm 3.25$, and $M=67.3\pm 76.3$ after taking into account the
errors. The values of $T$ in different centralities are shown in Tab. 1
(FOPb-$p$). According to Ref. [51], the values of $R_{0}$ and $\tau_{0}$ of
nucleons in different centralities can be calculated by using the equal source
volume $2\pi^{1.5}R^{3}=\pi R_{0}^{2}\tau_{0}$ and the fixed ratio
$R_{0}/\tau_{0}=19.7/15.5$, where $R_{0}$ and $\tau_{0}$ are the same as
“FOPb-$N$” in Ref. [50]. The source radius $R$ can be calculated from

$\displaystyle R=\frac{(3N_{n})^{1/3}}{[4C_{1}(mT)^{3/2}]^{1/3}},$ (9)

where $N_{n}$ is the neutron number which is the same as the proton number
[52] in collisions at the LHC energies, the factor $C_{1}=4.0\times 10^{-3}$
is obtained from the yield ratio d/p [51]. The $R_{0}$ and $\tau_{0}$ in
different centralities are shown in Tab. 1 (FOPb-$p$). With the values
$R_{0}$, $\tau_{0}$ and the temperature $T$, the other parameters can be
obtained by fitting the $p$ spectra [52] in the corresponding centrality and
shown in Tab. 1 (FOPb-$p$). Similar to “FOPb-p”, the $T$, $R_{0}$ and
$\tau_{0}$ of FOPb-$\Omega$ have same values with FOPb-$p$, then the other
parameters of FOPb-$\Omega$ can be obtained by fitting the $\Omega$ spectra
[49] in the corresponding centrality and shown in Tab. 1 (FOPb-$\Omega$).

The cluster yield also depends on the cluster size [50, 53, 54], we thus
calculate the yields of $p\Omega$ and $\Omega\Omega$ as a function of their
root-mean-square radii ($r_{\rm rms}$). Fig. 1 shows the $p_{T}$-integrated
yield in the mid-rapidity region ($-0.5\leq y\leq 0.5$) versus $r_{\rm rms}$
in the range of $0.5\sim 5$ fm through hadron coalescence (solid lines) in
central collision parameters. For FOPb-p and FOPb-$\Omega$ in centrality
$0\sim 10\%$, the yield (d$N/$d$y$) of $p\Omega$ ranges from $7.72\times
10^{-4}$ to $1.82\times 10^{-3}$, and the yield of $\Omega\Omega$ ranges from
$7.35\times 10^{-7}$ to $2.15\times 10^{-6}$.

In lattice QCD calculations, the binding energy of $N\Omega$ ($B_{N\Omega}$)
is 1.54 MeV and 2.46 MeV with and without taking into account the Coulomb
attraction [19], and the corresponding root mean square distance is
$\sqrt{\langle r_{M}^{2}\rangle}=3.77$ fm and $\sqrt{\langle
r_{M}^{2}\rangle}=3.24$ fm. The binding energy of $\Omega\Omega$ is
$B_{\Omega\Omega}=1.6$ MeV or $B_{\Omega\Omega}=0.7$ MeV with/without Coulomb
attraction [20] and the corresponding root mean square distance is about
$\sqrt{\langle r_{M}^{2}\rangle}=3.28$ fm and $\sqrt{\langle
r_{M}^{2}\rangle}=4.12$ fm. So we assume the binding energy of $p\Omega$ is
2.46 MeV, the binding energy of $\Omega\Omega$ is 0.7 MeV. The corresponding
$r_{\rm rms}$ of $p\Omega$ and $\Omega\Omega$ are 3.24 fm and 4.12 fm [28]
which are shown by diamonds in the solid lines in Fig. 1. The yields of
$p\Omega$ and $\Omega\Omega$ are $1.27\times 10^{-3}$ at $r_{\rm rms}$ = 3.24
fm and $9.87\times 10^{-7}$ at $r_{\rm rms}=4.12$ fm, respectively.

Table 1: Parameters of the blast-wave-like analytical parametrization for the
phase-space configurations at freeze-out for (anti-)nucleon (FOPb-$p$) [41],
(anti-)$\Omega$ (FOPb-$\Omega$) [50] for Pb+Pb central collisions at
$\sqrt{s_{NN}}=2.76$ TeV in different centralities.

FOPb-$p$ |  |  |  |  |  |
---|---|---|---|---|---|---
Centrality | T(MeV) | $\rho_{0}$ |  $R_{0}$(fm) |  $\tau_{0}$(fm/c) | $\Delta\tau$ | $\xi_{H}$
0-10% | 95.94 | 1.28 | 19.53 | 15.36 | 1.0 | 38.55
10-20% | 98.77 | 1.25 | 16.90 | 13.30 | 1.0 | 29.55
20-40% | 103.74 | 1.21 | 13.68 | 10.77 | 1.0 | 19.00
40-60% | 108.30 | 1.12 | 9.90 | 7.79 | 1.0 | 13.28
60-80% | 114.42 | 1.01 | 6.37 | 5.01 | 1.0 | 8.50
FOPb-$\Omega$ |  |  |  |  |  |
Centrality | T(MeV) | $\rho_{0}$ |  $R_{0}$(fm) |  $\tau_{0}$(fm/c) | $\Delta\tau$ | $\xi_{H}$
0-10% | 95.94 | 1.08 | 19.53 | 15.36 | 1.0 | 643.19
10-20% | 98.77 | 1.10 | 16.90 | 13.30 | 1.0 | 375.80
20-40% | 103.74 | 1.04 | 13.68 | 10.77 | 1.0 | 174.88
40-60% | 108.30 | 0.97 | 9.90 | 7.79 | 1.0 | 81.08
60-80% | 114.42 | 0.95 | 6.37 | 5.01 | 1.0 | 19.73

Table 2: Parameters of the blast-wave-like analytical parametrization for the
phase-space configurations at freeze-out for light quarks (FOPb-Q) [40] for
Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV in different centralities.

| Centrality | T(MeV) | $\rho_{0}$ |  $R_{0}$(fm) |  $\tau_{0}$(fm/c) | $\Delta\tau$ | $\xi_{u}$ | $\xi_{s}$
---|---|---|---|---|---|---|---|---
| $0\sim 10\%$ | 154 | 1.08 | 13.6 | 11.0 | 1.3 | 1.02 | 0.89
| $10\sim 20\%$ | 154 | 1.08 | 12.0 | 9.7 | 1.3 | 1.02 | 0.89
| $20\sim 40\%$ | 154 | 1.08 | 9.9 | 8.0 | 1.3 | 1.02 | 0.89
| $40\sim 60\%$ | 157 | 1.03 | 7.3 | 5.94 | 1.3 | 1.02 | 0.89
| $60\sim 80\%$ | 160 | 0.95 | 4.8 | 3.9 | 1.3 | 1.13 | 0.85

Figure 1: The yields of $p\Omega$ (top) and $\Omega\Omega$ (bottom) versus
root-mean-square radii through hadron coalescence and quark coalescence by
using the parameters of Table 1 and 2 in centrality $0-10\%$.

###  3.2 Six-quark states from quark coalescence

We now model $p\Omega$ and $\Omega\Omega$ as compact six-quark states, where
$p\Omega$ is a $uddsss$ bound state, and $\Omega\Omega$ is a $ssssss$ bound
state. Their yields are calculated using the quark coalescence model. By
fitting the measured spectra of $p(udd)$ and $\Omega^{-}(sss)$, as detailed in
our previous work [40], we obtain phase-space distributions of $u,d,s$ quarks
and their anti-partners at freeze-out. We assume $u$ and $d$ quarks have a
mass of 300 MeV, while the $s$-quark mass is 500 MeV. The freeze-out
parameters are listed in Tab. 2 (FOPb-Q). For instance, in the $0-10\%$
centrality, $T=154$ MeV, $\rho_{0}=1.08$, $R_{0}=13.6$ fm, $\tau_{0}=11.0$
fm/$c$, $\Delta\tau_{0}=1.3$ fm/$c$, $\xi_{u}=1.02$, and $\xi_{s}=0.89$ are
the temperature, transverse flow rapidity, transverse radius, mean
longitudinal proper freeze-out time, proper time dispersion, and fugacities
for $u$ and $s$ quarks, respectively (refer to Ref. [40] for more details).
Using the quark freeze-out configuration FOPb-Q, we calculate the yields of
$p\Omega$ and $\Omega\Omega$ through six-quark coalescence.

The $r_{\rm rms}$ dependence of the yields of $p\Omega$ and $\Omega\Omega$ in
$0-10\%$ centrality are depicted as dashed lines in Fig. 1. For quark
coalescence, the yield of $p\Omega$ varies from $3.89\times 10^{-5}$ to
$1.54\times 10^{-3}$, and for $\Omega\Omega$, from $2.49\times 10^{-7}$ to
$8.80\times 10^{-6}$. We assume the harmonic oscillator frequencies for
$p\Omega$ and $\Omega\Omega$ are $\omega_{s}=78$ MeV [40]. Consequently, the
calculated $r_{\rm rms}$ values for $p\Omega$ and $\Omega\Omega$ are 1.30 fm
and 1.12 fm, respectively, indicated by stars in Fig. 1, with corresponding
yields of $6.77\times 10^{-4}$ and $3.45\times 10^{-6}$.

In Fig. 1, for $p\Omega$, it is seen that its estimated yield from hadron
coalescence is 1.87 times of that from quark coalescence. For smaller $r_{\rm
rms}$, the yield from hadron coalescence is about one order larger than that
from quark coalescence. This discrepancy is attributed to the strong decays in
the hadron coalescence significantly contributing to proton multiplicities,
thereby enhancing $p\Omega$ yields through $p$-$\Omega$ coalescence [50].
Conversely, for $\Omega\Omega$, the yields predicted by quark coalescence are
marginally lower than hadron coalescence at $r_{\rm rms}<0.97$ fm but slightly
higher at $r_{\rm rms}>0.97$ fm. This difference is mainly due to that
$\Omega$ baryon consists of three valence strange quarks or anti-strange
quarks. So the strong decays rarely contribute to the $\Omega\Omega$
multiplicities, and the predicted yields by the hadron coalescence are near to
the yields of quark coalescence which are only produced directly from the
quark coalescence. As a result, it is difficult to distinguish whether the
$p\Omega$ and $\Omega\Omega$ are molecular states or compact six-quark states
merely by their yields.

Figure 2: The yields ratio of $p\Omega/\Omega$ (top) and $\Omega\Omega/\Omega$
(bottom) versus centrality. The yields of $p\Omega$ and $\Omega\Omega$ are
calculated through hadron coalescence and quark coalescence by using the
parameters of Table 1 and 2 with the fixed root-mean-square radii. The yields
of $\Omega$ comes form experiment data [46, 47, 48, 49].

Figure 2 shows the centrality dependence of yield ratios
$N_{p\Omega}/N_{\Omega}$ (top) and $N_{\Omega\Omega}/N_{\Omega}$ (bottom),
with fixed $r_{\rm rms}=1.0$ fm for quark coalescence and $r_{\rm rms}=4.0$ fm
for hadron coalescence. In Fig. 2 (a), the yellow circles represent the
$p\Omega$ yield ratio from hadron coalescence to $\Omega$ from experimental
data [49], showing a rapid decrease from central to peripheral collisions. The
black squares show the yield ratio of $p\Omega$ from 6-quark coalescence to
$\Omega$ from experimental data, with a milder trend change from central to
peripheral collisions. Fig. 2 (b) displays the $N_{\Omega\Omega}/N_{\Omega}$
ratio, indicating a rapid drop from central to peripheral collisions in hadron
coalescence but only a mild change in multi-quark coalescence. These different
centrality dependencies provide a promising tool for distinguishing molecular
and multi-quark states.

Figure 3: Same with Fig. 2 but in the predicted root-mean-square radii.

Figure 3 displays the yield ratios of $N_{p\Omega}/N_{\Omega}$ (top) and
$N_{\Omega\Omega}/N_{\Omega}$ (bottom) with the estimated root-mean-square
radii, i.e., $r_{\rm rms}=1.3$ fm for multi quark coalescence of $p\Omega$
($r_{\rm rms}=1.12$ fm for $\Omega\Omega$), and $r_{\rm rms}=3.24$ fm for
hadron coalescence of $p\Omega$ ($r_{\rm rms}=4.12$ fm for $\Omega\Omega$). It
is seen that the $N_{p\Omega}/N_{\Omega}$ and $N_{\Omega\Omega}/N_{\Omega}$
drops rapidly in hadron coalescence but remains same trend in multi-quark
coalescence as the centrality increasing. So the obviously different
centrality dependence of yield ratio for $p\Omega$ and $\Omega\Omega$ is a
unique way to distinguish whether a structure is molecular or multi-quark
state. The difference in yield ratio between molecular and multi-quark state
increases with the greater difference in root-mean-square radii.

##  4 Conclusion

In the present study, we have analyzed the production of $p\Omega$ and
$\Omega\Omega$ in central Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV,
employing a covariant coalescence model with blast-wave-like parametrization
for constituent particle phase-space configurations at freeze-out. We have
considered two scenarios: $p\Omega$ and $\Omega\Omega$ as either molecular or
compact six-quark states. For the former, the phase-space configurations of
constituent particles are determined from light nuclei and hypernuclei
spectra; for the latter, they are determined from experimental data on hadrons
such as $p$, $\Lambda$, $\phi$, $\Xi$, $\Omega^{-}$.

Our results indicate that in 0-10% centrality, $p\Omega$ yields range from
$7.72\times 10^{-4}$ to $1.82\times 10^{-3}$ in the hadron coalescence and
$3.89\times 10^{-5}$ to $1.54\times 10^{-3}$ in the quark coalescence, with
$r_{\rm rms}$ spanning 0.5 to 5.0 fm. $\Omega\Omega$ yields vary from
$7.35\times 10^{-7}$ to $2.15\times 10^{-6}$ in the hadron coalescence and
$2.49\times 10^{-7}$ to $8.80\times 10^{-6}$ in the quark coalescence. The
STAR Collaboration’s measurement of 4He yields ($8.6\times 10^{-9}$ [55]) in
$\sqrt{s_{NN}}=$200 GeV central Au+Au collisions suggests the feasibility of
experimentally measuring $p\Omega$ and $\Omega\Omega$ yields. Estimated root-
mean-square radii for $p\Omega$ are 3.24 fm (hadron coalescence) and 1.30 fm
(multi-quark coalescence), with yields of $1.27\times 10^{-3}$ and $6.77\times
10^{-4}$, respectively. For $\Omega\Omega$, these radii are 4.12 fm (hadron
coalescence) and 1.12 fm (multi-quark coalescence), with yields of $9.87\times
10^{-7}$ and $3.45\times 10^{-6}$. Hadron coalescence predicts $p\Omega$
yields 1.87 times higher than multi-quark coalescence, and $\Omega\Omega$
yields in quark coalescence are 3.49 times those from hadron coalescence.

The $p\Omega$ and $\Omega\Omega$ yields in different centralities have been
calculated by using the corresponding collision parameters. The yield ratios
of $N_{p\Omega}/N_{\Omega}$ and $N_{\Omega\Omega}/N_{\Omega}$ going from
central to peripheral collisions decrease for the hadron coalescence but
changes mildly for the multi-quark coalescence. This distinct centrality
dependence of yield ratios for molecular and multi-quark states offers a
promising way for distinguishing these two states. Our results thus suggest
that the measurements on the production of $p\Omega$ and $\Omega\Omega$ in
relativistic heavy-ion collisions can shed light on their internal structures
and advance our understanding of the nature of strong interaction.

## Acknowledgements

We thank Kai-Xuan Cheng and Yu-Ting Wang for reading the manuscript. This work
was supported by the National Key Research and Development Project of China
under Grant No. 2022YFA1602303, the National Natural Science Foundation of
China under Grant No. 12105079, No. 12235010, No. 12375121, and the National
SKA Program of China No. 2020SKA0120300.

## References

  * [1] H. Clement, On the History of Dibaryons and their Final Observation, Prog. Part. Nucl. Phys. 93 (2017) 195.  arXiv:1610.05591, doi:10.1016/j.ppnp.2016.12.004.
  * [2] S. Cho, et al., Exotic hadrons from heavy ion collisions, Prog. Part. Nucl. Phys. 95 (2017) 279–322.  arXiv:1702.00486, doi:10.1016/j.ppnp.2017.02.002.
  * [3] M. Oka, Flavor Octet Dibaryons in the Quark Model, Phys. Rev. D 38 (1988) 298\.  doi:10.1103/PhysRevD.38.298.
  * [4] A. Gal, Meson assisted dibaryons, Acta Phys. Polon. B 47 (2016) 471.  arXiv:1511.06605, doi:10.5506/APhysPolB.47.471.
  * [5] H. C. Urey, F. G. Brickwedde, G. M. Murphy, A hydrogen isotope of mass 2, Phys. Rev. 39 (1932) 164–165.  doi:10.1103/PhysRev.39.164.
  * [6] R. L. Jaffe, Perhaps a Stable Dihyperon, Phys. Rev. Lett. 38 (1977) 195–198, [Erratum: Phys.Rev.Lett. 38, 617 (1977)].  doi:10.1103/PhysRevLett.38.195.
  * [7] J. T. Goldman, K. Maltman, G. J. Stephenson, Jr., K. E. Schmidt, F. Wang, STRANGENESS -3 DIBARYONS, Phys. Rev. Lett. 59 (1987) 627.  doi:10.1103/PhysRevLett.59.627.
  * [8] F. Etminan, H. Nemura, S. Aoki, T. Doi, T. Hatsuda, Y. Ikeda, T. Inoue, N. Ishii, K. Murano, K. Sasaki, Spin-2 $N\Omega$ dibaryon from Lattice QCD, Nucl. Phys. A 928 (2014) 89–98.  arXiv:1403.7284, doi:10.1016/j.nuclphysa.2014.05.014.
  * [9] K. Morita, A. Ohnishi, F. Etminan, T. Hatsuda, Probing multistrange dibaryons with proton-omega correlations in high-energy heavy ion collisions, Phys. Rev. C 94 (3) (2016) 031901, [Erratum: Phys.Rev.C 100, 069902 (2019)].  arXiv:1605.06765, doi:10.1103/PhysRevC.94.031901.
  * [10] F. Wang, J. l. Ping, G. h. Wu, L. j. Teng, J. T. Goldman, Quark delocalization, color screening and dibaryons, Phys. Rev. C 51 (1995) 3411.  arXiv:nucl-th/9512014, doi:10.1103/PhysRevC.51.3411.
  * [11] H. Huang, J. Ping, F. Wang, Further study of the $N\Omega$ dibaryon within constituent quark models, Phys. Rev. C 92 (2015) 065202.  arXiv:1507.07124, doi:10.1103/PhysRevC.92.065202.
  * [12] T. Sekihara, Y. Kamiya, T. Hyodo, $N\Omega$ interaction: meson exchanges, inelastic channels, and quasibound state, Phys. Rev. C 98 (1) (2018) 015205.  arXiv:1805.04024, doi:10.1103/PhysRevC.98.015205.
  * [13] Z. Y. Zhang, Y. W. Yu, P. N. Shen, L. R. Dai, A. Faessler, U. Straub, Hyperon nucleon interactions in a chiral SU(3) quark model, Nucl. Phys. A 625 (1997) 59–70.  doi:10.1016/S0375-9474(97)00033-X.
  * [14] Z. Y. Zhang, Y. W. Yu, C. R. Ching, T. H. Ho, Z.-D. Lu, Suggesting a di-omega dibaryon search in heavy ion collision experiments, Phys. Rev. C 61 (2000) 065204\.  doi:10.1103/PhysRevC.61.065204.
  * [15] Q. B. Li, P. N. Shen, Z. Y. Zhang, Y. W. Yu, Dibaryon systems in chiral SU(3) quark model, Nucl. Phys. A 683 (2001) 487–509.  arXiv:nucl-th/0009038, doi:10.1016/S0375-9474(00)00458-9.
  * [16] L. R. Dai, Z. Y. Zhang, Y. W. Yu, Structure of di-Omega dibaryon, Chin. Phys. Lett. 23 (2006) 3215–3218.  doi:10.1088/0256-307X/23/12/026.
  * [17] H. Huang, X. Zhu, J. Ping, Reanalysis of the most strange dibaryon within constituent quark models, Phys. Rev. C 101 (3) (2020) 034004.  arXiv:1912.11256, doi:10.1103/PhysRevC.101.034004.
  * [18] M. I. Buchoff, T. C. Luu, J. Wasem, S-wave scattering of strangeness -3 baryons, Phys. Rev. D 85 (2012) 094511.  arXiv:1201.3596, doi:10.1103/PhysRevD.85.094511.
  * [19] T. Iritani, et al., $N\Omega$ dibaryon from lattice QCD near the physical point, Phys. Lett. B 792 (2019) 284–289.  arXiv:1810.03416, doi:10.1016/j.physletb.2019.03.050.
  * [20] S. Gongyo, et al., Most Strange Dibaryon from Lattice QCD, Phys. Rev. Lett. 120 (21) (2018) 212001.  arXiv:1709.00654, doi:10.1103/PhysRevLett.120.212001.
  * [21] L. Adamczyk, et al., $\Lambda\Lambda$ Correlation Function in Au+Au collisions at $\sqrt{s_{NN}}=$ 200 GeV, Phys. Rev. Lett. 114 (2) (2015) 022301.  arXiv:1408.4360, doi:10.1103/PhysRevLett.114.022301.
  * [22] L. Adamczyk, et al., Measurement of Interaction between Antiprotons, Nature 527 (2015) 345–348.  arXiv:1507.07158, doi:10.1038/nature15724.
  * [23] J. Chen, D. Keane, Y. G. Ma, A. Tang, Z. Xu, Antinuclei in Heavy-Ion Collisions, Phys. Rept. 760 (2018) 1–39.  arXiv:1808.09619, doi:10.1016/j.physrep.2018.07.002.
  * [24] J. Adam, et al., The Proton-$\Omega$ correlation function in Au+Au collisions at $\sqrt{s_{NN}}$=200 GeV, Phys. Lett. B 790 (2019) 490–497.  arXiv:1808.02511, doi:10.1016/j.physletb.2019.01.055.
  * [25] S. Acharya, et al., Unveiling the strong interaction among hadrons at the LHC, Nature 588 (2020) 232–238, [Erratum: Nature 590, E13 (2021)].  arXiv:2005.11495, doi:10.1038/s41586-020-3001-6.
  * [26] Z. W. Lin, L. Zheng, Further developments of a multi-phase transport model for relativistic nuclear collisions, Nucl. Sci. Tech. 32 (10) (2021) 113.  arXiv:2110.02989, doi:10.1007/s41365-021-00944-5.
  * [27] S. Cho, et al., Multi-quark hadrons from Heavy Ion Collisions, Phys. Rev. Lett. 106 (2011) 212001.  arXiv:1011.0852, doi:10.1103/PhysRevLett.106.212001.
  * [28] S. Cho, et al., Studying Exotic Hadrons in Heavy Ion Collisions, Phys. Rev. C 84 (2011) 064910.  arXiv:1107.1302, doi:10.1103/PhysRevC.84.064910.
  * [29] S. Zhang, Y. G. Ma, $\Omega$-dibaryon production with hadron interaction potential from the lattice QCD in relativistic heavy-ion collisions, Phys. Lett. B 811 (2020) 135867.  arXiv:2007.11170, doi:10.1016/j.physletb.2020.135867.
  * [30] M. Chen, F. Gong, H. Huang, J. Ping, Dibaryons with strangeness in quark models, Mod. Phys. Lett. A 25 (2010) 1603–1612.  doi:10.1142/S0217732310032913.
  * [31] S. T. Butler, C. A. Pearson, Deuterons from High-Energy Proton Bombardment of Matter, Phys. Rev. Lett. 7 (1961) 69–71.  doi:10.1103/PhysRevLett.7.69.
  * [32] H. Sato, K. Yazaki, On the coalescence model for high-energy nuclear reactions, Phys. Lett. B 98 (1981) 153–157.  doi:10.1016/0370-2693(81)90976-X.
  * [33] L. P. Csernai, J. I. Kapusta, Entropy and Cluster Production in Nuclear Collisions, Phys. Rept. 131 (1986) 223–318.  doi:10.1016/0370-1573(86)90031-1.
  * [34] A. Andronic, P. Braun-Munzinger, J. Stachel, H. Stocker, Production of light nuclei, hypernuclei and their antiparticles in relativistic nuclear collisions, Phys. Lett. B 697 (2011) 203–207.  arXiv:1010.2995, doi:10.1016/j.physletb.2011.01.053.
  * [35] J. Cleymans, K. Redlich, E. Suhonen, Canonical description of strangeness conservation and particle production, Z. Phys. C 51 (1991) 137–141.  doi:10.1007/BF01579571.
  * [36] J. Steinheimer, K. Gudima, A. Botvina, I. Mishustin, M. Bleicher, H. Stocker, Hypernuclei, dibaryon and antinuclei production in high energy heavy ion collisions: Thermal production versus Coalescence, Phys. Lett. B 714 (2012) 85–91.  arXiv:1203.2547, doi:10.1016/j.physletb.2012.06.069.
  * [37] J. Cleymans, S. Kabana, I. Kraus, H. Oeschler, K. Redlich, N. Sharma, Antimatter production in proton-proton and heavy-ion collisions at ultrarelativistic energies, Phys. Rev. C 84 (2011) 054916.  arXiv:1105.3719, doi:10.1103/PhysRevC.84.054916.
  * [38] P. Braun-Munzinger, J. Stachel, Production of strange clusters and strange matter in nucleus-nucleus collisions at the AGS, J. Phys. G 21 (1995) L17–L20.  arXiv:nucl-th/9412035, doi:10.1088/0954-3899/21/3/002.
  * [39] P. Braun-Munzinger, J. Stachel, The quest for the quark-gluon plasma, Nature 448 (2007) 302–309.  doi:10.1038/nature06080.
  * [40] K. J. Sun, L. W. Chen, Production of $\Lambda\Lambda$ and $\overline{\Lambda\text{n}}$ in central Pb+Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV within a covariant coalescence model, Phys. Rev. C 94 (6) (2016) 064908.  arXiv:1607.04037, doi:10.1103/PhysRevC.94.064908.
  * [41] K. J. Sun, L. W. Chen, Production of antimatter 5,6Li nuclei in central Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV, Phys. Lett. B 751 (2015) 272–277.  arXiv:1509.05302, doi:10.1016/j.physletb.2015.10.056.
  * [42] F. Retiere, M. A. Lisa, Observable implications of geometrical and dynamical aspects of freeze out in heavy ion collisions, Phys. Rev. C 70 (2004) 044907\.  arXiv:nucl-th/0312024, doi:10.1103/PhysRevC.70.044907.
  * [43] G. P. Lepage, A New Algorithm for Adaptive Multidimensional Integration, J. Comput. Phys. 27 (1978) 192.  doi:10.1016/0021-9991(78)90004-9.
  * [44] R. Mattiello, H. Sorge, H. Stoecker, W. Greiner, Nuclear clusters as a probe for expansion flow in heavy ion reactions at 10-A/GeV - 15-A/GeV, Phys. Rev. C 55 (1997) 1443–1454.  arXiv:nucl-th/9607003, doi:10.1103/PhysRevC.55.1443.
  * [45] L. W. Chen, C. M. Ko, phi and omega production from relativistic heavy ion collisions in a dynamical quark coalescence model, Phys. Rev. C 73 (2006) 044903\.  arXiv:nucl-th/0602025, doi:10.1103/PhysRevC.73.044903.
  * [46] B. Abelev, et al., Pion, Kaon, and Proton Production in Central Pb–Pb Collisions at $\sqrt{s_{NN}}=2.76$ TeV, Phys. Rev. Lett. 109 (2012) 252301\.  arXiv:1208.1974, doi:10.1103/PhysRevLett.109.252301.
  * [47] B. B. Abelev, et al., $K^{0}_{S}$ and $\Lambda$ production in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV, Phys. Rev. Lett. 111 (2013) 222301.  arXiv:1307.5530, doi:10.1103/PhysRevLett.111.222301.
  * [48] B. B. Abelev, et al., $K^{*}(892)^{0}$ and $phi(1020)$ production in Pb-Pb collisions at $\sqrt{s{NN}}$ = 2.76 TeV, Phys. Rev. C 91 (2015) 024609.  arXiv:1404.0495, doi:10.1103/PhysRevC.91.024609.
  * [49] B. B. Abelev, et al., Multi-strange baryon production at mid-rapidity in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV, Phys. Lett. B 728 (2014) 216–227, [Erratum: Phys.Lett.B 734, 409–410 (2014)].  arXiv:1307.5543, doi:10.1016/j.physletb.2014.05.052.
  * [50] K. J. Sun, L. W. Chen, Antimatter ${}^{4}_{\Lambda}$H Hypernucleus Production and the ${}^{3}_{\Lambda}$H/3He Puzzle in Relativistic Heavy-Ion Collisions, Phys. Rev. C 93 (6) (2016) 064909.  arXiv:1512.00692, doi:10.1103/PhysRevC.93.064909.
  * [51] K. J. Sun, C. M. Ko, B. Dönigus, Suppression of light nuclei production in collisions of small systems at the Large Hadron Collider, Phys. Lett. B 792 (2019) 132–137.  arXiv:1812.05175, doi:10.1016/j.physletb.2019.03.033.
  * [52] B. Abelev, et al., Centrality dependence of $\pi$, K, p production in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV, Phys. Rev. C 88 (2013) 044910.  arXiv:1303.0737, doi:10.1103/PhysRevC.88.044910.
  * [53] H. Zhang, J. Liao, E. Wang, Q. Wang, H. Xing, Deciphering the Nature of X(3872) in Heavy Ion Collisions, Phys. Rev. Lett. 126 (1) (2021) 012301.  arXiv:2004.00024, doi:10.1103/PhysRevLett.126.012301.
  * [54] K. X. Cheng, J. Pu, Y. T. Wang, Y. F. Guo, C. W. Ma, Non-frozen process of heavy-ion fusion reactions at deep sub-barrier energies, Nucl. Sci. Tech. 33 (10) (2022) 132.  doi:10.1007/s41365-022-01114-x.
  * [55] H. Agakishiev, et al., Observation of the antimatter helium-4 nucleus, Nature 473 (2011) 353, [Erratum: Nature 475, 412 (2011)].  arXiv:1103.3312, doi:10.1038/nature10079.