# Continuously Adapting Random Sampling (CARS) for Power Electronics Parameter
Design

Dominik Happel, Philipp Brendel, Andreas Rosskopf, Stefan Ditze  Fraunhofer
Institute for Integrated Systems and Device Technology IISB
91058 Erlangen, Germany
Email: {dominik.happel, philipp.brendel, andreas.rosskopf,
<EMAIL_ADDRESS>
###### Abstract

To date, power electronics parameter design tasks are usually tackled using
detailed optimization approaches with detailed simulations or using brute
force grid search grid search with very fast simulations. A new method, named
”Continuously Adapting Random Sampling” (CARS) is proposed, which provides a
continuous method in between. This allows for very fast, and / or large
amounts of simulations, but increasingly focuses on the most promising
parameter ranges. Inspirations are drawn from multi-armed bandit research and
lead to prioritized sampling of sub-domains in one high-dimensional parameter
tensor. Performance has been evaluated on three exemplary power electronic
use-cases, where resulting designs appear competitive to genetic algorithms,
but additionally allow for highly parallelizable simulation, as well as
continuous progression between explorative and exploitative settings.

###### Index Terms:

power electronics design, multi-armed bandit, genetic algorithms, focused
random sampling

##  I Introduction

Power electronics as research field is increasingly drawing the attention of
applied Artificial Intelligence researchers. The Overview of Artificial
Intelligence Applications in Power Electronics by Zhao et al [1], points that
the amount of research papers in this domain has been continuously increasing
over the last years. Additionally, they provide three application categories
corresponding to the lifecycle of power electronics systems: design process,
control, and maintenance.

The review by De León-Aldaco et al [2] notes that “classical” optimization
methods with optimality guarantees are often overwhelmed by the complexity of
power electronic parameter spaces. Instead “metaheuristic” methods aim at
”good enough” solutions, which enables them to effectively explore large
search spaces, but which prevents arguments about solution completeness or
optimality. Among metaheuristic methods, a separation can be performed into
“trajectory-based” and “population-based” methods. Trajectory-based methods
often perform local optimization, based on one solution, to quickly reach some
local optima. Population based methods, which also include Genetic Algorithms
(GAs) [3], on the other hand consider multiple solutions simultaneously, often
leading to better exploration.

A different argument is proposed by Guillod and Kolar [4], who suggest that
fast, vectorized simulations and brute force grid search should be preferred
for magnetics design. Using a priori fixed grid sizes, this will eventually
lead to many samples investigating non useful areas and only a few samples
providing relevant results. This leaves engineers with a rough idea of
interesting areas but requires another more fine-grained iteration to
investigate subtle differences and progressions between relevant samples.

A fast approach, considering all possible parameter settings, while
increasingly focusing on the most relevant parameter ranges would be
desirable. Similar problems are well known in reinforcement learning research,
typically referred to as ”exploration-exploitation dilemma” [5]. Dividing
parameter ranges in a fixed set of sub-domains, power electronics design tasks
can be considered as multi-armed bandit problem, where the goal is to
investigate those sub-domains providing the best results, while ensuring that
no even better sub-domain is ignored. A typical solution to this problem is
upper confidence bounds, which is often superior to other approaches when
attempting to find one most promising sub-domain [5]. However, since power
electronics design is rather concerned with investigating most, if not all
relevant sub-domains, a weighted softmax function [6] can be considered
instead. A similar stochastic sampling approach has been successfully applied
by Schaul et. al. [7] and additionally allows continuous progression between
random sampling (exploration) and greedy focusing on a very few promising sub-
domains (exploitation). Using fast, vectorized array operations [8], sampling
speed is primarily limited by the number of sub-domains on which the softmax
function is computed, which increases exponentially with the number of
sampling parameters.

The sketched procedure presents similarities with Enhanced Continuous Tabu
Search (ECTS) [9], where a search space is first divided into large global
search areas, and later more focused search is performed in sub-domains of the
most promising global search areas. Their comparison with other optimization
approaches showed that their approach is quite competitive, especially for
high dimensional parameter spaces (more than 10 parameters). However, this
approach considers two steps of granularity, where either global
diversification or local intensification steps are performed, which cannot be
easily specified for arbitrary design tasks in advance.

The proposed method, Continuously Adapting Random Sampling (CARS), is based on
the described multi-armed bandit analogy, which provides less granularity than
ECTS [9] but allows arbitrary and especially continuous progression from
exploration as random sampling, to exploitation of one single sub domain, just
as ECTS intensification [7]. Performance of CARS is evaluated on three
exemplary power electronic use-cases and compared to GAs.

The paper is structured as follows, section II describes the CARS method is
detail. A brief overview of GAs is provided in section III. A comparison of
both methods, focused on sample evaluation is performed in section IV, as well
as a CARS runtime study for increasing amounts of parameters and samples. The
three power electronic design tasks are introduced, and comparison results are
presented in section V. Conclusions are drawn in section VI, and finally
limitations and possible directions for further research are discussed in the
last section, VII.

##  II Continuously Adapting Random Sampling

The basic idea of the proposed CARS method is to use a scalar preference
metric for sub-domains of specified parameter ranges, to determine which sub-
domains are ”interesting” (should be investigated more often) and which are
”uninteresting” (should be investigated less often). In reinforcement learning
research this preference measure is referred to as ”reward” [5], for GAs, the
term ”fitness” is established [10]. This paper uses the term ”fitness”, since
GAs are considered for comparison. Fitness computation details for CARS and
GAs are provided in section IV.

Initially, before any CARS sampling can be performed, a specified number of
parameters $N_{param}$, is divided into a number of sub-domains
$N_{subdomain}$, resulting in a tensor with $N_{param}$ dimensions of
$N_{subdomain}$ elements. Inspired by optimistic initialization [5], a fitness
value of $0.75$ is assigned to all sub-domains, which represents an above
average fitness considering fitness normalization to approximately $[0,1]$, as
described in section IV-A.

The actual CARS procedure starts with computing a weighted softmax function
[6] (2) on the flattened sub-domain tensor [8] of size $N_{tensor}$ (1). This
represents probabilities for sampling any sub-domain, identified by its index
$j$. The softmax weighting factor $\alpha$ is used for shifting the softmax
focus from exploration to exploitation (or vice versa).

$N_{tensor}=N_{subdomain}^{N_{param}}.$ (1)

$P(\text{sub-
domain}=j)=\frac{e^{fitness_{j}\cdot\alpha}}{\sum^{N_{tensor}}_{k=1}e^{fitness_{k}\cdot\alpha}}.$
(2)

Since sub-domains represent parameter ranges instead of specific parameter-
values, samples are selected uniformly random within sub-domains.

After simulation, the fitness value (computed according to IV-A) is associated
with its initially selected tensor position $j$ and in case sub-domains have
been investigated multiple times, the respective maximum is chosen. The
assigned fitness values directly affect sub-domain probabilities after softmax
computation, where large fitness values are sampled more frequently, while
initial fitness values (unexplored sub-domains) are still more probable than
very poorly evaluated sub-domains.

The top row in figure 1 shows this procedure for an intermediate iteration,
which starts with assigning previous simulation results to their tensor sub-
domains according to positional index $j$. Then, the softmax function (2) is
computed for the whole tensor, converting fitness values to sub-domain
sampling probabilities. Samples for the next simulation iteration are drawn
uniformly random inside selected sub-domains and used for simulation.

Figure 1: CARS method overview: Top row shows the main process, where fitness
values computed in previous simulations are assigned to their sub-domain
(according to their positional index $j$). Here, variable $Z$ represents a
non-normalized fitness measure. A weighted softmax function (2) is computed,
to provide probabilities for sampling individual sub-domains for the next
simulation iteration (qualitative results for exemplary sub-domains are
shown). The bottom row shows extensions to distribute information to
neighboring areas (pooling) and prevent unnecessary simulations based on their
local fitness estimates (oversampling).

Two additional extensions are implemented, to improve the information gain by
distributing information to adjacent sub-domains, and to prevent unnecessary
simulations due to the fact of random sampling inside sub-domains.

###  II-A Extension: Max-Pooling

Caused by the exponential relation between number of parameters $N_{param}$,
number of sub-domains $N_{subdomain}$, and total number of tensor elements
$N_{tensor}$ (1), even minor problems do not necessarily provide fitness
information for each sub-domain. For example $6$ parameters with $10$ sub-
domains each, results in a tensor of one million ($10^{6}$) sub-domains.

Therefore, a max pooling [11] inspired pooling option is included, to combine
$N_{pool}$ adjacent sub-domains for each parameter to their maximum value,
resulting in a compressed tensor of size:

$N_{pooltensor}=\left(\frac{N_{subdomain}}{N_{pool}}\right)^{N_{param}}.$ (3)

This compressed information is added to the initial $N_{tensor}$ tensor, where
each fitness value is summed with the max-pooled value to which it
contributed. The intuition is to allow some simple interpolation of
neighboring sub-domains (if they are combined in the same max-pooled area),
while simultaneously maintaining differences between sub-domains.
Additionally, the max-pooling operation requires exponentially less operations
compared to averaging operations on the initial tensor, only $N_{pooltensor}$
maximization operations precisely. This is presented in the lower part of
figure 1, where initially $64$ $(8^{2})$ sub-domains are combined to only $4$
$(2^{2})$ max-pooled values and added to the initial tensor.

###  II-B Extension: Oversampling

Besides spreading information to neighboring sub-domains, also sample-
efficiency can be essential. Since not all simulations can be performed with
several thousand simulations per second, a fast approximation of expected
fitness values can reduce simulation requirements. Therefore, fitness values
for prospective samples are estimated using the scikit-learn implementation of
k-nearest neighbor regression [12]. Since nearest neighbor regression does not
build a predictive model, it can estimate fitness values solely based on
existing samples and their distances to samples under consideration.

This extension is especially useful for detailed circuit simulations, where
one circuit simulation consumes much more time than multiple nearest neighbor
fitness estimations. In this case, the CARS method does not just sample an
array of $N_{samples}$, but instead a matrix of $N_{samples}\times
N_{oversampling}$. These samples are estimated using nearest neighbor
regression, and samples for simulation are selected row-wise, where the sample
with highest fitness estimate out of $N_{oversampling}$ samples is chosen. For
example, when selecting two samples with oversampling of three, $2\times 3$
sample candidates are sampled and estimated using nearest neighbor regression.
Based on the fitness estimate, the most promising samples are used for
simulation, in this case each of the two rows selects one out of three
candidates.

###  II-C Heuristics

Since electronics design engineers are not necessarily familiar with
algorithmic details, additional experience driven heuristics have been
implemented to allow configuration of the CARS method based on abstract high-
level-settings.

To ensure fast sampling for potentially large tensors, and to reduce the
”overhead per sample” for computing the softmax function on large tensors,
CARS samples batches of $N_{samples}$ simultaneously. These details are not
necessarily relevant for design engineers, therefore only the number of total
samples $N_{total}$ is required, which is distributed in $N_{iteration}$
iterations of $N_{samples}$ samples each. To prevent extremely odd iteration
numbers, they are preferably assigned as a multiple of $5$.

$\displaystyle N_{iter}=floor(\frac{N_{total}^{\frac{2}{5}}}{5})\cdot 5$
$\displaystyle ifN_{iter}\geq 5$ (4) $\displaystyle
N_{iter}=floor(N_{total}^{\frac{2}{5}})$ $\displaystyle otherwise$ .

$N_{samples}=floor(\frac{N_{total}}{N_{iter}}).$ (5)

Similarly, parameter ranges ($N_{param}$) are by default divided into $9$ sub-
domains ($N_{subdomain}$) and pooling is performed by combining $3$ adjacent
sub-domains of each parameter ($N_{pool}$). By default, the softmax weighting
$alpha$ increases by $1$ for each iteration, starting with initially random
sampling ($\alpha=0$).

Oversampling, settings are also inferred based on the number of parameters
$N_{param}$, representing the dimensionality of the sampling space. If
sufficient samples are available, the number of neighbors for nearest neighbor
regression is computed as (6). Otherwise, all available samples are utilized.

$N_{neighbor}=(N_{param}\cdot 2)+1.$ (6)

The number of additional samples used for oversampling $N_{oversampling}$ is
computed according to:

$N_{oversampling}=N_{param}^{\frac{3}{2}}.$ (7)

This leaves engineers with a need to specify only the number of total samples
$N_{total}$ as well as a set of parameters $N_{param}$ and objectives, while
other settings can be inferred automatically.

##  III Genetic Algorithms

GAs are frequently used to solve many engineering problems [1], [13], [14],
[15] and are therefore used for benchmarking against the proposed CARS method.
GAs are based on the idea of evolving populations consisting of individuals
via the three basic principles (i) mutation, (ii) crossover and (iii)
selection. In the context of this paper, each individual represents a set of
circuit simulation parameters that are randomly mutated, or partially
exchanged (”crossover”) with another individual during evolution. After a
certain number of iterations, the algorithm selects certain individuals based
on their fitness-values, representing the goals of optimization. For this
selection process, we use an adaption of the Non-dominated Sorting Genetic
Algorithm ”NSGA”) [16], that was later referred to as NSGA-II by Deb et al.
[10]. NSGA’s selection is guided by non-dominated sorting, i.e. an individual
is preferably selected, if it is pareto-optimal (there is no other individual
that is better in one objective and not worse in all other objectives) and a
sharing-method promoting diversity in a population and preventing convergence
to single pareto-optimal points, a known caveat of GAs in general, cf. e.g.
[17]. For more details on the specific adaptions introduced by NSGA-II, we
refer to the original work in [16] and [10], respectively. This concept is
even more powerful, when running multiple parallel GAs on separate ”islands”
[3]. On the one hand, this allows efficient usage of High-Performance-
Computing (HPC) resources, to perform simulations for multiple islands in
parallel. On the other hand, it enables the exploration of different evolution
strategies that are ultimately re-combined by selecting the best individuals
among all islands. More Details about the GA setup as well as the used
mutation and crossover techniques are presented in appendix VII-A.

##  IV Method overview

###  IV-A Fitness computation

Both methods, CARS and GAs are designed to prioritize boundary conditions over
objectives [18]. Therefore, a penalty term is computed for parameter settings
that do not meet the boundary condition criteria according to the following
formula, based on the square root of the Canberra distance:

$penalty=\rho*\sqrt{\frac{|value-target|}{|value|+|target|}}.$ (8)

Priorization is achieved by applying a weighting factor $\rho=100$ for CARS
and $\rho=10,000$ for GAs. The square root of the Canberra distance is used,
to ensure relatively large penalties (derived from the distance as in (8))
even for minor deviations. This is shown in figure 2 for a target power value
(solid lines) and for a power range (dashed lines).

Figure 2: Square root Canberra Distance shows a larger increase of its
distance measure, exemplary presented for a ”target” value of 2300W (solid)
and a ”range” of 2600W to 2700W (dashed).

Individual objective fitness values depend on the objective type, where the
measurement value itself is used for ”maximization” objectives. For
”minimization” objectives, it is the negative value. The square root of the
Canberra distance (solid line in Figure 2) with a weight of $\rho=-1$ is used
for ”target” objectives. Finally, a ”minimize range” objective is introduced,
which considers potentially multiple operating points simultaneously. This
fitness measure computes the difference of the largest and the smallest
measurement values for $n$ operating points, according to (9).

$\begin{gathered}range=-max(meas_{1},...,meas_{n})\\\
+min(meas_{1},...,meas_{n}).\end{gathered}$ (9)

CARS considers a scalar fitness value for each sub-domain; therefore, the mean
fitness is computed for all objectives, which is subtracted by the mean
penalty for all boundary conditions. In contrast, GAs consider objectives
individually and therefore the penalty term for each individual boundary
condition is subtracted from all objective-fitness-values independently.

Since CARS is performed in batches, it offers the possibility to automatically
scale fitness values for objectives or penalties for boundary conditions to
approximately $[0,1]$ based on values encountered during the first batch.
Therefore, minimum, and maximum values of the first batch $vals_{first}$ are
stored and used for normalizing current batches $vals_{newest}$ (10). This
normalization is computed for every boundary condition and objective
independently, while all operating points are considered simultaneously. To
prevent overflows caused by extremely large softmax exponents, the resulting
scalar fitness value is normalized as well. However, the maximum fitness value
is usually slightly larger than $1.0$, since the overall best samples are
seldom discovered during the first batch.

$normalized=\frac{vals_{newest}-min(vals_{first})}{max(vals_{first})-min(vals_{first})}.$
(10)

###  IV-B CARS runtime

Since CARS aims to improve on brute force grid search for very fast
simulations [4], the CARS runtime is investigated for batches of up to one
million samples. This includes assigning previous simulations to their
respective sub-domain, softmax probability computation, max-pooling, and
fitness computation, for different amounts of parameters ($N_{param}$) and
batch sizes ($N_{samples}$). The simulation time and the oversampling
extension are excluded for this test, since they would scale linearly with the
number of samples, and nearest neighbor estimation provides little benefit
when performing very fast simulations. The simulation itself could also be
distributed to multiple HPC machines, allowing highly parallelizable
execution.

For each test, $20$ consecutive sampling iterations, performed on one core of
an Intel Xeon Gold 6134 Server CPU, are investigated and reported with mean
(solid line), as well as min and max values (shaded areas) in figure 3. As
expected, especially for many parameters, larger batch-sizes show less
distinct increases in runtime, since the essential softmax operation on the
very large sub-domain tensors needs to be computed only once. For $9$
parameters, the resulting 9-dimensional tensor, with by default $9$ sub-
domains each, contains more than $387$ million sub-domains ($9^{9}$). These
results suggest CARS can serve fast simulations with thousands or even one
million samples in terms of a few seconds. However, especially for many
parameters it might be beneficial to reduce the default number of sub-domains
per parameter, where a reduction to just $8$ sub-domains would reduce the
tensor size to about $134$ million ($8^{9}$).

Figure 3: CARS runtime for different batch-sizes and numbers of sampling
parameters, with each parameter divided into $9$ sub-domains. The CARS method
is able to suggest up to one million samples in the time span of a few
seconds, allowing focusing sampling in combination with very fast simulations.

##  V Power Electronic Use-Cases

To determine the performance of the proposed CARS method in section II, a
parameter design evaluation for three exemplary power electronic use cases is
performed, and compared to GAs. A first proof-of-concept application is a
textbook boost converter V-A, the second is a reproduction of a previous LLC
Resonant Converter optimization using GAs [3] V-B, and the third is a recent
research question for a LLC Resonant Converter with large power range V-C.
These use cases are modelled as LTspice netlists and simulations are performed
using the SchoKI framework [19] for parallel computation of LTspice [20]
simulations.

###  V-A Boost Converter

The circuit for the textbook Boost Converter is shown in Figure 4. Parameters,
objectives and boundary conditions are listed in Table I. This application
includes only one operating point and a total of $5,000$ samples are
investigated for each method. $log$ indicates sampling in logarithmic scale.
Further details on measurements are provided in the appendix VII-B. Note that
this first proof of concept example does not include any component loss
models.

Figure 4: Schematic of Boost Converter. Parameters
---
Name | Type | Values | Operating Points
C1 | log | [1e-9, 1e-3] | 1
L1 | log | [1e-6, 100e-3] | 1
fsw | log | [100, 1e6] | 1
Objectives
Name | Type | Values | Operating Points
vmean | target | 12 | 1
eff_tot | max | - | 1
Boundary Conditions
Name | Type | Values | Operating Points
vmean | range | [11.5, 12.5] | 1
vrip | range | [0, 2] | 1
TABLE I: Parameters, Objectives, and Boundary Conditions for Boost Converter.

The CARS method found in total $2,477$ samples satisfying all boundary
conditions, while GAs found $3,367$. Valid parameters in figure 5 show for
CARS primarily one large, densely sampled parameter range. The resulting,
valid GA samples cover a wider parameter range, but samples tend to be grouped
in distinct clusters. Regarding measurements, presented in Figure 6, highest
”efficiencies” were achieved primarily by GA samples, the desired output
voltage of $12V$ was reached by both methods, and a CARS samples achieved the
lowest voltage ripple / overshoot.

Figure 5: Parameter values of samples that satisfied Boost Converter boundary
conditions (CARS blue, GA orange). Figure 6: Simulation results for samples
that satisfied Boost Converter boundary conditions (CARS blue, GA orange)

###  V-B LLC-2018

Figure 7 shows a simplified LLC Resonant Converter, highlighting the most
relevant components. The complete circuit for simulation is provided in the
appendix, figure 15. Parameters are listed in Table II. This application
considers five operating points, where hardware components are fixed, but five
different switching frequencies are applied to reach different power levels
for different output voltages. The output voltage is defined as categorical
$grid$ parameter with only one possible value for each operating point. A
total of $100,000$ samples are investigated for each method. Further details
on measurements are provided in the appendix VII-B.

Figure 7: Exemplary schematic of an LLC resonant converter, detailed
schematics including static loss models are provided in the appendix figures
15 and 16. Parameters
---
Name | Type | Values | Operating Points
C_res | log | [1e-9, 1.6e-6] | 1
L_ser | linear | [5e-6, 40e-6] | 1
L1 | log | [50e-6, 1e-3] | 1
n | linear | [0.6, 2.2] | 1
f_sw | linear | [50e3, 320e3] | 5
V_out | grid | 300 | 5
350
400
450
500
Objectives
Name | Type | Values | Operating Points
i_pri_res_rms | min | - | 5
i_sec_res_rms | min | - | 5
p_out_avg | target | 2660 | 5
3,160
3,660
3,660
3,660
Boundary Conditions
Name | Type | Values | Operating Points
psi_grad | range | [10, 90] | 5
p_out_avg | range | [1995, 3325] | 5
[2,370, 3,950]
[2,745, 4,575]
[2,745, 4,575]
[2,745, 4,575]
TABLE II: Parameters, Objectives, and Boundary Conditions for LLC-2018.

In contrast to the previous test, for LLC-2018, CARS achieved $48,501$ samples
satisfying all boundary conditions, and GAs $11,684$. Regarding valid
parameter ranges, shown in figure 8, CARS samples again tend to be grouped
into one large parameter range, while GA samples form multiple sample
clusters. Measurements, presented in figure 9, do not allow for a simple
comparison, since overall performance would require an averaged measure
considering all operating points. But the results indicate that ”locally” best
performance for one of the operating points was partially achieved by CARS
samples and partially by GA samples. It is however worth noting, that GAs
discovered additional clusters for large $L\\_ser$ and small $f\\_sw$ values
that were not considered by CARS.

Figure 8: Parameter values of samples that satisfied LLC-2018 boundary
conditions (CARS blue, GA orange) Figure 9: Simulation results for samples
that satisfied LLC-2018 boundary conditions (CARS blue, GA orange)

###  V-C LLC-2022

The third use-case, LLC-2022, shown in detail in Figure 16 in the appendix,
considers three operating points simultaneously, and a large output power
range, as listed in Table III. Again, $100,000$ samples are evaluated for each
method. The third objective $fsw$ accounts for the fact that large power
ranges shall be covered and encourages small switching frequency deviations,
which is computed according to (9) for all operating points simultaneously.

Parameters
---
Name | Type | Values | Operating Points
L_st1 | log | [2e-6, 200e-6] | 1
A_L | log | [20e-9, 500e-9] | 1
np | log | [2, 200] | 1
ns | log | [2, 200] | 1
C_res1 | log | [1e-9, 100e-9] | 1
fsw | linear | [100e3, 1e6] | 3
V_output | grid | 850 | 3
840
600
Objectives
Name | Type | Values | Operating Points
eta | max | - | 3
mag_current_rms | min | - | 3
fsw | $\text{min}_{\text{range}}$ | - | all
Boundary Conditions
Name | Type | Values | Operating Points
psi_deg | range | [5, 85] | 3
eta | range | [0.75, 1.0] | 3
output_power_avg | range | [100, 300] | 3
[3,600, 4,000]
[3,000, 3,400]
i_swprimoff | larger | 0 | 3
TABLE III: Parameters, Objectives, and Boundary Conditions for LLC-2022.

This third use case presents $719$ samples satisfying all boundary conditions
for CARS, in comparison to $14,358$ for GAs. Valid parameters, presented in
figure 10 show two smaller parameter ranges of valid samples for CARS and
multiple separate clusters for GAs. Measurements, presented in figure 11
suggest superior performance for GAs, with larger efficiency values $eta$
usually achieved by GAs, as well as lower transformer currents
$mag\\_current\\_rms$.

Figure 10: Parameter values of samples that satisfied LLC-2022 boundary
conditions (CARS blue, GA orange) Figure 11: Simulation results for samples
that satisfied LLC-2022 boundary conditions (CARS blue, GA orange)

###  V-D CARS extension study

To investigate performance benefits due to the extensions described in
sections II-A and II-B, all three use cases were investigated without
extensions for comparison. The minimum, maximum and mean fitness value for the
default setting described in section II-C with both extensions, using
oversampling but no pooling, using pooling but no oversampling and using no
extension are presented in figures 12, 13 and 13 respectively. Note that
fitness values were normalized based on the maximum and minimum values sampled
during the first iteration and may therefore exceed the previously best
samples with fitness $1.0$. Similarly, fitness values itself do not allow for
performance comparisons since normalization constants vary between different
CARS iterations.

Especially for the more complex use cases LLC-2018 and LLC-2022 the advantages
of the extensions are well observable. As intended, pooling leads to faster
increase in sampled fitness values since approximated fitness information is
distributed to adjacent sub-domains. Similarly, oversampling leads to less
variance, since randomly selected samples with poor neighbor-based fitness
estimates are omitted from simulation. Finally, even tests using no extension
suggest that softmax based sampling with increasing $\alpha$ prioritization
can suggest useful samples. Assuming representative samples during the first
batch (which is performed as random sampling, due to constant sub-domain
fitness initialization), the increase in mean fitness indicates more emphasis
on relevant parameter areas compared to random sampling.

Figure 12: Effect of extensions (pooling and oversampling) on sampled fitness
for Boost Converter. Figure 13: Effect of extensions (pooling and
oversampling) on sampled fitness for LLC-2018. Figure 14: Effect of extensions
(pooling and oversampling) on sampled fitness for LLC-2022.

##  VI Conclusion

Overall, the proposed CARS method is able to suggest valid parameter settings,
satisfying all boundary conditions for all three power electronic use cases.
Arguments regarding superiority of any method seem not justified based on the
exemplary applications. GAs might reach slightly better objective values,
while CARS suggests more diverse and continuous ranges of valid samples.
However, there are some recurring properties of both methods: CARS tends to
thoroughly investigate large parameter ranges, likely related to highly rated
and therefore frequently sampled sub-domains of the parameter tensor. In
contrast, GAs tend to find valid samples in small clusters, which is likely
related to multiple parallel GAs sampling independently on distinct ”islands”.
GAs occasionally provides samples in parameter ranges that were not
successfully investigated by CARS, but CARS indicates that valid samples exist
also between clusters found by GAs. This highlights the fact that in contrast
to globally optimal approaches, neither of both methods can find all valid
samples.

An additional advantage of the proposed CARS method is fast sampling of large
sample batches (up to a million samples) in a few seconds, which is
advantageous for very fast simulations or in cases where time consuming
simulations can be performed on multiple HPC machines in parallel. Additional
extensions allow CARS to investigate large parameter ranges faster (pooling)
and more focused (oversampling), while even softmax sampling itself seems to
improve on random sampling. Finally, sub-domain indices of CARS samples can be
stored alongside simulated samples, to re-establish previous parameter tensors
and to perform additional simulations using for example specialized
explorative or exploitative softmax weighting settings.

##  VII Discussion

As previous conclusions highlight different advantages and drawbacks for CARS
and GAs, further research might investigate a combination of both methods.
Well performing parameter settings discovered by GAs could be associated with
their respective CARS parameter tensor sub-domain. In this way the ability of
GAs to find many locally clustered solutions can be combined with the ability
of CARS to investigate large areas around well-performing samples thoroughly.

The suggested square-root Canberra distance measure has nice properties for
the considered use-cases, but it will saturate at a distance of $1$ for
negative values, which leads to no information regarding closeness.
Applications including negative target values should consider this fact and
use another distance metric.

To ease the application for electronic engineers, both optimization
strategies, CARS and GAs used experience driven heuristics for various
hyperparameters. This led to reasonable results for the considered use cases,
but it is also likely that application specific hyperparameter settings would
yield different and eventually better results. The extension study V-D
highlights an interesting situation, where especially for the more complex
use-cases (LLC-2018 and LLC-2022) there are periods of larger ”fitness slope”.
Since the prioritization parameter $\alpha$ corresponds by default to the
iteration, this increase in average fitness seems related to a combination of
the prioritization parameter $\alpha$ and the information density available in
the parameter tensor, which is additionally increased by pooling and
oversampling.

The three power electronic use-cases described in Section V were chosen to
represent a variety of different parameter design applications, but it cannot
be guaranteed that either of the investigated methods performs well on other
applications. Additionally, the considered netlists consider simplified,
static loss models or no loss models at all. Therefore, parameters suggested
by any of the methods are unlikely ”ideal” parameters for any of the use
cases, but instead rather some initial suggestion for more detailed
investigations.

Finally, the CARS method used a constant fitness value of $0.75$ as
”optimistic initialization” of the parameter sub-domain tensor. It is
therefore required to investigate all sub-domains individually (or approximate
them by max pooling or oversampling). However, the parameter sub-domain tensor
could also be initialized based on engineering knowledge about desirable and
undesirable parameter ranges, which reduces the need for evaluating irrelevant
areas through circuit simulation and might improve sample efficiency
additionally.

## Acknowledgment

This work was supported by the Federal Ministry of Education and Research in
the CODAPE ”Kollaborative Entwicklungsumgebung für die Leistungselektronik”
project via grant 16ME0356.

## References

  * [1] S. Zhao, F. Blaabjerg, H. Wang (2021). An Overview of Artificial Intelligence Applications for Power Electronics. In: IEEE Trans. Power Electron. 36 (4), S. 4633–4658. DOI: 10.1109/TPEL.2020.3024914.
  * [2] S. E. De Leon-Aldaco, H. Calleja, J. A. Alquicira, Jesus (2015). Metaheuristic Optimization Methods Applied to Power Converters: A Review. In: IEEE Trans. Power Electron. 30 (12), S. 6791–6803. DOI: 10.1109/TPEL.2015.2397311.
  * [3] A. Rosskopf, S. Volmering, S. Ditze, C. Joffe and E. Baer (2018). Autonomous circuit design of a resonant converter (LLC) for on-board chargers using genetic algorithms. In: 2018 IEEE Transportation and Electrification Conference and Expo (ITEC). S. 96–101.
  * [4] T. Guillod, J. W. Kolar (2020). From brute force grid search to artificial intelligence: Which algorithms for magnetics optimization. Ind. Session at IEEE Energy Convers. Congr. and Expo. (ECCE USA), Mar. 2020.
  * [5] R. S. Sutton, A. Barto (2020). Reinforcement learning. An introduction. Second edition. Cambridge, Massachusetts, London, England: The MIT Press (Adaptive computation and machine learning).
  * [6] C. Bishop (2009): Pattern recognition and machine learning. Corrected at 8th printing 2009. New York, NY: Springer (Information science and statistics).
  * [7] T. Schaul, J. Quan, I. Antonoglou, D. Silver (2015). Prioritized Experience Replay.
  * [8] C. R. Harris, K. J. Millman, ; S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau et al. (2020): Array programming with NumPy. In: Nature 585 (7825), S. 357–362. DOI: 10.1038/s41586-020-2649-2 .
  * [9] R. Chelouah, P. Siarry (2000). Tabu Search applied to global optimization. In: European Journal of Operational Research 123 (2), S. 256–270. DOI: 10.1016/S0377-2217(99)00255-6 .
  * [10] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, vol. 6, no. 2, pp. 182-197. DOI: 10.1109/4235.996017 .
  * [11] D. Scherer, A. Müller, S. Behnke (2010). Evaluation of Pooling Operations in Convolutional Architectures for Object Recognition. In: David Hutchison, Takeo Kanade, Josef Kittler, Jon M. Kleinberg, Friedemann Mattern, John C. Mitchell et al. (Hg.): Artificial Neural Networks – ICANN 2010, Bd. 6354. Berlin, Heidelberg: Springer Berlin Heidelberg (Lecture Notes in Computer Science), S. 92–101.
  * [12] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel et al. (2011). Scikit-learn: Machine Learning in Python. In: Journal of Machine Learning Research 12 (85), S. 2825–2830. http://jmlr.org/papers/v12/pedregosa11a.html.
  * [13] T. Fühner (2007): Direct optimization approach for lithographic process conditions. In: J. Micro/Nanolith. MEMS MOEMS 6 (3), S. 31006. DOI: 10.1117/1.2785031 .
  * [14] K. F. Man, K. S. Tang, S. Kwong (1996): Genetic algorithms: concepts and applications [in engineering design]. In: IEEE Trans. Ind. Electron. 43 (5), S. 519–534. DOI: 10.1109/41.538609 .
  * [15] J. M. Johnson, V. Rahmat-Samii (1997): Genetic algorithms in engineering electromagnetics. In: IEEE Antennas Propag. Mag. 39 (4), S. 7–21. DOI: 10.1109/74.632992 .
  * [16] N. Srinivas, K. Deb (1994). Muiltiobjective Optimization Using Nondominated Sorting in Genetic Algorithms. Evolutionary Computation, vol. 2, no. 3, pp. 221-248, Sept. 1994. DOI: 10.1162/evco.1994.2.3.221.
  * [17] J. Schaffer (1985). Multiple Objective Optimization with Vector Evaluated Genetic Algorithms. Proceedings of the First Int. Conference on Genetic Algortihms, Ed. G.J.E Grefensette, J.J. Lawrence Erlbraum. 93-100.
  * [18] A. Rezaee Jordehi (2015): A review on constraint handling strategies in particle swarm optimisation. In: Neural Comput & Applic 26 (6), S. 1265–1275. DOI: 10.1007/s00521-014-1808-5 .
  * [19] Fraunhofer Institute for Integrated Systems and Device Technology IISB, SchOKI – AI based circuit optimization for dc/dc converters, https://schoki.app/, LZE GmbH, 2023.
  * [20] Analog Devices, LTspice, https://www.analog.com/en/design-center/design-tools-and-calculators/ltspice-simulator.html, Analog Devices, 2023.
  * [21] K. Deb, K. Sindhya, and T. Okabe (2007). Self-adaptive simulated binary crossover for real-parameter optimization. In Proceedings of the 9th annual conference on Genetic and evolutionary computation (GECCO ’07). Association for Computing Machinery, 1187–1194. DOI: 10.1145/1276958.1277190 .

## Appendix

###  VII-A GA settings

Important parameters for genetic algorithms are the chosen number of islands,
their population sizes, as well as the number of generations across which the
populations are evolved. Table IV lists the corresponding choices for the
three considered use cases in this paper.

GA Settings
---
Use-Case | # Islands | Population size | Generations
Boost Converter | 5 | 20 | 50
LLC-2018 | 50 | 20 | 100
LLC-2022 | 50 | 20 | 100
TABLE IV: GA setups for the three power electronic use cases.

Furthermore, the choice and probabilities for mutations and crossovers are
decisive for a good trade-off between exploration of new areas in the
parameter space and exploitation in the close surroundings of already
evaluated individuals. To allow partial evolution, mutation happens based on
two probabilities $P_{mutate}$, which determines whether an individual is
mutated at all and $P_{mutate\\_val}$, which determines the likelihood of each
individual value being mutated in case the individual is mutated at all.
Mutation is based on a Gaussian distribution with standard deviation
$\sigma_{mutate}$. Furthermore, we use simulated binary crossover with
hyperparamaters $\alpha_{crossover}$ and $\eta_{crossover}$ determining the
general resemblance of the product with respect to the two individuals used
for crossover. We refer to [21] for more information on this crossover
technique. Table V lists the parameters for mutation and crossover, which are
promoting diversity within the populations and have proven superior for the
considered use cases over less explorative configurations.

Parameter | Value
---|---
$P_{\text{mutate}}$ | 0.1
$P_{\text{mutate\\_val}}$ | 0.3
$\sigma_{\text{mutate}}$ | 0.4
$P_{\text{crossover}}$ | 0.95
$\alpha_{\text{crossover}}$ | 0.2
$\eta_{\text{crossover}}$ | 1
TABLE V: Hyperparameters used for Mutation and Crossover

###  VII-B Circuit details

For the textbook Boost Converter, presented in Section V-A, the mean output
voltage $vmean$ is obtained as mean voltage on the output load resistor after
a startup phase of $10ms$. The ripple / overshoot measure $vrip$ subtracts the
mean output voltage $vmean$ from the maximum output voltage, also measured
after a startup phase of $10ms$. The simplified efficiency measure $eff\\_tot$
divides total output power by the total input power. Since the schematic does
not include component loss models, this efficiency measure primarily considers
stored energy, e.g., in the output capacitor. This way there are competing
objectives, where the $vrip$ measure requires large capacitors for smoothing,
while the efficiency measure $eff\\_tot$ suggests smaller capacitors to
achieve large output powers.

Figure 15 provides the schematic used for LLC-2018 simulations. The
measurements $i\\_pri\\_res\\_rms$ and $i\\_sec\\_res\\_rms$ are RMS current
measurements on the respective components. Similarly output power
$p\\_out\\_avg$ is the current and voltage product in the output sink. The
phase shift measure $psi\\_grad$ is used to determine whether zero voltage
switching is possible. Therefore, the primary side zero-crossing is measured
as $psi$ and $psi\\_grad$ computed as:

$psi\\_grad=360\cdot f\\_sw\cdot psi.$ (11)

Figure 15: Schematic used for LLC-2018 simulations.

Finally, figure 16 provides the schematic used for LLC-2022 simulations. The
output power $output\\_power\\_avg$ is again computed as product of current
and voltage in the output sink. The efficiency measure $eta$ is computed based
on the output power $output\\_power\\_avg$ and static component losses.

$eta=\frac{output\\_power\\_avg}{output\\_power\\_avg+losses}.$ (12)

The current $i\\_swprimoff$ is measured in the respective component. The phase
shift measure $psi\\_deg$ is obtained by measuring the zero crossing of the
primary voltage as $psi$ and computed as:

$psi\\_deg=360\cdot fsw\cdot psi.$ (13)

The magnetic current $mag\\_current\\_rms$ is measured as current in the
primary stray inductance subtracted by the turns-ratio factored current of the
secondary stray inductance.

$mag\\_current\\_rms=I(Lst1)-\frac{ns}{np}\cdot I(Lst2).$ (14)

Figure 16: Schematic used for LLC-2022 simulations.


# Hunting for tetraquarks in ultra-pheripheral heavy ion collisions

Angelo Esposito<EMAIL_ADDRESS>Theoretical Particle Physics
Laboratory (LPTP), Institute of Physics, EPFL, 1015 Lausanne, Switzerland
School of Natural Sciences, Institute for Advanced Study, Princeton, New
Jersey 08540, USA  Claudio Andrea Manzari<EMAIL_ADDRESS>Physik-Institut, Universität Zürich, Winterthurerstrasse 190, CH-8057 Zürich,
Switzerland  Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland
Alessandro Pilloni  INFN Sezione di Roma, I-00185, Rome, Italy  Dipartimento
di Scienze Matematiche e Informatiche, Scienze Fisiche e Scienze della Terra,
Università degli Studi di Messina, I-98122 Messina, Italy  INFN Sezione di
Catania, I-95123 Catania, Italy  Antonio Davide Polosa  INFN Sezione di Roma,
I-00185, Rome, Italy  Dipartimento di Fisica, Sapienza Università di Roma,
I-00185 Roma, Italy

###### Abstract

Ultra-peripheral heavy ion collisions constitute an ideal setup to look for
exotic hadrons because of their low event multiplicity and the possibility of
an efficient background rejection. We propose to look for four-quark states
produced by photon-photon fusion in these collisions at the center-of-mass
energy per nucleon pair $\sqrt{s_{\text{NN}}}=5.5~{}{\rm Te\kern-1.00006ptV}$.
In particular, we focus on those states that would represent a definite
smoking gun for the compact tetraquark model. We show that the $X(6900)$, a
likely $cc\bar{c}\bar{c}$ compact state, is a perfect candidate for this
search, and estimate a production cross section ranging from around $250$ nb
to $1150$ nb, depending on its quantum numbers. Furthermore, we discuss the
importance of ultra-peripheral collisions to the search for the scalar and
tensor partners of the $X(3872)$ predicted by the diquarkonium model, and not
yet observed. The completion of such a flavor-spin multiplet would speak
strongly in favor of the compact tetraquark model.

††preprint: ZU-TH 40/21, PSI-PR-21-20

##  I Introduction

The existence of hadrons with more than three valence constituents is now well
assessed Esposito  _et al._ (2017); Guo  _et al._ (2018); Olsen  _et al._
(2018); Brambilla  _et al._ (2020), but the understanding of their nature
remains a long standing problem of low-energy QCD. Are these states extended
hadronic molecules arising from color neutral interactions? Or are they rather
compact tetraquarks generated by short distance forces, analogs to mesons and
baryons?

The solution to this issue requires the identification of some smoking guns,
able to clearly discriminate between the two models. One such possibility is
the recent observation by LHCb of a narrow resonance in the di-$J/\psi$ mass
spectrum Aaij  _et al._ (2020a), dubbed $X(6900)$ and compatible with a
$cc\bar{c}\bar{c}$ structure. The possibility of such a state was already
anticipated by several studies, and later further investigated (see, e.g.,
Chao (1981); Heller and Tjon (1985); Barnea  _et al._ (2006); Vijande  _et
al._ (2007); Ebert  _et al._ (2007); Berezhnoy  _et al._ (2012); Wu  _et al._
(2018); Chen  _et al._ (2017); Wang (2017); Debastiani and Navarra (2019);
Richard  _et al._ (2017); Anwar  _et al._ (2018); Bedolla  _et al._ (2020);
Karliner  _et al._ (2017); Becchi  _et al._ (2020); Dong  _et al._ (2021a);
Cao  _et al._ (2021); Liang  _et al._ (2021); Wang  _et al._ (2021); Dong  _et
al._ (2021b)). Crucially, no single light hadron can mediate the interaction
between charmonia to generate a loosely bound molecule Maiani (2020); Dong
_et al._ (2021b). The $X(6900)$ seems likely to be a compact tetraquark.

Another compelling indication of the tetraquark nature of the exotic states
would be the observation of a complete flavor-spin multiplet, as predicted in
Maiani  _et al._ (2014). In the hidden charm sector, the $J^{PC}=1^{+-}$
resonances—the so-called $Z_{c}(3900)$ and $Z_{c}^{\prime}(4020)$—have been
observed in three charge states, while the $1^{++}$ one—the famous
$X(3872)$—has only been observed in a single neutral component. Besides the
charged partners of the $X(3872)$, to complete the $c\bar{c}q\bar{q}$
multiplet, one would have to observe the predicted scalar and tensor states
Maiani  _et al._ (2014).

In this work we propose to look for the above-mentioned smoking guns in ultra-
peripheral heavy ion collisions (UPCs) at the LHC. In these events the impact
parameter is much larger than the ions’ radii, which then scatter off each
other elastically Baur  _et al._ (2002); Bertulani  _et al._ (2005); Baltz
(2008). This causes a lack of additional calorimetric signals and a large
rapidity gap between the particles produced and the outgoing beams, which can
be used for an efficient background rejection. For this reason they are an
optimal environment for exotic searches, ranging from hadronic states to extra
dimensions (see, e.g., Grabiak  _et al._ (1989); Drees  _et al._ (1989);
Greiner  _et al._ (1993); Ahern  _et al._ (2000); Bertulani (2009); Goncalves
_et al._ (2013); Moreira  _et al._ (2016); Knapen  _et al._ (2017); Goncalves
and Moreira (2019); Gonçalves and Moreira (2021)). These collisions are
particularly amenable to search for states, like the ones of interest to us,
that can be produced by photon-photon fusion. Indeed, the large charge of lead
ions ($Z=82$) induces a huge $Z^{4}$ enhancement in the coherent photon–photon
luminosity, consequently boosting the production cross section for these
states.

The results we find are very encouraging. At the center-of-mass energy per
nucleon pair $\sqrt{s_{\text{NN}}}=5.5~{}{\rm Te\kern-1.00006ptV}$, both the
$X(6900)$ and the scalar and tensor $c\bar{c}q\bar{q}$ states are expected to
be copiously produced in UPCs. In particular, due to its likely large width
into vector charmonia, the $X(6900)$ should be produced with cross sections of
the order of fractions of microbarn, or even more. The scalar and tensor
states of the $X(3872)$ multiplet should instead be produced with cross
sections larger than the measured one of the $X(3872)$ in prompt $pp$
collisions Chatrchyan  _et al._ (2013); Aaij  _et al._ (2021a). The
observation of these states in UPCs would be another indication of the
existence of compact tetraquarks in the spectrum of short distance QCD,
alongside with a recently emerging pattern which includes the observation of
the hidden charm and strange states Ablikim  _et al._ (2021); Aaij  _et al._
(2021b); Maiani  _et al._ (2021) and the study of the lineshape of the
$X(3872)$ Aaij  _et al._ (2020b); Esposito  _et al._ (2021).

##  II Photon–photon interaction

When two ions pass each other at distances larger than their radii they
interact solely via their electromagnetic fields. For relativistic ions with
$Z\gg 1$, the electric and magnetic fields are perpendicular, and the
configuration may be represented as a flux of almost-real photons following
the Wizsäcker–Williams method von Weizsacker (1934); Williams (1934). In
particular, the number of photons per unit area and energy emitted by an ion
with boost factor $\gamma\gg 1$ is given by Klein  _et al._ (2017)

$\displaystyle\frac{dN_{\gamma}(k,\bm{b})}{dk\,d^{2}b}=\frac{Z^{2}\alpha}{\pi^{2}}\frac{k}{\gamma^{2}}K_{1}^{2}\left(\frac{k|\bm{b}|}{\gamma}\right)\,,$
(1)

where $k$ is the photon energy, $\bm{b}$ is the transverse distance from the
moving ion, $\alpha$ is the electromagnetic fine structure constant and
$K_{1}$ is the modified Bessel functions. Since the photons are quasireal, in
Eq. (1) only the flux of transversely polarized photons has been considered.

In an UPC the two-photon luminosity is given by

$\displaystyle\begin{split}\frac{d^{2}N_{\gamma\gamma}(k_{1},k_{2})}{dk_{1}dk_{2}}=&\int
d^{2}b_{1}d^{2}b_{2}\,P_{\rm NOHAD}\big{(}|\bm{b}_{1}-\bm{b}_{2}|\big{)}\\\
&\times\frac{dN_{\gamma}(k_{1},\bm{b}_{1})}{dk_{1}\,d^{2}b_{1}}\frac{dN_{\gamma}(k_{2},\bm{b}_{2})}{dk_{2}\,d^{2}b_{2}}\,,\end{split}$
(2)

which evidently features a $Z^{4}$ enhancement—see Eq. (1). The requirement
that the two nuclei do not interact hadronically is imposed by $P_{\rm
NOHAD}(b)$, which is the probability of having no hadronic interactions at
impact parameter $b$. In what follows we use the STARlight code Klein  _et
al._ (2017), where

$\displaystyle P_{\rm NOHAD}(b)=\mathrm{e}^{-\sigma_{\rm NN}T_{\rm AA}(b)}\,,$
(3)

with $\sigma_{\rm NN}$ the nucleon–nucleon interaction cross section, and
$T_{\rm AA}(b)=\int d^{2}b_{1}\,T_{A}(b_{1})\,T_{A}(|\bm{b}_{1}-\bm{b}|)$ the
nuclear overlap function determined from the Woods–Saxon nuclear density
distributions of the two nuclei, $T_{A}(b)$.

We are interested in processes where the two photons produce a state, $X$,
with invariant mass $W=\sqrt{4k_{1}k_{2}}$ and rapidity
$Y=\frac{1}{2}\ln(k_{1}/k_{2})$. The cross section for such a process
factorizes in two terms: the two-photon luminosity associated to the incoming
nuclei, and the cross section, $\sigma_{\gamma\gamma}(W)$, for the creation of
$X$ from two photons, i.e.

$\displaystyle\sigma(\ce{Pb}\,\ce{Pb}\to\ce{Pb}\,\ce{Pb}\,X)=\int
dYdW\frac{d^{2}N_{\gamma\gamma}}{dWdY}\,\sigma_{\gamma\gamma}(W)\,$

The cross section to produce a single meson in a photon–photon interaction is
given by Klein  _et al._ (2017)

$\displaystyle\begin{split}\sigma_{\gamma\gamma}(W)&=8\pi(2J+1)\frac{\Gamma_{\gamma\gamma}\Gamma}{(W^{2}-m^{2})^{2}+m^{2}\Gamma^{2}}\\\
&\simeq
8\pi^{2}(2J+1)\frac{\Gamma_{\gamma\gamma}}{2m^{2}}\delta(W-m)\,,\end{split}$
(4)

where $m$ is the meson mass, $\Gamma_{\gamma\gamma}$ is its width in two
photons, $\Gamma$ is the total width, $J$ is its spin. The last step realizes
the narrow width approximation. From here we see that lighter and higher spin
particles are produced more copiously.

###  II.1 Partial widths into $\gamma\gamma$

It is clear that the central quantity in this formalism is the partial width
of the state $X$ in two photons. In what follows we will compute it using the
vector meson dominance model Sakurai (1960). In this picture, the radiative
decay of a hadron happens first via its decay into vector mesons, which then
mix with photons—see Figure 1.

Figure 1: Decay of the state $X$ into $\gamma\gamma$ via two vectors. | $\rho$ | $\omega$ | $J/\psi$
---|---|---|---
$\Gamma_{ee}$ (${\rm Ge\kern-1.00006ptV}$) |  $\simeq 7.0\times 10^{-6}$ |  $\simeq 6.4\times 10^{-7}$ |  $\simeq 5.5\times 10^{-6}$
$f_{V}$ (${\rm Ge\kern-1.00006ptV}^{2}$) |  $\simeq 0.16$ |  $\simeq 0.17$ |  $\simeq 1.3$
Table 1: Electronic widths as taken from PDG Zyla  _et al._ (2020) and the
decay constants extracted from them. The results are consistent with what was
found in Casalbuoni  _et al._ (1997); Deandrea  _et al._ (2002).

In particular, the vector–photon mixing is given by

$\displaystyle=\kappa_{V}\,e\,f_{V}\,\eta^{\mu\nu}\,,$

with
$\kappa_{V}=\big{(}\frac{1}{\sqrt{2}},\frac{1}{3\sqrt{2}},\frac{2}{3}\big{)}$
if $V=(\rho,\omega,J/\psi)$.111 The vector-photon mixing is obtained from the
standard electromagnetic Lagrangian,
$\mathcal{L}=A_{\mu}\sum_{a}Q_{a}e\bar{q}_{a}\gamma^{\mu}q_{a}$, together with
the meson states
$|\rho\rangle=\frac{|u\bar{u}\rangle-|d\bar{d}\rangle}{\sqrt{2}}$,
$|\omega\rangle=\frac{|u\bar{u}\rangle+|d\bar{d}\rangle}{\sqrt{2}}$ and
$|J/\psi\rangle=|c\bar{c}\rangle$. The decay constants are defined through the
matrix element $\langle
0|\sum_{a}Q_{a}e\bar{q}_{a}\gamma^{\mu}q_{a}|V\rangle=\kappa_{V}f_{V}\epsilon^{\mu}$.
The decay constants $f_{V}$ can instead be extracted from the electronic width
of the corresponding vector, $\Gamma(V\to
e^{+}e^{-})=4\pi\alpha^{2}\kappa_{V}^{2}f_{V}^{2}/(3m_{V}^{3})$. In Table 1 we
report the electronic widths and the corresponding mixing constants.

The most general matrix elements for the decay of the scalar and tensor exotic
mesons in two vectors can be written as

$\displaystyle\langle X(0^{++})|V_{1}V_{2}\rangle$
$\displaystyle=\alpha_{0}\,\epsilon_{1}\cdot\epsilon_{2}+\beta_{0}\left(\epsilon_{1}\cdot
k_{2}\right)\left(\epsilon_{2}\cdot k_{1}\right)\,,$ (8a)
$\displaystyle\langle X(2^{++})|V_{1}V_{2}\rangle$
$\displaystyle=\pi_{\mu\nu}\big{[}\alpha_{2}\epsilon_{1}^{\mu}\epsilon_{2}^{\nu}+\beta_{2}\left(\epsilon_{1}\cdot
k_{2}\right)\epsilon_{2}^{\mu}k_{1}^{\nu}$ (8b)
$\displaystyle+\beta_{2}^{\prime}\left(\epsilon_{2}\cdot
k_{1}\right)\epsilon_{1}^{\mu}k_{2}^{\nu}+\gamma_{2}\left(\epsilon_{1}\cdot\epsilon_{2}\right)k_{1}^{\mu}k_{2}^{\nu}$
$\displaystyle+\delta_{2}\left(\epsilon_{1}\cdot
k_{2}\right)\left(\epsilon_{2}\cdot
k_{1}\right)k_{1}^{\mu}k_{2}^{\nu}\big{]}\,,$

where $\epsilon_{i}$ and $k_{i}$ are the polarization and momentum of the
vector $V_{i}$, and $\pi_{\mu\nu}$ is the polarization of the tensor.222The
sum over spin-2 polarizations is given by Faccini  _et al._ (2012); Gleisberg
_et al._ (2003)
$\sum_{\text{pol}}\pi_{\mu\nu}(k)\pi_{\rho\sigma}(k)=\frac{1}{2}\mathcal{P}_{\mu\rho}\mathcal{P}_{\nu\sigma}+\frac{1}{2}\mathcal{P}_{\mu\sigma}\mathcal{P}_{\nu\rho}-\frac{1}{3}\mathcal{P}_{\mu\nu}\mathcal{P}_{\rho\sigma}$,
with $\mathcal{P}_{\mu\nu}=-\eta_{\mu\nu}+k_{\mu}k_{\nu}/m^{2}$ and
$k^{2}=m^{2}$. In absence of further information it is impossible to determine
all the above couplings from the data. We will therefore adopt a minimal
model, somewhat inspired from an EFT approach, and neglect all terms
proportional to the particle momenta (see, e.g., Mathieu  _et al._
(2020)).333For the scalar case, we checked that including the $\beta_{0}$
coefficient and letting it vary around its natural value, $\beta_{0}\sim
1/m_{X}$, does not change the order of magnitude estimates of Table 2.

###  II.2 Production of the $X(6900)$

As already mentioned, the $X(6900)$ is, in all likelihood, a compact
$cc\bar{c}\bar{c}$ state. Its mass and width are $m_{X}=6886\pm 2~{}{\rm
Me\kern-1.00006ptV}$ and $\Gamma_{X}=168\pm 102~{}{\rm Me\kern-1.00006ptV}$
Aaij  _et al._ (2020a), while its quantum numbers are yet to be determined.
Were it to have $J^{PC}=0^{++}$ or $2^{++}$, it could be produced from
photon–photon fusion in UPCs, as also discussed in Gonçalves and Moreira
(2021).

To provide an order of magnitude estimate of its partial width in two photons
we make the assumption that its coupling to vector mesons is dominated by the
di-$J/\psi$ one Karliner  _et al._ (2017). Indeed, with four heavy quarks
involved, the coupling to light vector mesons involves annihilation processes,
and are thus OZI-suppressed by powers of $\alpha_{s}(4m_{c})$. The
contribution to the vector meson dominance from excited charmonia is also
suppressed by their greater spatial extent Redlich  _et al._ (2000).444In
Gonçalves and Moreira (2021) the partial width of the $X(6900)$ in two photons
is taken to be the same as the $\chi_{cJ}$ quarkonium with the same quantum
numbers. This underlines the somewhat strong assumption that the short
distance dynamics of the two states is the same, which is not guaranteed. Here
we take a more conservative approach and keep the branching ratio unspecified.

Starting from the matrix elements in Eqs. (8), and using the vector meson
dominance as in Figure 1, we can obtain the partial width in two photons. Note
that, since the amplitude is not gauge invariant, one must restrict oneself to
the transverse photon polarizations. The results for the scalar and tensor
case are

$\displaystyle\Gamma_{\gamma\gamma}^{0^{++}}$
$\displaystyle=\frac{16\pi\alpha^{2}}{81}\frac{\alpha_{0}^{2}}{m_{X}}\left(\frac{f_{\psi}}{m^{2}_{\psi}}\right)^{4}\,,$
(9a) $\displaystyle\Gamma_{\gamma\gamma}^{2^{++}}$
$\displaystyle=\frac{56\pi\alpha^{2}}{1215}\frac{\alpha_{2}^{2}}{m_{X}}\left(\frac{f_{\psi}}{m_{\psi}^{2}}\right)^{4}\,.$
(9b)

The couplings $\alpha_{J}$ can be extracted from the partial width of the
$X(6900)$ in di-$J/\psi$. Since the corresponding branching ratio is yet
unknown, we will keep it general, bearing in mind that it is likely that this
channel will dominate the total width Karliner  _et al._ (2017). For the
scalar and tensor cases one gets, respectively

$\displaystyle\mathcal{B}_{\psi}\,\Gamma_{X}$
$\displaystyle=\frac{\alpha_{0}^{2}\,p}{16\pi
m^{2}_{X}}\left(3-\frac{m_{X}^{2}}{m_{\psi}^{2}}+\frac{1}{4}\frac{m_{X}^{4}}{m_{\psi}^{4}}\right)\,,$
(10a) and $\displaystyle\mathcal{B}_{\psi}\,\Gamma_{X}$
$\displaystyle=\frac{\alpha_{2}^{2}\,p}{16\pi
m_{X}^{2}}\left(\frac{7}{15}+\frac{1}{10}\frac{m_{X}^{2}}{m_{\psi}^{2}}+\frac{1}{120}\frac{m_{X}^{4}}{m_{\psi}^{4}}\right)\,,$
(10b)

where $\mathcal{B}_{\psi}$ is the branching ratio of the di-$J/\psi$ final
state, and $p=\lambda^{1/2}(m_{X}^{2},m_{\psi}^{2},m_{\psi}^{2})/(2m_{X})$ the
decay momentum, with $\lambda$ the Källén function.

In Table 2 we report the partial widths in two photons and the corresponding
cross sections for production in UPCs as obtained from the STARlight code
Klein  _et al._ (2017).

State, $J^{PC}$ |  $\Gamma_{\gamma\gamma}/\mathcal{B}_{\psi}$ (eV) |  $\sigma(\text{Pb}\text{Pb}\to\text{Pb}\text{Pb}X)/\mathcal{B}_{\psi}$ (nb)
---|---|---
$X(6900)$, $0^{++}$ | $\sim 104$ | $\sim 282$
$X(6900)$, $2^{++}$ | $\sim 86$ | $\sim 1165$
Table 2: Partial widths in two photons and corresponding production cross
sections in UPCs for the $X(6900)$, obtained for
$\sqrt{s_{\text{NN}}}=5.5~{}{\rm Te\kern-1.00006ptV}$ and normalized by the
di-$J/\psi$ braching ratio. The latter is unknown, but expected to be close to
one.

In Figure 2 we report the momentum distributions of the two $J/\psi$’s
produced by the decay of the $X(6900)$. We apply the pseudorapidity cuts
corresponding to the LHCb and ALICE acceptances. As one can see, both
experiments should be sensitive to energetic final states.

Figure 2: Momentum distributions of the two $J/\psi$’s produced by the decay
of the $X(6900)$, for the pseudorapidity cuts corresponding to the LHCb
detector (left), the ALICE electron detector (center) and the ALICE muon
detector (right). While for the ALICE electron detector the two decay products
are rather soft, the forward detectors should collect energetic ones. The
pseudorapidity cuts have efficiencies of 21%, 8%, and 10%, respectively.

###  II.3 Production of scalar and tensor $c\bar{c}q\bar{q}$ states

Contrary to the $X(6900)$, the scalar and tensor states of the
$c\bar{c}q\bar{q}$ multiplet are yet to be observed. The diquarkonium model
predicts two $J^{PC}=0^{++}$ states, dubbed $X_{0}$ and $X_{0}^{\prime}$ and
with masses around $m_{0}\simeq 3770~{}{\rm Me\kern-1.00006ptV}$ and
$m_{0^{\prime}}\simeq 4000~{}{\rm Me\kern-1.00006ptV}$, respectively Maiani
_et al._ (2014). It also predicts one $2^{++}$ state, dubbed $X_{2}$, and is
degenerate with the $X_{0}^{\prime}$, $m_{2}\simeq 4000~{}{\rm
Me\kern-1.00006ptV}$. They all have the right quantum numbers to be produced
via photon–photon fusion in UPCs.

As for the $X(3872)$, the states above can in principle decay into both
$J/\psi\,\rho$ and $J/\psi\,\omega$. Indeed, the isospin breaking mechanism
for tetraquarks holds regardless of the mass splitting in the multiplet Rossi
and Veneziano (2004). We therefore assume that the scalars and tensor share
the same isospin breaking pattern as the $X(3872)$. In terms of the spins of
the $c\bar{c}$ and $q\bar{q}$ pairs, one has Maiani  _et al._ (2014)

$\displaystyle|X_{0}\rangle=\frac{1}{2}|0_{c\bar{c}},0_{q\bar{q}}\rangle_{0}+\frac{\sqrt{3}}{2}|1_{c\bar{c}},1_{q\bar{q}}\rangle_{0}\,,$
(14a)
$\displaystyle|X_{0}^{\prime}\rangle=\frac{\sqrt{3}}{2}|0_{c\bar{c}},0_{q\bar{q}}\rangle_{0}-\frac{1}{2}|1_{c\bar{c}},1_{q\bar{q}}\rangle_{0}\,,$
(14b)
$\displaystyle|X_{1}\rangle\equiv|X(3872)\rangle=|1_{c\bar{c}},1_{q\bar{q}}\rangle_{1}\,,$
(14c) $\displaystyle|X_{2}\rangle=|1_{c\bar{c}},1_{q\bar{q}}\rangle_{2}\,,$
(14d)

where $|J_{c\bar{c}},J_{q\bar{q}}\rangle_{J}$ is a state with spin
$J_{q_{1}\bar{q}_{2}}$ for the quark-antiquark pairs and total spin $J$.
Considering again the matrix elements in Eqs. (8) under the minimal
model,555In the case of the $X(3872)$, it has been shown Maiani  _et al._
(2005) that the minimal model is able to reproduce the observed partial widths
in $J/\psi\rho$ and $J/\psi\omega$. The dynamics of the members of the same
multiplet will likely be similar, hence justifying restricting oneself to the
minimal model for the scalar and tensor states as well. one can relate their
couplings to that of the $X(3872)$. The exact relations between the couplings
to $J/\psi V$ ($V=\rho,\omega$) of the $X_{0}$, $X_{0}^{\prime}$, $X_{2}$, and
$X(3872)$ depend on the dynamics of the multiplet. Being this level of
precision negligible for the scope of this work, and referring to them
respectively as $\alpha_{0,V}$, $\alpha_{0^{\prime},V}$, $\alpha_{2,V}$ and
$\alpha_{1,V}$, we expect

$\displaystyle\frac{\alpha_{1,V}}{m_{1}}\sim\frac{\alpha_{0,V}}{m_{0}}\sim\frac{\alpha_{0^{\prime},V}}{m_{0^{\prime}}}\sim\frac{\alpha_{2,V}}{m_{2}}\,,$
(15)

where $m_{1}$ is the mass of the $X(3872)$. The matrix element for the
$X(3872)\to J/\psi V$ decay can be written as Maiani  _et al._ (2005)

$\displaystyle\langle X(3872)|J/\psi
V\rangle=\alpha_{1,V}\big{(}\bm{\epsilon}_{X}\times\bm{\epsilon}_{J/\psi}\big{)}\cdot\bm{\epsilon}_{V}\,.$
(16)

The couplings $\alpha_{1,\rho}$ and $\alpha_{1,\omega}$ can be extracted from
the branching ratios $\mathcal{B}(X(3872)\to J/\psi\pi\pi)\simeq 3.8\%$ and
$\mathcal{B}(X(3872)\to J/\psi\pi\pi\pi)\simeq 4.3\%$ Zyla  _et al._ (2020);
Brazzi  _et al._ (2011); *Faccini:2012zv; *Albaladejo:2020tzt. In particular,
the partial widths for these decays can be computed as

$\displaystyle\begin{split}&\Gamma(J/\psi f)=\mathcal{B}(V\to f)\\\
&\quad\times\int_{s_{\text{min}}}^{s_{\text{max}}}ds\,\frac{\Gamma_{V\to
f}(s)}{\Gamma_{V\to f}(m^{2}_{V})}\,\text{BW}(s)\Gamma(J/\psi
V)\,.\end{split}$ (17)

Here $f=\pi\pi$ or $\pi\pi\pi$, $\mathcal{B}(V\to f)$ is the branching ratio
for the decay of the light vector into the final state $f$, and $\Gamma(J/\psi
V)$ is the decay width of the $X(3872)$ into $J/\psi$ and a light vector of
invariant mass $s$, as computed from Eq. (16). Moreover,
$s_{\text{max}}=(m_{X}-m_{J/\psi})^{2}$, $s_{\text{min}}=(2m_{\pi})^{2}$ if
$V=\rho$ and $(3m_{\pi})^{2}$ if $V=\omega$, and $\Gamma_{V\to f}$ are the
decay rates reported in Appendix A. Using the Breit–Wigner width of the
$X(3872)$ as recently measured by LHCb, $\Gamma_{X}=1.39\pm 0.34~{}{\rm
Me\kern-1.00006ptV}$ Aaij  _et al._ (2020b), one finds

$\displaystyle\alpha_{1,\rho}\simeq 342~{}{\rm
Me\kern-1.00006ptV}\,,\quad\text{ and }\quad\alpha_{1,\omega}\simeq
1119~{}{\rm Me\kern-1.00006ptV}\,.$ (18)

Diagrams with an intermediate $\rho$ and $\omega$ will now both contribute
coherently to the total width in two photons. Starting again from Eqs. (8),
one finds, for the scalars and tensor mesons,

$\displaystyle\Gamma_{\gamma\gamma}^{X_{0}^{(\prime)}}$
$\displaystyle=\frac{\pi\alpha^{2}\kappa_{\psi}^{2}f_{\psi}^{2}}{m_{0^{(\prime)}}m_{\psi}^{4}}\left|\sum_{V=\omega,\rho}\frac{\kappa_{V}f_{V}\alpha_{0,V}}{m_{V}^{2}-im_{V}\Gamma_{V}}\right|^{2}\,,$
(19a) $\displaystyle\Gamma_{\gamma\gamma}^{X_{2}}$
$\displaystyle=\frac{7\pi\alpha^{2}\kappa_{\psi}^{2}f_{\psi}^{2}}{30m_{2}m_{\psi}^{4}}\left|\sum_{V=\omega,\rho}\frac{\kappa_{V}f_{V}\alpha_{2,V}}{m_{V}^{2}-im_{V}\Gamma_{V}}\right|^{2}\,,$
(19b)

with $\Gamma_{V}$ the width of the light vector. Putting everything found so
far together, one obtains the partial widths and production cross sections in
UPCs reported in Table 3.

State, $J^{PC}$ |  $\Gamma_{\gamma\gamma}$ (eV) |  $\sigma(\text{Pb}\,\text{Pb}\to\text{Pb}\,\text{Pb}\,X)$ (nb)
---|---|---
$X_{0}(\sim 3770)$, $0^{++}$ | $\sim 6.3$ | $\sim 185$
$X_{0}^{\prime}(\sim 4000)$, $0^{++}$ | $\sim 6.7$ | $\sim 156$
$X_{2}(\sim 4000)$, $2^{++}$ | $\sim 1.6$ | $\sim 187$
Table 3: Partial widths in two photons and corresponding production cross
sections in UPCs for the scalar and tensor elements of the $X(3872)$
multiplet. Again, the values are computed for $\sqrt{s_{\text{NN}}}=5.5~{}{\rm
Te\kern-1.00006ptV}$.

Due to the small widths of $X(3872)\to J/\psi V$, the resulting production
cross sections are smaller than that of the $X(6900)$. Nonetheless, they are
still larger than that of the $X(3872)$ as observed produced promptly in $pp$
collisions Artoisenet and Braaten (2010). Moreover, the decay of the exotic in
its final state is dominated by the $S$-wave component just like for the
$X(6900)$. For this reason, we expect similar distributions as in Figure 2.

##  III Conclusion

We proposed to look for compact tetraquarks in ultra-peripheral heavy ion
collisions. In particular, we focus on those resonances whose observation
represents a clear indication of a compact tetraquark nature.

The first is the $X(6900)$, recently discovered by LHCb and having a
$cc\bar{c}\bar{c}$ valence structure. Since there is no known mechanism that
can bind together two charmonia in a loosely bound molecule, this state is
likely compact. We find that, due to its strong coupling to a di-$J/\psi$
final state, this resonance is expected to be produced copiously in ultra-
peripheral collisions. Its study in this context would allow to shed further
light into its properties.

The other possible direction that would demonstrate the existence of four-
quark objects in short distance QCD is the observation of a complete flavor-
spin multiplet, very much analogously to what happened for standard mesons and
baryons. In particular, beside the famous charged partners of the $X(3872)$,
the missing pieces of the $S$-wave diquarkonia are the scalar and tensor
states. These too are expected to be produced in ultra-peripheral collisions,
with cross sections larger than the (large) prompt production cross section of
the $X(3872)$ in proton-proton collision.

Ultra-peripheral heavy ion collisions are an ideal setup for different sorts
of exotic searches, and they could provide a key insight into a yet unanswered
question of strong interactions.

###### Acknowledgements.

We are grateful to F. Antinori, G. M. Innocenti and A. Uras for very fruitful
discussions and for encouraging the present study. A.E. is a Roger Dashen
Member at the Institute for Advanced Study, whose work is also supported by
the U.S. Department of Energy, Office of Science, Office of High Energy
Physics under Award No. DE-SC0009988. A.E. has also received funding by the
Swiss National Science Foundation under Contract No. 200020-169696 and through
the National Center of Competence in Research SwissMAP. The work of C.A.M. is
supported by the Swiss National Science Foundation (PP00P2_176884). A.P. has
received funding from the European Union’s Horizon 2020 research and
innovation program under the Marie Skłodowska-Curie Grant Agreement No.
754496.

##  Appendix A $\rho\to 2\pi$ and $\omega\to 3\pi$

The process $\rho\to\pi^{+}\pi^{-}$ is a $P$-wave decay. The matrix element is
given by

$\displaystyle\mathcal{M}=C_{\rho}\epsilon_{\rho}^{\mu}(p_{+}-p_{-})_{\mu}\,,$
(24)

where $\epsilon_{\rho}^{\mu}$ is the polarization vector of the $\rho$ and
$p_{+(-)}$ is the four momentum of the $\pi^{+(-)}$. For the decay rate, in
the rest frame of the $\rho$, one finds

$\displaystyle\Gamma_{\rho\to
2\pi}(s)\propto\int_{-1}^{1}d\cos\theta\,\frac{p}{\sqrt{s}}\sum_{\text{pol}}\left|\mathcal{M}\right|^{2}\,,$
(25)

where $s$ is the $\rho$ invariant mass, $\theta$ is the angle between the
$\rho$ quantization axis and the $\pi^{+}$ direction of flight in the $\rho$
rest frame, and
$p=\lambda^{1/2}\\!\left(s,m_{+}^{2},m_{-}^{2}\right)\big{/}2\sqrt{s}$ with
$m_{0,\pm}$ the pion masses. Performing the integral, we have

$\displaystyle\Gamma_{\rho\to 2\pi}(s)\propto\frac{p^{3}}{\sqrt{s}}\,.$ (26)

In the decay $\omega\to\pi^{+}\pi^{-}\pi^{0}$, the pair $\pi^{+}\pi^{-}$ has
an angular momentum $\ell=1$ and a relative angular momentum $\ell=1$ with the
$\pi^{0}$. The matrix element for this process is given by

$\displaystyle\mathcal{M}=C_{\omega}\varepsilon_{\mu\nu\rho\sigma}P_{\omega}^{\mu}\epsilon_{\omega}^{\nu}p_{+}^{\rho}p_{-}^{\sigma}\,,$
(27)

and for the decay rate one finds

$\displaystyle\begin{split}\Gamma_{\omega\to
3\pi}(s)&\propto\int_{(m_{+}+m_{-})^{2}}^{(\sqrt{s}-m_{0})^{2}}d\sigma\int_{-1}^{1}d\cos\theta^{+}\frac{\sqrt{\sigma}p}{s}\\\
&\quad\times\frac{q}{\sqrt{\sigma}}\sum_{\text{pol}}\left|\mathcal{M}\right|^{2}\,.\end{split}$
(28)

where $\sigma$ is the $\pi^{+}\pi^{-}$ invariant mass, $s$ is the $\omega$
invariant mass and $\theta^{+}$ is the angle between $\omega$ and $\pi^{+}$
directions of flight in the $\pi^{+}\pi^{-}$ rest frame. In addition,
$p=\lambda^{1/2}\\!\left(s,m_{0}^{2},\sigma\right)\big{/}2\sqrt{\sigma}$ and
$q=\lambda^{1/2}\\!\left(\sigma,m_{+}^{2},m_{-}^{2}\right)\big{/}2\sqrt{\sigma}$.

Neglecting irrelevant constants which cancel in Eq. (17), we find

$\displaystyle\Gamma_{\omega\to
3\pi}(s)\propto\int_{(m_{+}+m_{-})^{2}}^{(\sqrt{s}-m_{0})^{2}}d\sigma\frac{(\sqrt{\sigma}p)^{3}}{s}\frac{q^{3}}{\sqrt{\sigma}}\,.$
(29)

## References

  * Esposito  _et al._ (2017) A. Esposito, A. Pilloni, and A. D. Polosa, Phys. Rept. 668, 1 (2017), arXiv:1611.07920 [hep-ph] .
  * Guo  _et al._ (2018) F.-K. Guo, C. Hanhart, U.-G. Meißner, Q. Wang, Q. Zhao, and B.-S. Zou, Rev. Mod. Phys. 90, 015004 (2018), arXiv:1705.00141 [hep-ph] .
  * Olsen  _et al._ (2018) S. L. Olsen, T. Skwarnicki, and D. Zieminska, Rev. Mod. Phys. 90, 015003 (2018), arXiv:1708.04012 [hep-ph] .
  * Brambilla  _et al._ (2020) N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C.-P. Shen, C. E. Thomas, A. Vairo, and C.-Z. Yuan, Phys. Rept. 873, 1 (2020), arXiv:1907.07583 [hep-ex] .
  * Aaij  _et al._ (2020a) R. Aaij _et al._ (LHCb), Sci. Bull. 65, 1983 (2020a), arXiv:2006.16957 [hep-ex] .
  * Chao (1981) K.-T. Chao, Z. Phys. C 7, 317 (1981).
  * Heller and Tjon (1985) L. Heller and J. A. Tjon, Phys. Rev. D 32, 755 (1985).
  * Barnea  _et al._ (2006) N. Barnea, J. Vijande, and A. Valcarce, Phys. Rev. D 73, 054004 (2006), arXiv:hep-ph/0604010 .
  * Vijande  _et al._ (2007) J. Vijande, A. Valcarce, and J. M. Richard, Phys. Rev. D 76, 114013 (2007), arXiv:0707.3996 [hep-ph] .
  * Ebert  _et al._ (2007) D. Ebert, R. N. Faustov, V. O. Galkin, and W. Lucha, Phys. Rev. D 76, 114015 (2007), arXiv:0706.3853 [hep-ph] .
  * Berezhnoy  _et al._ (2012) A. V. Berezhnoy, A. V. Luchinsky, and A. A. Novoselov, Phys. Rev. D 86, 034004 (2012), arXiv:1111.1867 [hep-ph] .
  * Wu  _et al._ (2018) J. Wu, Y.-R. Liu, K. Chen, X. Liu, and S.-L. Zhu, Phys. Rev. D 97, 094015 (2018), arXiv:1605.01134 [hep-ph] .
  * Chen  _et al._ (2017) W. Chen, H.-X. Chen, X. Liu, T. G. Steele, and S.-L. Zhu, Phys. Lett. B 773, 247 (2017), arXiv:1605.01647 [hep-ph] .
  * Wang (2017) Z.-G. Wang, Eur. Phys. J. C 77, 432 (2017), arXiv:1701.04285 [hep-ph] .
  * Debastiani and Navarra (2019) V. R. Debastiani and F. S. Navarra, Chin. Phys. C 43, 013105 (2019), arXiv:1706.07553 [hep-ph] .
  * Richard  _et al._ (2017) J.-M. Richard, A. Valcarce, and J. Vijande, Phys. Rev. D 95, 054019 (2017), arXiv:1703.00783 [hep-ph] .
  * Anwar  _et al._ (2018) M. N. Anwar, J. Ferretti, F.-K. Guo, E. Santopinto, and B.-S. Zou, Eur. Phys. J. C 78, 647 (2018), arXiv:1710.02540 [hep-ph] .
  * Bedolla  _et al._ (2020) M. A. Bedolla, J. Ferretti, C. D. Roberts, and E. Santopinto, Eur. Phys. J. C 80, 1004 (2020), arXiv:1911.00960 [hep-ph] .
  * Karliner  _et al._ (2017) M. Karliner, S. Nussinov, and J. L. Rosner, Phys. Rev. D 95, 034011 (2017), arXiv:1611.00348 [hep-ph] .
  * Becchi  _et al._ (2020) C. Becchi, J. Ferretti, A. Giachino, L. Maiani, and E. Santopinto, Phys. Lett. B 811, 135952 (2020), arXiv:2006.14388 [hep-ph] .
  * Dong  _et al._ (2021a) X.-K. Dong, V. Baru, F.-K. Guo, C. Hanhart, and A. Nefediev, Phys. Rev. Lett. 126, 132001 (2021a), arXiv:2009.07795 [hep-ph] .
  * Cao  _et al._ (2021) Q.-F. Cao, H. Chen, H.-R. Qi, and H.-Q. Zheng, Chin. Phys. C 45, 093113 (2021), arXiv:2011.04347 [hep-ph] .
  * Liang  _et al._ (2021) Z.-R. Liang, X.-Y. Wu, and D.-L. Yao, (2021), arXiv:2104.08589 [hep-ph] .
  * Wang  _et al._ (2021) J.-Z. Wang, D.-Y. Chen, X. Liu, and T. Matsuki, Phys. Rev. D 103, 071503 (2021), arXiv:2008.07430 [hep-ph] .
  * Dong  _et al._ (2021b) X.-K. Dong, V. Baru, F.-K. Guo, C. Hanhart, A. Nefediev, and B.-S. Zou, (2021b), arXiv:2107.03946 [hep-ph] .
  * Maiani (2020) L. Maiani, Sci. Bull. 65, 1949 (2020), arXiv:2008.01637 [hep-ph] .
  * Maiani  _et al._ (2014) L. Maiani, F. Piccinini, A. D. Polosa, and V. Riquer, Phys. Rev. D 89, 114010 (2014), arXiv:1405.1551 [hep-ph] .
  * Baur  _et al._ (2002) G. Baur, K. Hencken, D. Trautmann, S. Sadovsky, and Y. Kharlov, Phys. Rept. 364, 359 (2002), arXiv:hep-ph/0112211 .
  * Bertulani  _et al._ (2005) C. A. Bertulani, S. R. Klein, and J. Nystrand, Ann. Rev. Nucl. Part. Sci. 55, 271 (2005), arXiv:nucl-ex/0502005 .
  * Baltz (2008) A. J. Baltz, Phys. Rept. 458, 1 (2008), arXiv:0706.3356 [nucl-ex] .
  * Grabiak  _et al._ (1989) M. Grabiak, B. Muller, W. Greiner, G. Soff, and P. Koch, J. Phys. G 15, L25 (1989).
  * Drees  _et al._ (1989) M. Drees, J. R. Ellis, and D. Zeppenfeld, Phys. Lett. B 223, 454 (1989).
  * Greiner  _et al._ (1993) M. Greiner, M. Vidovic, and G. Soff, Phys. Rev. C 47, 2288 (1993).
  * Ahern  _et al._ (2000) S. C. Ahern, J. W. Norbury, and W. J. Poyser, Phys. Rev. D 62, 116001 (2000), arXiv:gr-qc/0009059 .
  * Bertulani (2009) C. A. Bertulani, Phys. Rev. C 79, 047901 (2009), arXiv:0903.3174 [nucl-th] .
  * Goncalves  _et al._ (2013) V. P. Goncalves, D. T. Da Silva, and W. K. Sauter, Phys. Rev. C 87, 028201 (2013), arXiv:1209.0701 [hep-ph] .
  * Moreira  _et al._ (2016) B. D. Moreira, C. A. Bertulani, V. P. Goncalves, and F. S. Navarra, Phys. Rev. D 94, 094024 (2016), arXiv:1610.06604 [hep-ph] .
  * Knapen  _et al._ (2017) S. Knapen, T. Lin, H. K. Lou, and T. Melia, Phys. Rev. Lett. 118, 171801 (2017), arXiv:1607.06083 [hep-ph] .
  * Goncalves and Moreira (2019) V. P. Goncalves and B. D. Moreira, Eur. Phys. J. C 79, 7 (2019), arXiv:1809.08125 [hep-ph] .
  * Gonçalves and Moreira (2021) V. P. Gonçalves and B. D. Moreira, Phys. Lett. B 816, 136249 (2021), arXiv:2101.03798 [hep-ph] .
  * Chatrchyan  _et al._ (2013) S. Chatrchyan _et al._ (CMS), JHEP 04, 154 (2013), arXiv:1302.3968 [hep-ex] .
  * Aaij  _et al._ (2021a) R. Aaij _et al._ (LHCb), (2021a), arXiv:2109.07360 [hep-ex] .
  * Ablikim  _et al._ (2021) M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 126, 102001 (2021), arXiv:2011.07855 [hep-ex] .
  * Aaij  _et al._ (2021b) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 127, 082001 (2021b), arXiv:2103.01803 [hep-ex] .
  * Maiani  _et al._ (2021) L. Maiani, A. D. Polosa, and V. Riquer, Sci. Bull. 66, 1460 (2021), arXiv:2103.08331 [hep-ph] .
  * Aaij  _et al._ (2020b) R. Aaij _et al._ (LHCb), Phys. Rev. D 102, 092005 (2020b), arXiv:2005.13419 [hep-ex] .
  * Esposito  _et al._ (2021) A. Esposito, L. Maiani, A. Pilloni, A. D. Polosa, and V. Riquer, (2021), arXiv:2108.11413 [hep-ph] .
  * von Weizsacker (1934) C. F. von Weizsacker, Z. Phys. 88, 612 (1934).
  * Williams (1934) E. J. Williams, Phys. Rev. 45, 729 (1934).
  * Klein  _et al._ (2017) S. R. Klein, J. Nystrand, J. Seger, Y. Gorbunov, and J. Butterworth, Comput. Phys. Commun. 212, 258 (2017), arXiv:1607.03838 [hep-ph] .
  * Sakurai (1960) J. Sakurai, Annals of Physics 11, 1 (1960).
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Casalbuoni  _et al._ (1997) R. Casalbuoni, A. Deandrea, N. Di Bartolomeo, R. Gatto, F. Feruglio, and G. Nardulli, Phys. Rept. 281, 145 (1997), arXiv:hep-ph/9605342 .
  * Deandrea  _et al._ (2002) A. Deandrea, G. Nardulli, and A. D. Polosa, in  _3rd Workshop on Hard Probes in Heavy Ion Collisions: 3rd Plenary Meeting_ (2002) arXiv:hep-ph/0211431 .
  * Faccini  _et al._ (2012) R. Faccini, F. Piccinini, A. Pilloni, and A. D. Polosa, Phys. Rev. D 86, 054012 (2012), arXiv:1204.1223 [hep-ph] .
  * Gleisberg  _et al._ (2003) T. Gleisberg, F. Krauss, K. T. Matchev, A. Schalicke, S. Schumann, and G. Soff, JHEP 09, 001 (2003), arXiv:hep-ph/0306182 .
  * Mathieu  _et al._ (2020) V. Mathieu, A. Pilloni, M. Albaladejo, L. Bibrzycki, A. Celentano, C. Fernandez-Ramirez, and A. P. Szczepaniak (JPAC), Phys. Rev. D 102, 014003 (2020), arXiv:2005.01617 [hep-ph] .
  * Redlich  _et al._ (2000) K. Redlich, H. Satz, and G. M. Zinovjev, Eur. Phys. J. C 17, 461 (2000), arXiv:hep-ph/0003079 .
  * Rossi and Veneziano (2004) G. C. Rossi and G. Veneziano, Phys. Lett. B 597, 338 (2004), arXiv:hep-ph/0404262 .
  * Maiani  _et al._ (2005) L. Maiani, F. Piccinini, A. D. Polosa, and V. Riquer, Phys. Rev. D 71, 014028 (2005), arXiv:hep-ph/0412098 .
  * Brazzi  _et al._ (2011) F. Brazzi, B. Grinstein, F. Piccinini, A. D. Polosa, and C. Sabelli, Phys. Rev. D 84, 014003 (2011), arXiv:1103.3155 [hep-ph] .
  * Albaladejo  _et al._ (2020) M. Albaladejo, A. N. H. Blin, A. Pilloni, D. Winney, C. Fernandez-Ramirez, V. Mathieu, and A. Szczepaniak (JPAC), Phys. Rev. D 102, 114010 (2020), arXiv:2008.01001 [hep-ph] .
  * Artoisenet and Braaten (2010) P. Artoisenet and E. Braaten, Phys. Rev. D 81, 114018 (2010), arXiv:0911.2016 [hep-ph] .