# Baryons and tetraquarks using instanton-induced interactions

Nicholas Miesch<EMAIL_ADDRESS>Center for Nuclear Theory,
Department of Physics and Astronomy, Stony Brook University, Stony Brook, New
York 11794–3800, USA  Edward Shuryak<EMAIL_ADDRESS>Center for
Nuclear Theory, Department of Physics and Astronomy, Stony Brook University,
Stony Brook, New York 11794–3800, USA  Ismail Zahed
<EMAIL_ADDRESS>Center for Nuclear Theory, Department of Physics
and Astronomy, Stony Brook University, Stony Brook, New York 11794–3800, USA

###### Abstract

We analyze some aspects of the perturbative and non-perturbative interactions
in the composition of heavy quarkonia, heavy and light baryons ($ccc$ and
$uuu$ ones), as well as all charm tetraquarks ($cc\bar{c}\bar{c}$). Using the
hyper-spherical approximation and effective radial potentials (in 6 and 9
dimensions, respectively) we derive their spectra and wave functions. In all
of the cases, we focus on the splittings between the s-shell levels, which are
remarkably insensitive to the quark masses, but proportional to the effective
interaction potentials. We use the traditional Cornell-like potentials, and
the non-perturbative instanton-induced static potentials, from correlators of
two, three and four Wilson lines, and find rather satisfactory description of
spectra in all cases.

##  I Introduction

One of the goals of the present paper is to reassess the information one can
extract from hadronic spectroscopy, on quark-quark and quark-antiquark
interactions in multi-quark system, baryons and tetraquarks. Of course, we do
not attempt to cover this vast field, but rather focus on some particular
(flavor-symmetric) objects. Furthermore, we focus on specific features of
spectroscopy that are most sensitive to the inter-quark interactions (and less
sensitive to the quark masses and additive constants in the potentials).

As examples of these observables, we use the splittings between the $1S,2S,3S$
states ( S-shell, as usual means zero angular momentum), and evaluate the
splittings between the p-shell states (negative parity). We start in section
II with a very traditional issue, in relation to the central quarkonium
potential in charmonium and bottomonium. In spite of the large mass difference
between $c,b$ quarks, these splittings in are very similar. While this fact by
itself is widely known, we elaborate in more detail on its sensitivity to the
exact shape of the potentials. More specifically, we use comparatively the
Cornell (Coulomb-plus-linear) potential, the Martin potential, and the
instanton-induced potential we derived in Shuryak and Zahed (2023a) (our first
paper of the light-front series) .

In sections II.1 and III.1 we discuss the flavor-symmetric baryons ($ccc,uuu$)
and $cc\bar{c}\bar{c}$ tetraquarks, respectively. As a method, we use the so
called hyperdistance (or hypercentral, or K-harmonics) approximation, by
reducing quantum mechanics (in 6 dimensions for baryons and 9 for tetraquarks)
to a single radial Schrodinger equation. The method has been developed in the
1960’s, in nuclear physics, for the quantum mechanical treatment of light
nuclei, see e.g. Simonov (1966). One well-known four-body example is
${}^{4}He$, for which this approach reproduces its binding $\approx 28\,MeV$
from conventional nuclear forces. Furthermore, in a relatively recent paper
Shuryak and Torres-Rincon (2019) on 4-nucleon clustering in heavy ion
collisions at finite temperatures, it was found that the same hyperdistance
equation also correctly predicts the $second$ (2S) level of ${}^{4}He$, with
small binding and close to its experimentally known location. This work
increased our trust in this approach, at least for the splittings of the
s-shell states.

The applications of the hyperdistance approximation to various multiquark
states have also been done over the years, see e.g. Badalian  _et al._ (1987);
Chernyshev  _et al._ (1996); Ferraris  _et al._ (1995); Gandhi  _et al._
(2018); Brambilla  _et al._ (2022) and references therein. With the recent
revival of hadronic spectroscopy due to the discovery of multiple tetra- and
penta-quark states, there has been renewed interest in few-body quantum
mechanics in general, and in the hyperdistance approximation in particular.

In sections III.1 and C, we start our analyses using Cornell-like potentials,
and focus chiefly on the level splittings of the S-shell states. Much like in
quarkonia, we found that the splittings are insensitive to the quark mass
values: for $ccc$ and $uuu$ baryons, these splittings are nearly identical,
see Fig.3. What they are sensitive to is the overall strength of the confining
force in multiquark systems, to which they are roughly proportional. (This is
highly nontrivial, since the wave functions and even sizes of these systems
are vastly different!)

The main physics issues this paper is aimed at are the following questions:
Can the interactions inside few-body hadrons (such as baryons and tetraquarks
specifically) be approximated by a sum of two-body ones? Do they only depend
on the hyper-distance variable (related to the sum of distances between all
constituents squared)? If so, are the shapes of the corresponding effective
potentials similar to the Coulomb+linear potentials used for quarkonia? What
are the magnitude of their respective coefficients, for 3 and 4 quarks? What
is the order of magnitude of the corrections stemming from the lowest hyper-
spherical approximation?

Static interquark potentials have well known expressions in terms of the
correlators of Wilson lines. Lattice gauge theories in fact started with the
evaluation of correlated Wilson lines $\langle WW^{\dagger}\rangle$, and the
shape of the $\bar{q}q$ static potential. The literature on this subject is
vast and cannot be covered in this work. Yet the lattice studies of
correlators of three or four Wilson lines are surprisingly limited, with
results limited to simple geometries for forces in baryons Koma and Koma
(2017) and in tetraquarks Bicudo  _et al._ (2017). For the baryons, three
quarks are set at the corners of triangles of different shapes and sizes,
while for tetraquarks four quarks are set at the corners of rectangles
$r\times R$ with variable $r,R$. Unfortunately, no attempts to project the
results on the hyperdistance variable were made so far.

While numerical lattice simulations are first-principle-based approaches to
vacuum ensemble of gauge fields, some simple models for the vacuum have been
developed over the years. Instantons are 4-d spherical solutions in Euclidean
gauge theory vacua, describing tunneling between configurations with different
Chern-Simons numbers. They are the basis of semiclassical approximation and
trans-series in QFT. Here, we will not go into the theoretical details of the
model, but just just state that we will use what we refer to as the “dense
instanton vacuum model". Unlike earlier version, the “dilute instanton vacuum"
Shuryak (1982) (which included only isolated instantons and focused on their
fermionic zero modes and chiral symmetry breaking, for review see Schäfer and
Shuryak (1998)), it also includes close instanton-antiinstanton pairs which
contribute to Wilson lines by their fields. The two parameters of the model
are the 4d density $n_{I+\bar{I}}\approx 7\,{\rm fm}^{-4}$ of instantons and
antiinstantons, and their mean size $\rho\approx\frac{1}{3}\,{\rm fm}$. For
discussion of these parameters and its relation to static quarkonium potential
see Shuryak and Zahed (2023a). Throughout the same parameters will be used in
the evaluation of the effective potentials in all systems considered, e.g.
$\bar{q}q,qqq,qq\bar{q}\bar{q}$. So the resulting potentials are essentially
parameter free. All we do is to generalize the calculation of Wilson line
correlators $\langle WW^{\dagger}\rangle$, to $WWW$ and
$WWW^{\dagger}W^{\dagger}$ pertinent correlators. We are satisfied to find
that the resulting potentials do indeed reproduce the pertinent hadron spectra
in all cases.

##  II Central potential in quarkonia

Since the early 1970’s, when heavy $c,b$ quarks were discovered, the
nonrelativistic description of quarkonia $\bar{c}c,\bar{b}b$ states became the
main pillars of hadronic spectroscopy. And yet, we decided to start again from
quarkonia, to test how well the proposed potentials work.

A good motivation is to also include current data on bottomonium states, as
listed in the review of particle physics 2022 Workman  _et al._ (2022). We
focus below on the S-level $splittings$, insensitive to quark masses and
overall constant terms in the potentials. One such input is the splitting
between the spin-averaged charmonium masses

$M_{spin-average}=(3M_{\psi}+M_{\eta_{c}})/4$ (1)

in the 2S and 1S shell

$M_{spin-average}^{2S}-M_{spin-average}^{1S}\approx 605\,MeV$ (2)

In the bottomonium family we now have $\Upsilon$ states listed up to 4S, but
not enough $\eta_{c}$ states to perform spin average for all of them. Hence,
we will use the Upsilon mass splittings instead. (Spin-dependent forces are
not yet included, but those are $\sim 1/m_{b}^{2}$ and rather small in
bottomonia.) These data points are shown as red pentagons in Fig.1.

Figure 1: Mass splittings in GeV, of the $\Upsilon(nS)-\Upsilon(1S),n=2,3,4$
(upper plots) and spin-averaged splittings in $\bar{c}c$ systems (lower plots)
shown versus $n-1$ by red pentagons. The blue circles correspond to solution
of the Schrodinger equation using the Cornell potentials, the green circles
are for the Martin potential, and the black disks are for the instanton-
induced potential explained in the text.

The precise shape of static potential, describing known quarkonium masses, is
an old problem, discussed many times since e.g. the first reviews such as
Grosse and Martin (1980). The classic Cornell potential is the sum of the one-
gluon exchange Coulomb term, and the linear confining term. For definiteness,
we use its version with coefficients

$V_{Cornell}=-{0.641\over r}+0.113\cdot r$ (3)

(Throughout, all dimensions are in GeV units, e.g. the distance $r$ is in
$GeV^{-1}$, and the string tension in $GeV^{2}$.) The two parameters in the
potential are fitted so as to reproduce the three known bottomonium splittings
$2S-1S,3S-1S,4S-1S$, shown in the upper plot by the blue circles. Note that it
does not reproduce well the charmonium splittings as shown in the lower pannel
of Fig. 1.

Let us remind the reader that potentials which grow as smaller power of $r$,
display a smaller growth of energy with $n$. ( Only an oscillator $V\sim
r^{2}$ has $E_{n}\sim n$ linearly.) Early studies argued that static
potentials $V\sim r^{\alpha}$ with $\alpha<1$ can describe spectra better. The
so called Martin potential fitted to (then available) data, see e.g. Grosse
and Martin (1980)

$V_{Martin}(r)=-8.054+6.870\cdot r^{0.1}$ (4)

(all in GeV) has a very small power, and therefore is close to a logarithmic
potential, for which it is known that the 1S-2S splittings of the charmonium
and bottomonium systems would be equal. And indeed, they are almost the same.
More specifically, the spin-weighted average (1S, 2S) masses are (3068.65,
3673.95) MeV for charmonium and (9444.9, 10017.2) MeV for bottomonium, so
M(2S)-M(1S) = (605.3,572.3) MeV for ($c\bar{c},b\bar{b}$).

As one can see from the lower pannel of Fig.1, with the Martin potential
predictions for splittings are shown by green circles, it indeed does give a
very accurate number for the first charmonium splitting (to which the
potential parameter was originally fitted). Yet its predictions for $\bar{b}b$
Upsilon states (upper plot) is rather far from the data points.

To summarize our preliminary discussion: the Cornell potential is still
perhaps the best choice, reproducing well these particular observables both
for $\bar{c}c$ and $\bar{b}b$ systems. The potentials with smaller power, such
as the Martin potential, do fit charmonium data but are not very good for
bottomonia.

###  II.1 Instanton-induced forces in mesons

Understanding the hadronic spectroscopy can be viewed as only a step towards
the even more ambitious goal of understanding the structure of the (Euclidean)
gauge fields in the QCD vacuum. It may be approximate and schematic, yet it
allows for practical evaluation of all kinds of non-local correlations
functions of nonperturbative gauge fields, without supercomputers. The static
potentials are the nonlocal correlators of fields inside the Wilson lines, the
path-ordered exponents of color fields.

Instantons are nontrivial gauge configurations, being the basis of the
semiclassical theory. Their properties in the QCD vacuum, chiral symmetry
breaking and physical effects on mesonic and baryonic point-to-point
correlators were reviewed in Schäfer and Shuryak (1998). Recent discussion
including the role of instanton-antiinstanton configurations on the lattice
and for effective potentials can be found in the first paper of this series
Shuryak and Zahed (2023a).

The contribution of instantons to the quark potential is an old subject,
pioneered by Callan  _et al._ (1978), but its practical application to central
and spin-dependent forces between quarks attracted very little attention
Chernyshev  _et al._ (1996). In our recent series of papers Shuryak and Zahed
(2023a, b, 2021), a novel version of the instanton model was proposed for
that. It incorporates lattice data on close instanton-antiinstanton pairs and
was shown to generate the $\bar{q}q$ static potential which, while not
confining per se, is close to the Cornell potential up to about 1 fm
distances.

We then extended this model to calculation of the spin-dependent $\bar{q}q$
forces, and static quark-quark potentials in diquarks and baryons. Those can
be expressed as Wilson lines, decorated by two field strengths. In the cases
of spin-spin and tensor forces, those fields are gluo-magnetic. (Note that
attributing confinement to gluo-electric flux tubes, one is faced with the
problematic of the spin-dependent forces. Indeed, there is no place as such
for magnetic fields, and the only contributution is via the Thomas’
precession. The instantons are self-dual, with gluomagnetic and gluoelectric
components being equal, and thus they do generate significant spin-dependent
forces.)

We now briefly recall some general points in the derivation of the $\bar{q}q$
potentials. The Wilson lines go along the world-path of static charges along
the Euclidean time. The produced $W$ operators are unitary color matrices,
describing rotations of quark color vectors by certain angles and in certain
directions, around which these rotations are to be made. Let the 3-vectors
representing distances from the quark locations $\vec{r}_{i},i=1...N$ to the
instanton center be $\vec{\gamma}_{i}=\vec{r}_{i}-\vec{y}$ (amusingly it acts
as a vertex center). The color rotation angles are

$\alpha_{i}=\pi\bigg{(}1-{\gamma_{i}\over\sqrt{\gamma_{i}^{2}+\rho^{2}}}\bigg{)}$
(5)

where $\rho$ is the instanton size. Far from the instanton,
$\gamma_{i}\gg\rho$, this angle vanishes. At the instanton center
$\gamma_{i}=0$, the angle is maximal $\alpha_{i}=\pi$. The $SU(2)$ part of the
Wilson line can be expressed using Pauli (rather than Gell-Mann) matrices
$\tau^{A},A=1,2,3$

$\displaystyle{\bf W}^{a}_{lb}$ $\displaystyle=$ $\displaystyle{\rm exp}(\pm
i\vec{\tau}\cdot\vec{n}_{l}\alpha_{l})$ $\displaystyle=$
$\displaystyle\big{(}cos(\alpha_{l})\,{\bf\hat{1}}-i(\vec{\tau}\cdot\vec{n}_{l})\,sin(\alpha_{l})\big{)}^{a}_{b}$

where the unit vectors are $n_{i}=\vec{\gamma}_{i}/\gamma_{i}$, and
$a,b=1,2;l=1...N_{q}$. The third color is not affected by the instanton
fields, so this $2\times 2$ matrix should be complemented trivially by the
unit value in the third row.

The integral over (large) Euclidean time $\tau\in[-T/2,T/2]$ may include some
number of instantons $N_{inst}$. Summing over $N_{inst}$ leads to the
exponentiation of the correlator. The exponent is extensive in the Euclidean
time $T$ times the potential

$V(r)=n_{I+\bar{I}}\int d^{3}y{1\over N_{c}}{\rm Tr}\big{(}\hat{\bf
1}-W(0)W^{\dagger}(\vec{r})\big{)}$ (7)

For the plot of it see Fig.5 of Shuryak and Zahed (2023a). Below these
formulae will be generalized for 3 and 4 quark lines. These calculations can
later be straightforwardly generalized to other multi-quark hadrons. We plan
to calculate the spin-dependent forces in the same settings in forthcoming
publications.

The instanton-induced potential together with the Coulomb term, can be put
into the Schrodinger equation, to derive the quarkonium states. Let us now
return to Fig.1 for the mass differences between the excited $2S,3S,4S$
states, from the $1S$ ground state of $\bar{b}b,\bar{c}c$ quarkonia. The
splittings for the instanton-based potential are shown in Fig.1 by black
discs. For the 2S-1S and 3S-1S splittings of bottomonia the results match the
data quite well, and only for the 4S state they visibly deviate from the data,
by about 10% or so. For charmonia this potential leads to smaller level
splittings. Overall, its accuracy is in between the Cornell and Martin
potentials. Note however, that unlike those, fitted to reproduce quarkonia
levels, the instanton parameters (sizes and density) were not fitted to
spectroscopy, the instanton density and mean size taken directly from lattice
data about vacuum topology.

Figure 2: Correlator of three Wilson lines in baryons, connected by a vertex
with Levi-Cevita symbol.

##  III Baryons: hyper-spherical approximation and instanton-induced potential

###  III.1 Baryons in hyperspherical approximation

We started this paper discussing heavy quarkonia, and we will end discussing
full-charm tetraquarks. However, it is also natural to address heavy-quark
baryons, $ccc$ or $bbb$, and anything in between. Unfortunately, these hadrons
still elude experiments, so we will also include calculations for baryons
composed of light quarks. In doing so, we will assign the same flavor for all
three quarks, i.e. $uuu=\Delta^{++}$ states. (There is not enough data on
$sss=\Omega^{-}$ spectra.)

We will use the hyper-spherical approximation, also known as the method of
K-harmonics. The main idea is that the lowest states depend mostly on the
radial hyperspherical coordinate in 6-d space of Jacobi coordinates

$R_{6}^{2}=\vec{\rho}^{2}+\vec{\lambda}^{2}$ (8)

with a small admixture of components with angular dependences. In standard
terms, we study 6-d s-shell states. Their wave function can be obtained by
solving a single $radial$ Schroedinger equation. The method has been used for
baryons over the years, see e.g. Guimarães  _et al._ (1981) for heavy baryons
of interest. The kinematical details of this approximation are given in
Appendix A,B.

The radial Schrodinger equation for the $S$ states has the form

$\bigg{(}-u"+{15\over 4R_{6}^{2}}u\bigg{)}{1\over 2M}+(V_{6}-E_{Ns})u=0$ (9)

where the reduced wave function $u(R_{6})$ is related to the radial wave
function by

$\psi(R_{6})=u(R_{6})/R_{6}^{5/2}$ (10)

Note, that $R_{6}^{5}$ gets absorbed in the volume element, so that $u$ has
the usual normalization, as in one dimension $\sim\int dR_{6}|u(R_{6})|^{2}$.
This reduction eliminates the first derivative from the Laplacian, but adds an
extra “quasi-centrifugal" term. The radial projection of the potential on
$R_{6}$ is defined in the Appendix B.

The mass splittings of the lowest s- and p-states, $ccc$ and $uuu$ baryons,
following numerically from (9) are shown in Fig.3. The binary $QQ$ potential
is taken as $1/2$ of the Cornell potential (used above for $\Upsilon$ family).
We do not show the absolute masses but focus instead on the level splittings.
Experience shows that the splittings have relatively little dependence on the
choice of the quark mass, but depend primarily on the interaction potential.
Indeed, as shown in Fig.3, the plots for $ccc$ and $uuu$ baryons look nearly
identical, in spite of the quark masses differing by more than a factor 3.
(This is a reminder of the nearly identical splittings in the charmonium and
bottomonium families emphasized above: those are also different by a similar
ratio of 3 in absolute mass scale.)

More specifically, the splittings for the s-shell $ccc$ baryons are calculated
to be

$M_{ccc}^{2S}-M_{ccc}^{1S}\approx 391,\,\,M_{ccc}^{3S}-M_{ccc}^{1S}\approx
685\,MeV$

and for the $uuu=\Delta^{++}$ states they are

$M_{uuu}^{2S}-M_{uuu}^{1S}\approx 438,\,\,M_{uuu}^{3S}-M_{uuu}^{1S}\approx
804\,MeV$

The splittings for $uuu=\Delta(3/2^{+})$ should be compared to experimental
values of these splittings

$M_{uuu}^{2S}-M_{uuu}^{1S}\approx 278,\,\,M_{uuu}^{3S}-M_{uuu}^{1S}\approx
688\,MeV,$

(from the particle data tables 2022, using masses from the Breit-Wigner fits).
While the deviations between the calculated and experimental values are
comparable to the shifts expected from the spin-spin and ’t Hooft forces,
which are so far ignored in the calculation, its sign would be opposite,
shifting 1S down more than higher states. (Also, these are relativistic
corrections in $\sim 1/M^{2}$, and would of course be much less important for
$c$ quarks.) We therefore conclude that the effective potential for baryons
used (Ansatz A), is somewhat stronger than it is in reality. We will elaborate
on this point in the next section.

Figure 3: Masses of $ccc$ baryons (upper plot) and $uuu=Delta^{++}$ (lower
plot) $M(N,L)$, for $L=0$ and $L=1$ shells, counted from the ground states.

We also show the corresponding wave functions in Fig.4. Unlike the splittings,
those show radical difference between the $uuu$ and $ccc$ baryons. Indeed, the
corresponding hadrons are of quite different sizes, with the $ccc$ ones much
smaller. Also, all of them are somewhat unusual for the ground states. In
spite of the enhanced Coulomb forces in the 6d systems, the quasi-centrifugal
potential leads to wave functions very strongly suppressed at small hyper-
radial distances.

Figure 4: Wave functions $u(R_{6})$ of $ccc$ (upper) and $uuu$ baryons (lower)
for $L=0$ (solid lines) and $L=1$ (dashed lines), versus the hyper-radial
distance $R_{6}\,(GeV^{-1})$.

These wave functions can be used for getting matrix elements of various
operators. R.m.s radii of three s-states of $\Delta$ are

$\sqrt{\langle R_{6}^{2}\rangle}=4.97,\,7.53,\,9.72\,GeV^{-1}$ (11)

To summarize this section: we have tested a popular assumption, that the non-
perturbative forces between quarks in baryons can be described by a sum of
binary forces reduced by the same factor 1/2 (as perturbative ones), as
compared to mesonic $\bar{Q}Q$ potential (the so called Ansatz A). We found
that it leads to splittings between s-shell states that are somewhat larger
than observed in the $uuu=\Delta^{++}$ baryons. A possible fix is to introduce
an arbitrary factor f for the confining force, and then fit it experimentally
to the splittings if $f<1/2$.

###  III.2 Instanton-induced potentials for baryons

We now derive the three-quark effective potential along the radial 6d
coordinate $V_{3}(R_{6})$, using the instanton model of the QCD vacuum.

We start by briefly recalling the approach we took in Ref.Shuryak and Zahed
(2023b). If no assumptions about the colors of the three quarks are made, the
color indices of $W$ remain open. Furthermore, to model a color-isotropic
vacuum in which instantons are at random 2-d planes in the $SU(3)$ color
space, we introduced random matrices $U\in SU(3)$ which rotate the instanton
fields from their standard $SU(2)$ plane, $W\rightarrow UWU^{\dagger}$. All
six matrices $U$ were then averaged using the invariant Haar measure of the
SU(3) group,

$\displaystyle\int dU(U^{a_{1}}_{i_{1}}{\bf W}^{i_{1}}_{1j_{1}}U^{\dagger
j_{1}}_{b_{1}})(U^{a_{2}}_{i_{2}}{\bf W}^{i_{2}}_{2j_{2}}U^{\dagger
j_{2}}_{b_{2}})(U^{a_{3}}_{i_{3}}{\bf W}^{i_{3}}_{3j_{3}}U^{\dagger
j_{3}}_{b_{3}})$

before convolution of the external color indices. In Shuryak and Zahed (2023b)
we used the so called 6-$U$ “Weingarten formula" for this integral, and the
Wilson lines running through instanton fields, leading to a “generic triple-
quark potential" of the form

$\displaystyle V$ $\displaystyle=$ $\displaystyle n_{I+\bar{I}}\int
d^{3}z\bigg{[}(1-c_{1}c_{2}c_{3})\delta^{a_{1}}_{b_{1}}\delta^{a_{2}}_{b_{2}}\delta^{a_{3}}_{b_{3}}$
(13) $\displaystyle+$ $\displaystyle\frac{9}{8}\,c_{1}s_{2}s_{3}n_{2}\cdot
n_{3}\delta^{a_{1}}_{b_{1}}\bigg{(}\frac{1}{2}\lambda_{2}^{B}\bigg{)}^{a_{2}}_{b_{2}}\bigg{(}\frac{1}{2}\lambda_{3}^{B}\bigg{)}^{a_{3}}_{b_{3}}$
$\displaystyle+$ $\displaystyle{\rm 2\,perm.}\bigg{]}$

with trigonometric functions involving color rotation angles, that depend on
the 3-dimensional distances $\vec{\gamma}_{i},i=1,2,3$ between the location of
the Wilson lines $\vec{r}_{i}$, and the instanton center $\vec{y}$

$\displaystyle c_{i}$ $\displaystyle\equiv$ $\displaystyle{\rm
cos}(\alpha_{i}),\,\,s_{i}\equiv{\rm sin}(\alpha_{i})$
$\displaystyle\gamma^{2}_{i}$ $\displaystyle=$
$\displaystyle(\vec{r}_{i}-\vec{y})^{2},\,\,\,\vec{n}_{i}=\vec{\gamma}_{i}/|\vec{\gamma}_{i}|$
(14)

Since the standard instanton is an $SU(2)$ solution, it does not interact with
a quark with the “third" color, orthogonal to the $SU(2)$ plane. The third
color was ignored in the formula above. (Below we will extend the $2\times 2$
matrix by adding $W_{A}(3,3)=1$ and other zeros to the $3\times 3$ matrix.)
Yet, for random orientations of the instanton induced by a rotation matrix
$U$, and for arbitrary colors of the three quarks, the potential is
generically a three-body potential, $not$ a sum of two-body ones.

Here, we perform a different calculation, limiting it to the case of three
quarks making baryons. The corresponding color wave functions are
$C_{ijk}=\epsilon_{ijk}/\sqrt{6}$. Putting those before and after the Wilson
lines with the instanton-induced rotations (III.2), the result simplifies
considerably. The randomizing matrices $U$ acting on the Levi-Civita tensor
cancel out (unit determinant). The resulting potential is therefore reduced to
what we will call a “determinantal form"

$\displaystyle V_{3}$ $\displaystyle=$
$\displaystyle(n_{I+\bar{I}}\rho^{3})\int{d^{3}y\over\rho^{3}}$
$\displaystyle\times$
$\displaystyle\bigg{(}1-\sum{\epsilon_{a_{1}a_{2}a_{3}}\over\sqrt{6}}W_{1}^{a_{1}b_{1}}W_{2}^{a_{2}b_{2}}W_{3}^{a_{3}b_{3}}{\epsilon_{b_{1}b_{2}b_{3}}\over\sqrt{6}}\bigg{)}$

Note that far from the instanton the color rotation angles of the Wilson lines
(5) vanish. All $W_{ab}\rightarrow\delta_{ab}$ and therefore the integrand
goes to zero, with a converging integral. Note also that we rewrote the
integral over the instanton location $d^{3}y$ as dimensionless, with a
coefficient $(n_{I+\bar{I}}\rho^{3})\approx 51\,MeV$.

Although this expression has the form of a generic three-body interaction of
three Wilson lines, for Wilson lines in the field of an instanton, we (with
some surprise) found that it does become the sum of three $binary$
interactions for 12,13,23 pairs, each proportional to the “relative rotation"

$cos(\alpha_{i})cos(\alpha_{j})+sin(\alpha_{i})sin(\alpha_{j})(\vec{n}_{i}\cdot\vec{n}_{j})$
(16)

The same combination appears in the $\bar{Q}Q$ potentials.

In Ref.Shuryak and Zahed (2023b) we evaluated the potential for different
triangles, and compared the results to similar triangles studied on the
lattice. Here we analyse baryons in the hyperspherical approximation, so we
just take one equilateral triangle with variable size. The quarks are put at
the locations $(r,0,0),(0,r,0),(0,0,r)$, so the distances between them are
$r_{ij}=\sqrt{2}r$. Using the definitions of the Jacobi coordinates, one can
see that $R_{6}=\sqrt{2}r$ as well.

The effective three-quark potential $V_{3}$ in (III.2) for this equilateral
triangle, is shown in Fig. 5.

Figure 5: The points show the instanton-induced effective potential in baryons
$V(R_{6}),(MeV)$ vs $R_{6}\,(fm).$. The tension force is of the order of 1
GeV/fm (shown by a line for comparison).

Here we used $n_{I+\bar{I}}=7\,{\rm fm}^{-4},\rho=1/3\,{\rm fm}$,
corresponding to the “dense instanton vacuum model" we used for the meson
potential. Remarkably, the effective $QQQ$ potential along $R_{6}$, is quite
similar in shape to the meson $Q\bar{Q}$ potential.

Using this potential (with the same one-gluon exchange as in previous section)
in a hyperspherical approximation, gives the following splittings in the $ccc$
and $uuu$ baryons (GeV)

$\displaystyle M_{ccc}^{2S}-M_{ccc}^{1S}$ $\displaystyle=$ $\displaystyle
450,\,\,\,M_{ccc}^{3S}-M_{ccc}^{1S}=675,MeV$ $\displaystyle
M_{uuu}^{2S}-M_{uuu}^{1S}$ $\displaystyle=$ $\displaystyle
232,\,\,\,M_{uuu}^{3S}-M_{uuu}^{1S}=394\,MeV$

While the results for the $ccc$ baryons are very close to those found in the
previous section for “ansatz A", those for the $uuu$ baryons are reduced
significantly, and are now much closer to the experimental value of the first
splitting. This means that the instanton-induced potential is a better
representation of the force at $R_{6}\sim 1\,fm$. The splitting of the 3S
state is too small because the instanton-induced potential at $R_{6}>1{\rm
fm}$ is flat. In order to get it in agreement with the data one would need a
potential in between linear and constant, perhaps growing as a slower function
of the distance.

##  IV The all-charm tetraquarks

###  IV.1 Masses and models

Experimentally these type of tetraquarks have been discovered by LHCb Aaij
_et al._ (2020) and then by CMS Zhang and Yi (2022) and ATLAS Xu (2023) in
$J/\psi,J/\psi$ and related channels. These states are named X(6550), X(6900),
X(7280). For definiteness, we will use the ATLAS masses

$M_{1}=6620,\,\,\,M2=6870,\,\,M3=7220\,(MeV)$ (17)

all with errors $\pm 30\,MeV$. These states provide a good opportunity to
check our understanding of both the multi-particle dynamics and inter-quark
forces.

Among the many models put forth for the tetraquarks, we note the diquark-
antidiquark formulation Karliner and Rosner (2020); Giron and Lebed (2020), in
which a 4-quark system is represented as a binary state of diquarks with
binary potentials, equal to that in quarkonium $\bar{Q}Q$. This would be true
if distances between quarks in diquarks (to be called $r$ below) are much
smaller than distances between diquarks $R$.

Since they are supposed to be comparable for $r\sim R$ in the lowest
tetraquark states, we cannot neglect any of the 6 pair-wise interactions in
the four-body compounds. In this section we use the four-body approach, with
all six binary potentials included. Lattice studies of four-quark static
potentials, with quenched and dynamical quarks, have been performed in Bicudo
_et al._ (2017). (Again, the dependence on the hyper-distance could be
determined, but was not.)

###  IV.2 The lowest all-charm tetraquark states in the hyper-distance
approximation

The use of the hyper-distance approximation for all-charm teraquarks was
pioneered in Badalian  _et al._ (1987) (see also the recent study in Badalian
(2023)). The main idea is that for the lowest states, one may assume that they
are mostly $L=0$ s-shell, with only a small d-shell admixture. If so, we only
need to solve a one-dimensional radial Schroedinger equation.

For four particles, there are three Jacobi coordinates, or 9-dimensions. We
note that in our definition of these coordinates and hyper-distance (see
Appendix D.1), as well as the ensuing radial equation, we differ with Badalian
(2023). Also, our 9-dimensional wave function is written with a different
power

$\psi(R_{9})=u(R_{9})/R_{9}^{4}$ (18)

which allows the elimination of the first derivative term in the 9-d
Laplacian, and also puts the normalization into the 1-dimensional form
$\sim\int dR_{9}|u(R_{9})|^{2}$. The corresponding radial Hamiltonian is then

$\bigg{(}-\frac{d^{2}}{dR_{9}^{2}}+{12\over R_{9}^{2}}\bigg{)}{1\over
2M}+(4M+V)$ (19)

in which the second term in the bracket is another “quasi-centrifugal" term,
following from the 9-dimensional Laplacian.

The potential includes the color factors derived in Appendix D.1, which happen
to be the same for $\bar{3}3$ and $6\bar{6}$ diquark color structures. The
other factors stemming from the angular projection of the interparticle
distances $r_{ij}$ to the hyper-distance $R_{9}$, are discussed in Appendix
D.1 for the Cornell-type potentials. Here we use the Cornell potential as used
for the Upsilons above, with

$V=2*(0.773*R_{9}*0.113-1.55*0.64/R_{9})$ (20)

As for the mesons and baryons we discussed above, and in order not to get
bogged down by the issue of the quark mass value and the constant in the
potentials, we focus on the level $splittings$. Solving the Schroedinger
equation (19), we find for the of the s-shell splittings

$\displaystyle M(2S)-M(1S)$ $\displaystyle\approx$ $\displaystyle 370\,MeV$
$\displaystyle M(3S)-M(1S)$ $\displaystyle\approx$ $\displaystyle 689\,MeV$
$\displaystyle M(4S)-M(1S)$ $\displaystyle\approx$ $\displaystyle 975\,MeV$
(21)

and for the p-shell splittings

$\displaystyle M(1P)-M(1S)$ $\displaystyle\approx$ $\displaystyle 587\,MeV$
$\displaystyle M(2P)-M(1S)$ $\displaystyle\approx$ $\displaystyle\,879\,MeV$
(22)

although the latters are not good candidate states for decaying into $J/\psi$
pair.

From experimentallly reported (ATLAS fit) masses of the three
$cc\bar{c}\bar{c}$ resonances, we find for the splittings

$\displaystyle M(2)-M(1)$ $\displaystyle\approx$ $\displaystyle 250\,MeV$ (23)
$\displaystyle M(3)-M(1)$ $\displaystyle\approx$ $\displaystyle 600\,MeV$

They are in fair agreement with the S-shell levels we calculated. As
experience shows, this means that the overall magnitude of the forces in all
the 6 quark pairs, is in reasonable agreement with the reported observations.

The corresponding wave functions are shown in Fig. 6. We note again, the very
strong suppression at small hyper-distance, which is due to the quasi-
centrifugal term in the Laplacian. (In Badalian (2023) it is noted that it
coincides with the centrifugal term for the $L=3$ or F-shell quarkonia.) The
r.m.s. sizes of these states are therefore rather large in $\,GeV^{-1}$,
respectively.

Figure 6:  The wave functions $u(\rho)$ for 1S,2S,3S states of all-charm
tetraquarks in hyperspherical approximation, shown by solid black,blue and red
lines, respectively. The dahsed lines are p-shell states.

##  V Instanton-induced forces in tetraquarks

The problem with the correlators of four Wilson lines, like for three ones,
can either be set for (i) arbitrary colors of the quarks, or (ii) for those in
a particular color wave function. In the former case, the correlator has color
indices open and each Wilson line is rotated by random matrices
$UWU^{\dagger}$, with subsequent averaging over the $U$ using the invariant
measure of the $SU(3)$ group. The general case of eight $U$ is not used, but
still discussed in Appendix E.

The color wave functions built from diquarks are discussed in Appendix D.1.
Applying those at early and late times, we get certain products of delta
indices. The convolution of the color indices is shown by red lines in Fig.7.
The analytic expression for the $\bar{3}3$ correlator is

$\displaystyle K_{\bar{3}3}(r,R)=$ (24) $\displaystyle\bigg{<}\sum
C_{\bar{3}3}^{c^{\prime}_{1},c^{\prime}_{2},c^{\prime}_{3},c^{\prime}4}W_{c^{\prime}_{1}}^{c_{1}}W_{c^{\prime}_{2}}^{c_{2}}W^{c^{\prime}_{3}}_{c_{3}}W^{c^{\prime}_{4}}_{c_{4}}C_{\bar{3}3,c_{1},c_{2},c_{3},c_{4}}\bigg{>}_{y}$

where the color wave functions are given in (D.1). The sum is assumed over all
color indices, and the averaging corresponds to the 3-dimensional integration
over the instanton center position $\vec{y}$. Note that the anti-quark $W$
(assumed here to be 3,4) are shown as transposed, yet by including the
opposite color charges this is equivalent to Hermitian transposition of the
corresponding unitary matrices. In effect, the correlator corresponds to
matrix products in the order indicated in Fig. 7. The same expression with the
sextet wave-function is used for the $6\bar{6}$ correlator.

Although the diquarks on which the classification is based are pairs of $qq$
or $\bar{q}\bar{q}$, the delta functions connectors are between quarks and
antiquarks, like in mesonic case. Therefore, the resulting structures consists
of product of two traces (upper row) or a single trace (lower row) of Wilson
lines products. Like in the baryon case, these traces do not change if $W$ are
all rotated by some matrices $U$. Again there is no need to use $SU(3)$
averaging of $U$ (via Weingarten-style formulae).

Figure 7: The correlators of four Wilson lines in tetraquarks, connected
according to the product of color wave-functions (D.1) at times $\pm T$. The
arrow directions show quarks (1 and 2) and antiquarks (3,4). In the $\bar{3}3$
case the convolutions in the second row appears with the minus sign, and in
the $6\bar{6}$ case all terms have plus sign.

While four $W$ can be located arbitrarily, we have only calculated the
instanton contributions to correlators of four Wilson lines, placed at the
corners of a $r\times R$ rectangle. The quarks are set to be 1,2 at distance
$r$ from each other, and the distance between the diquark and the antidiquark
is set to be $R$. The potential, as usual, is defined as

$V_{\bar{3}3}(r,R)=n_{I+\bar{I}}\int
d^{3}y\bigg{(}1-K_{\bar{3}3}(r,R)\bigg{)}$ (25)

where, we recall that the calculation is basically the 3-dimensional
integration over the location of the instanton center $\vec{y}$. For the
$6\bar{6}$ case, the same formula applies, with the pertinent color wave-
function. In Fig. 8 we show only the results for the potential, from the
symmetric “square setting" of charges, $r=R=R_{9}/\sqrt{2}$, its dependence on
the overall size of the tetraquark, or the hyper-distance in 9 dimensions.

A key question is whether the $\bar{3}3$ and $6\bar{6}$ channels, support the
same effective radial potentials. Our results in Fig.8, show that $6\bar{6}$
and $\bar{3}3$ potentials are not identical. An inherent repulsion in the
$6\bar{6}$ channel generates an extra increase in the potential, especially at
large $R_{9}>1\,fm$. We conclude that the instanton-based forces approximately
support equality of the potentials in both channels, as anticipated in
Badalian  _et al._ (1987).

As in all similar calculations at small $R_{9}\rightarrow 0$, the strong
cancellations between the quark and antiquark contributions, lead to a
vanishing effective potential . The overall shape of the potential is similar
(but not equal) to other instanton-induced potentials (e.g. $qqq$ one in
Fig.5).

At large distances, the potentials tend to constants, since the instanton
vacuum does not confine. However, the approach to these constants is from
below, and does not happen even at $R_{9}=2\,fm$, the maximal value shown in
the figure.

Finally, we show by a straight line the non-perturbative part of the Cornell
potential used in the previous section. The overall force (potential gradient)
is quite close to the instanton-induced force.

Figure 8: Four-quark effective static potentials (in MeV) versus the
hyperdistance $R_{9}$ (in fm), calculated from the correlators of four Wilson
lines in tetraquarks. The color indices are connected according to products of
color wave functions at $T\rightarrow\pm\infty$(D.1) at times $\pm
T\rightarrow\pm\infty$. The results for the $\bar{3}3$ channel are shown by
triangles, and that for the $6\bar{6}$ channel by closed points. The line
corresponds to the non-perturbative part of the Cornell potential projected to
the hyper-distance $R_{9}$, as described in Appendix B.

Using the instanton-induced potentials into the radial Schroedinger equation
in the $\bar{3}3$ channel, we have for the s-level splittings

$\displaystyle M(2)-M(1)$ $\displaystyle\approx$ $\displaystyle 362\,MeV$ (26)
$\displaystyle M(3)-M(1)$ $\displaystyle\approx$ $\displaystyle 578\,MeV$
$\displaystyle M(4)-M(1)$ $\displaystyle\approx$ $\displaystyle 657\,MeV$

The comparison to the empirical splittings in (23), shows satisfactory
agreement. We emphasize that the the calculation of the radial potential does
not involve any free parameters, as we made use of the same instanton ensemble
for also the $\bar{q}q$ and $qqq$ channels.

##  VI Summary and discussion

We conclude, by first recalling some of the questions of hadronic spectroscopy
we have tried to address. A central one is the nature of the effective static
potential for quarkonia, including the perturbative contribution one-gluon-
exchange at short distance, and the non-perturbative and about linear
contribution at large distance. It has approximately the same tension in
quarkonia as deduced from Regge trajectories of light mesons and baryons.

In this paper we have addressed the effective interaction potentials in two
types of few-body hadrons – baryons and tetraquarks. While for the one-gluon-
exchange modifications by color factors can be easily calculated (resulting
e.g. in the known factor $\frac{1}{2}$ in baryons), the corresponding factors
in the non-perturbative part is quite nontrivial. It is even unclear whether
it can or cannot be dominated by a sum of binary interactions.

For simplicity and symmetry, we usde the so called “hyperdistance
approximation", projecting the relevant multidimensional potentials to a
single radial variable $R_{6}$ or $R_{9}$, in the space of Jacobi coordinates.
The squared hyper-distances are proportional to the sum of the squares of all
the binary distances, in a very symmetric way. Solving radial Schrodinger
equation one finds spectra and wave functions, in 6 and 9 dimensions
respectively. We focused on the splittings between the subsequent s-shell
levels, known to be rather insensitive to the quark masses, but roughly
proportional to the force between the constituents. We found that the Cornell-
type potential, averaged over solid angles in these spaces and with proper
color factors, does indeed generate reasonable values for $uuu=\Delta^{++}$
baryons, and for the recently discovered all-charm tetraquarks.

The static spin-independent potentials are known to be related to correlators
of few (3,4 etc) Wilson lines, connected in ways depending on the external
color wave functions of the given hadron. Such objects are very non-local
correlators of vacuum fields, with a non-trivial dependence on the locations
of the quarks. A priori it not possible to use a traditional notion of the
potential energy being a sum of binary potentials. For corresponding lattice
studies, see e.g. Koma and Koma (2017) for baryons and Bicudo  _et al._ (2017)
for tetraquarks.

The Euclidean gauge fields in the QCD vacuum can be, in a semiclassical
approximation, approximated by $instantons$, 4-dimensional solitons (pseudo-
particles) describing tunneling events through topological barriers which
separate configurations with different Chern-Simons numbers. Models of their
ensembles were developed, both using phenomenology of hadronic correlation
functions (for review see e.g.Schäfer and Shuryak (1998)) and direct lattice
studies of the fields. The instanton fields are special because they are
giving simple analytic answers for Wilson lines, and therefore one can derive
rather simple analytic expressions for their correlators Callan  _et al._
(1978). Basically, the Wilson lines are certain rotation in the quark color
spaces, by an angle that depends on the distance to the instanton center.

In Shuryak and Zahed (2023a) we introduced and used the “dense instanton
liquid model" to evaluate static and spin-dependent potentials in quarkonia.
In this paper we continue this work, by evaluating the instanton-induced
potentials via correlators of Wilson lines, for baryons and tetraquarks.

The instanton-induced potential for baryons is derived in section III.2.
Somewhat surprisingly, we found that the triple-W correlators can be written
as a sum of binary potentials. Furthermore, its shape happens to be the same
as for quarkonia, differing just by an overall factor. We think that the main
reason for this, is the fact that the instanton fields are always inside a
certain $SU(2)$ subspace of the $SU(3)$ color group, while the color wave
functions of baryons require that the three quarks carry always three mutually
orthogonal colors.

We have evaluated the instanton-induced potential for tetraquarks. In this
case the color wave functions are different. We have considerec two options,
which we referred to as $\bar{3}3$ and $6\bar{6}$ according to their diquark
color content. The convolution with four Wilson lines leads to certain
diagrams shown in Fig.7. Apparently in this case, no splitting of the 4-body
interactions into pairwise binary terms takes place. Yet the effective
potentials for the $\bar{3}3$ and $6\bar{6}$ channels as a function of the
hyperdistance $R_{9}$, shown in Fig.8, are similar in shape to those in
quarkonia and $qqq$ cases.

Before we made use of these potentials in spectroscopy, in section II.1 we
returned to the (rather well-studied) case of heavy quarkonia, and in
particular to the relation between their spectra and the effective potentials.
In order not to deal with the subtleties related to the quark masses and
additive terms in the potentials, we focused on the $splittings$ between the
s-shell levels, using the standard 3-dimensional Schroedinger equation. We
have seen that while the Cornell-type linear+Coulomb potentials do a very good
job, our instanton-induced potential can also be used, with a precision even
better than the fitted Martin potential.

Proceeding to the baryons and tetraquarks, we used a very similar approach,
based on the so called “hyper-spherical approximation", or “method of
K-harmonics" known since 1960’s. It is based on the assumption that the lowest
few-body states are “hyper-spherically symmetric", in the 6-dimensional space
for baryons, and 9-dimensional space for tetraquarks, in terms of pertinent
Jacobi coordinates. Solving the radial hyper-spherical Schroedinger equation,
we obtained the spectra, in particular the splittings between s-shell levels.
For baryons we did so for the $ccc$ and $uuu$ flavor-symmetric baryons. While
the formers are not yet experimentally observed, the latters belong to the
family of $\Delta^{++}$ for which we found a reasonable agreement between the
derived effective potential, and the experimentally observed separations of
the (s-shell) states.

In section IV.2, following the pioneering paper Badalian  _et al._ (1987), we
applied the hyper-spherical approximation to all-charm tetraquarks. This
method assumes that lowest the states are “hyper-spherical" in 9-dimensional
space of Jacobi coordinates. The interaction assumed to be the sum of 6 binary
potentials of the Cornell type times $\lambda_{i}^{A}\lambda_{j}^{A}$ color
matrices, times some extra effective strength factor. Focusing not on masses
themselves, but on the $splittings$ between the states, we simplified the
problem by eliminating the overall constants, a conspiracy between quark
masses and constants in the potentials. Our instanton-based potentials in
$\bar{3}3$ or $6\bar{6}$ channels are very close, so in a way we confirmed the
conjecture in Badalian  _et al._ (1987) that they are the same, in spite of
certain repulsive terms in the $6\bar{6}$ case. Furthermore, this potential
leads to $1S-2S-3S$ splittings (26) in fair agreement with the observed
splittings between the three recently reported resonances (23). We also
compared our results to the recent analysis in Badalian (2023), based on the
hyper-distance approximation as well.

Finally, we stress once again that the hyper-radial $qqq$ and
$cc\bar{c}\bar{c}$ effective potentials we derived, Figs.5 and 8, originate
without change from the instanton model of the QCD vacuum. Its two parameters
are fixed by lattice observation of instanton, and therefore one may say that
the present paper has $no$ free parameters whatsoever. The hadronic masses and
wave functions we derived, follow directly from a model of the QCD vacuum.

Acknowledgements:
This work is supported by the Office of Science, U.S. Department of Energy
under Contract No. DE-FG-88ER40388.

##  Appendix A Hyper-spherical harmonics

The non-relativistic kinetic energy of $A$ particles using $N$ Jacobi
coordinates, takes the simple form Fabre de la Ripelle and Navarro (1979)

$-\frac{1}{2}\sum_{i=1}^{A}\,\nabla^{2}_{i}=-\frac{1}{2}\bigg{(}\partial_{\rho}^{2}+\frac{3N-1}{\rho}\,\partial_{\rho}-\frac{1}{\rho^{2}}K_{N}^{2}\bigg{)}$
(27)

The hyper-spherical harmonics (HHs) are the eigenstates of the grand-angular
momentum

$K_{N}^{2}\,{\cal Y}_{[K]}^{KLM_{L}}(\Omega_{N})=(K(K+3N-2))\,{\cal
Y}_{[K]}^{KLM_{L}}(\Omega_{N})$ (28)

The $2^{N}+N$ angles
$\Omega_{N}=(\hat{\xi}_{1},..,\hat{\xi}_{N};\varphi_{1},..,\varphi_{N})$, are
those of the individual Jacobi coordinates $\hat{\xi}_{i}$, and the
hyperangles ${\rm cos}\varphi_{j}=\xi_{j}/\rho$. $L,M_{L}$ are the standard
quantum numbers of the total orbital angular momentum $L^{2},L_{z}$. The
specific form of the HHs follows by recoupling the individual angular momenta
$l_{i}$. They are normalized as

$\int d\Omega_{N}\,{\cal Y}_{[K]}^{KLM_{L}\ *}(\Omega_{N})\,{\cal
Y}_{[K^{\prime}]}^{K^{\prime}L^{\prime}M^{\prime}_{L}}(\Omega_{N})=\delta_{[K],[K^{\prime}]}$
(29)

and their total number is Fabre de la Ripelle and Navarro (1979)

$d_{K}=(2K+3N-2)\frac{(K+3N-3)!}{K!(3N-2)!}$ (30)

For instance for $A=4$ particles with $N=3$, the $K=0$ HH has degeneracy
$d_{0}=1$, and the $K=1$ HHs have degeneracy $d_{1}=9$.

##  Appendix B Quantum mechanics of three quarks in hyperspherical coordinates

To describe the specifics of the hyper-spherical approximation, we first
consider three quarks in a baryon. The three coordinate vectors of the quarks
are compressed into two (modified) Jacobi vectors

$\displaystyle\vec{\rho}=\frac{1}{\sqrt{2}}\vec{r}_{12}\qquad\vec{\lambda}=\frac{1}{\sqrt{6}}(\vec{r}_{13}+\vec{r}_{23})$
(31)

with $\vec{r}_{ij}=\vec{r}_{i}-\vec{r}_{j}$ as 3-dimensional vectors. The
6-dimensional radial coordinate denoted by $R_{6}$ is given by

$R_{6}^{2}=\vec{\rho}^{2}+\vec{\lambda}^{2}={1\over
3}[(\vec{r}_{1}-\vec{r}_{2})^{2}+(\vec{r}_{1}-\vec{r}_{3})^{2}+(\vec{r}_{3}-\vec{r}_{2})^{2}]$
(32)

so it is related to the sum of all the squared interquark distances, in a
symmetric manner. The coefficients in the definitions of these coordinates
also ensure that the kinetic energy takes the form of a Laplacian in
6-dimensional space $\\{\vec{\rho},\vec{\lambda}\\}$, divided by $2m$. (The
total momentum has a different coefficient, but it is presumed to be zero
anyway.)

The hyper-spherical approximation in its most symmetric form does not
distinguish $\vec{\rho},\vec{\lambda}$, but treats all 6 coordinates on equal
footing. Standard spherical coordinates can be introduced with 5 angles
$\theta_{1},\theta_{2},\theta_{3},\theta_{4}\in[0,\pi],\phi\in[0,2\pi]$, and a
6-dimensional solid angle

$\displaystyle\Omega_{6}$ $\displaystyle=$ $\displaystyle\int
d\theta_{1}d\theta_{2}d\theta_{3}d\theta_{4}d\phi$ $\displaystyle\times$
$\displaystyle
sin(\theta_{1})^{4}sin(\theta_{2})^{3}sin(\theta_{3})^{2}sin(\theta_{4})=\pi^{3}$

(Alternative a 6-dimensional coordinate space parametrization with five
angular variables can be defined as
$\vec{\rho}=R_{6}cos(\Phi)\hat{\vec{n}}_{1},\vec{\lambda}=R_{6}sin(\Phi)\hat{\vec{n}}_{2}$
with two unit vectors having independent solid angle integrals.)

The effective potential $V_{6}$ in general depends on all 6 coordinates, but
can be “projected onto the 6-dimensional radial coordinate $R_{6}$" by angular
average. If the binary potentials are of Cornell-type, a combination of linear
and inverse terms, e.g. $V_{12}=\kappa_{2}/r_{12}+\sigma_{2}r_{12}$, then
their projections are done via the solid angle integral

$\displaystyle\langle r_{12}\rangle$ $\displaystyle=$
$\displaystyle\sqrt{2}\langle|\vec{\rho}|\rangle=\sqrt{2}R_{6}{\int
d\Omega_{6}d_{12}\over\Omega_{6}}$ (34) $\displaystyle\approx$ $\displaystyle
0.960R_{6}$

where

$d_{12}=\sqrt{C_{1}^{2}+S_{1}^{2}C_{2}^{2}+S_{1}^{2}S_{2}^{2}C_{3}^{2}}=\sqrt{1-S_{1}S_{2}S_{3}}$
(35)

with the short-hand notations $C_{i}=cos(\theta_{i}),S_{i}=sin(\theta_{i})$.

A similar averaging of the inverse power of $r_{12}$ has $1/\sqrt{2}d_{12}$ in
the integral

$\displaystyle\bigg{<}{1\over r_{12}}\bigg{>}$ $\displaystyle=$
$\displaystyle{1\over\Omega_{6}}\int d\Omega_{6}{1\over\sqrt{2}d_{12}}$ (36)
$\displaystyle\approx$ $\displaystyle 1.20/R_{6}$

Note that the convergence of this integral is due to the fact that $d_{12}$
can only vanish if all three angles
$\theta_{1},\theta_{2},\theta_{3}\approx\pi/2$. As a result, the effective
potential in hyper-spherical coordinate $R_{6}$ takes the form

$V(R_{6})=3\big{[}0.96*R_{6}*(0.113/2)-1.2R_{6}*(0.64/2)\big{]}$ (37)

where the factor 3 appears due to the inclusion of interactions of $12,13,23$
quark pairs. The numerical values $0.113,0.64$ (in GeV2 units) are from our
version of the Cornell potential (3), and the division by two comes from
assumption of Ansatz A for the $qq$ relative to the $\bar{q}q$ potential.

##  Appendix C Quantum mechanics of four quarks in hyper-spherical coordinates

The generalized Jacobi coordinates for 4 particles are defined by

$\displaystyle\vec{\xi}_{1}$ $\displaystyle=$ $\displaystyle\sqrt{1\over
2}(\vec{r}_{1}-\vec{r}_{2})$ (38) $\displaystyle\vec{\xi}_{2}$
$\displaystyle=$ $\displaystyle\sqrt{1\over
6}(\vec{r}_{1}+\vec{r}_{2}-2\vec{r}_{3})$ $\displaystyle\vec{\xi}_{3}$
$\displaystyle=$ $\displaystyle{1\over
2\sqrt{3}}(\vec{r}_{1}+\vec{r}_{2}+\vec{r}_{3}-3\vec{r}_{4})$

In this case the hyper-distance

$\displaystyle R_{9}^{2}$ $\displaystyle=$
$\displaystyle\vec{\xi}_{1}^{2}+\vec{\xi}_{2}^{2}+\vec{\xi}_{3}^{2}$
$\displaystyle=$ $\displaystyle{1\over
4}[(\vec{r}_{1}-\vec{r}_{2})^{2}+(\vec{r}_{1}-\vec{r}_{3})^{2}+(\vec{r}_{1}-\vec{r}_{4})^{2}$
$\displaystyle+$
$\displaystyle(\vec{r}_{3}-\vec{r}_{2})^{2}+(\vec{r}_{4}-\vec{r}_{2})^{2}+(\vec{r}_{3}-\vec{r}_{4})^{2}]$

is connected to the sum of the squared distances of all six pairs of quarks.
When supplemented by 8 angles, it describes the 9-dimensional space in which
quantum mechanics is performed. Again, with the total momentum set to zero,
the kinetic energy is given by the 9-dimensional Laplacian with the standard
coefficient

$K=-{1\over 2m_{c}}\Delta_{\xi}$ (40)

The corresponding calculations are in section IV.2.

(To our surprise, in Badalian (2023) $non$-Jacobi coordinates are used, namely

$\displaystyle\vec{\xi}_{1}$ $\displaystyle=$ $\displaystyle\sqrt{1\over
2}(\vec{r}_{1}-\vec{r}_{2})$ (41) $\displaystyle\vec{\xi}_{2}$
$\displaystyle=$ $\displaystyle{1\over
2}(\vec{r}_{3}+\vec{r}_{4}-\vec{r}_{1}-\vec{r}_{2})$
$\displaystyle\vec{\xi}_{3}$ $\displaystyle=$ $\displaystyle\sqrt{1\over
2}(\vec{r}_{3}-\vec{r}_{4})$

for which neither the radial coordinate is the symmetric sum of all distances
as above, nor the kinetic energy is the Laplacian. )

The dependence on the hyper-distance only means that we deal with the lowest
$K=0$ or s-shell. Only the radial Schroedinger equation needs to be solved.
Note that after changing to the reduced wave-function, it differs from the
familiar 3-dimensional case, by only the “quasi-centrifugal" term
$12/R_{9}^{2}$.

The angular averaging is performed as in the previous section, except now
there is a different volume element, namely

$\displaystyle\Omega_{9}$ $\displaystyle=$
$\displaystyle\int(\prod_{i=1}^{7}d\theta_{i})d\phi$ $\displaystyle\times$
$\displaystyle
sin(\theta_{1})^{7}*sin(\theta_{2})^{6}*sin(\theta_{3})^{5}*sin(\theta_{4})^{4}$
$\displaystyle\times$ $\displaystyle
sin(\theta_{5})^{3}*sin(\theta_{6})^{2}*sin(\theta_{7})={32\pi^{4}\over 105}$

The angular integrations when averaging using the Cornell potential, are done
as in the previous section, with obvious changes

$\displaystyle\langle r_{12}\rangle$ $\displaystyle=$
$\displaystyle\sqrt{2}\langle|\vec{\rho}|\rangle={\sqrt{2}R_{9}\over\Omega_{9}}\int
d\Omega_{9}d_{12}$ (44) $\displaystyle\approx$ $\displaystyle 0.773R_{9}$
$\displaystyle\bigg{<}{1\over r_{12}}\bigg{>}$ $\displaystyle=$
$\displaystyle{1\over R_{9}\Omega_{9}}\int{d\Omega_{9}\over\sqrt{2}d_{12}}$
(45) $\displaystyle\approx$ $\displaystyle 1.55/R_{9}$

##  Appendix D All-charm tetraquarks

###  D.1 Color-spin wave functions

The natural way to define the wave function of the $\bar{c}\bar{c}cc$
tetraquark, is to start with $cc$ and $\bar{c}\bar{c}$ diquarks. Fermi
symmetry requires wave functions to be antisymmetric under quark permutations,
leaving two possibilities

$\displaystyle color$ $\displaystyle=$
$\displaystyle[\bar{3}]\,\,(antysym),\,\,\,spin=1(sym)$ (46) $\displaystyle
color$ $\displaystyle=$ $\displaystyle[6]\,\,(sym),\,\,\,spin=0\,\,(antisym)$
(47)

The total color is zero, and therefore the wave function should take the form

$\psi=F_{\bar{3}3}(\xi_{i})C_{\bar{3}3}(c_{i})S_{1}(s_{i})+F_{6\bar{6}}(\xi_{i})C_{6\bar{6}}(c_{i})S_{0}(s_{i})$
(48)

where the two spatial wave functions $F$ depend in general on 9 Jacobi
coordinates, $C(c_{i}),i=1..4$, on color indices, and $S(s_{i}),i=1..4$ quark
spin variables. If particles 1,2 are quarks and 3,4 antiquarks, one can define
the color and spin structures as

$\displaystyle C_{\bar{3}3}(c_{1}c_{2}c_{3}c_{4})$ $\displaystyle=$
$\displaystyle\sum_{c}\epsilon_{c_{1},c_{2},c}\epsilon_{c_{3},c_{4},c}/\sqrt{12}$
$\displaystyle=$
$\displaystyle\big{(}\delta_{c_{1}c_{3}}\delta_{c_{2}c_{4}}-\delta_{c_{1}c_{4}}\delta_{c_{2}c_{3}}\big{)}/\sqrt{12}$
$\displaystyle C_{6\bar{6}}(c_{1}c_{2}c_{3}c_{4})$ $\displaystyle=$
$\displaystyle\big{(}\delta_{c_{1}c_{3}}\delta_{c_{2}c_{4}}+\delta_{c_{1}c_{4}}\delta_{c_{2}c_{3}}\big{)}/\sqrt{24}$
$\displaystyle S_{0}$ $\displaystyle=$
$\displaystyle(\uparrow\downarrow-\downarrow\uparrow)/\sqrt{2}$ $\displaystyle
S_{1}$ $\displaystyle=$
$\displaystyle(\uparrow\downarrow+\downarrow\uparrow)/\sqrt{2}$

with the normalization

$\langle\psi|\psi\rangle=\int d^{9}\xi|F_{\bar{3}3,6\bar{6}}|^{2}$ (50)

Note that, while the color structures in (48) are not orthogonal to each
other, the spin structures are. So, neglecting the spin-dependent forces
(small for charm quarks), we expect the Hamiltonian to be diagonal, in terms
of $F_{\bar{3}3},F_{6\bar{6}}$.

Following Badalian  _et al._ (1987), we assume the color dependence of the
forces to be similar to those from one-gluon-exchange, namely

$V=\sum_{i>j}w_{ij}(r_{ij})\sum_{A=1,8}\frac{1}{4}\lambda^{A}_{i}\lambda^{A}_{j}$
(51)

The lower indices of Gell-Mann matrices indicate on which quark color index
they act.

An additional condition is that for antiquarks (3,4 in current setting), we
have to conjugate

$\lambda^{A}\rightarrow\tilde{\lambda}^{A}=-\lambda^{A*}$

whereby five contributions change sign and three do not. Averaging $V$ over
$\psi$ explicitly, and summing over all indices (we did it using Mathematica)
we obtained two expressions (which agree with those in Badalian  _et al._
(1987))

$\displaystyle\langle\psi|V|\rangle=\begin{bmatrix}-{1\over
3}F_{\bar{3}3}^{2}(w_{13}+w_{14}+w_{24}+w_{32}+2w_{12}+2w_{34}),\\\ -{1\over
6}F_{6\bar{6}}^{2}\big{(}-2w_{12}-2w_{34}+5w_{24}+5w_{13}+5w_{32}+5w_{14}\big{)}\end{bmatrix}$

Note that in the $6\bar{6}$ channel, there are repulsive quark-quark
contributions. Furthermore, $assuming$ that all $w_{ij}$ are the same, one
arrives at the conclusion in Badalian  _et al._ (1987), namely that under
these assumptions both color-spin tetraquark structures have the same
interactions. It is $twice$ stronger than in the $Q\bar{Q}$ channel (while
including all six $i,j$ pairs).

##  Appendix E Averaging four Wilson lines over the $SU(3)$ rotations

In this appendix we present expressions for 8 $SU(N_{c})$ matrices, averaged
over the invariant Haar measure. We do not use them in this paper, since both
for the $qqq$ and $qq\bar{q}\bar{q}$ cases considered it turned out that the
color wave functions are such that these complicated formulae are not needed.
Yet they would be needed for four Wilson lines with arbitrary colors. Since we
derived such formulae anyway, we decided to show them in this appendix.

###  E.1 Weingarten-style formula

One way to carry the color averaging over four pairs of unitary matrices can
be carried using the Weingarten method, producing the following expressions

$\displaystyle\bigg{<}U^{a_{1}}_{c_{1}}U^{\dagger
b_{1}}_{d_{1}}U^{a_{2}}_{c_{2}}U^{\dagger
b_{2}}_{d_{2}}U^{a_{3}}_{c_{3}}U^{\dagger
b_{3}}_{d_{3}}U^{a_{4}}_{c_{4}}U^{\dagger b_{4}}_{d_{4}}\bigg{>}_{U}=$
$\displaystyle+$
$\displaystyle\frac{6-8N^{2}_{c}+N_{c}^{4}}{N_{c}^{2}(N_{c}^{2}-1)(N_{c}^{2}-4)(N_{c}^{2}-9)}\sum_{n=1}^{4!}\delta^{a_{1}a_{2}a_{3}a_{4}}_{P_{n}(d_{1}d_{2}d_{3}d_{4})}\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{1}b_{2}b_{3}b_{4})}$
$\displaystyle+$
$\displaystyle\frac{4N_{c}-N_{c}^{3}}{N_{c}^{2}(N_{c}^{2}-1)(N_{c}^{2}-4)(N_{c}^{2}-9)}\sum_{n=1}^{4!}\delta^{a_{1}a_{2}a_{3}a_{4}}_{P_{n}(d_{1}d_{2}d_{3}d_{4})}$
$\displaystyle\times\bigg{(}\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{2}b_{1}b_{3}b_{4})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{1}b_{3}b_{2}b_{4})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{1}b_{2}b_{4}b_{3})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{4}b_{2}b_{3}b_{1})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{1}b_{4}b_{3}b_{2})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{3}b_{2}b_{1}b_{4})}\bigg{)}$
$\displaystyle+$
$\displaystyle\frac{6+N_{c}^{2}}{N_{c}^{2}(N_{c}^{2}-1)(N_{c}^{2}-4)(N_{c}^{2}-9)}\sum_{n=1}^{4!}\delta^{a_{1}a_{2}a_{3}a_{4}}_{P_{n}(d_{1}d_{2}d_{3}d_{4})}\bigg{(}\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{2}b_{1}b_{4}b_{3})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{3}b_{4}b_{1}b_{2})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{4}b_{3}b_{2}b_{1})}\bigg{)}$
$\displaystyle+$
$\displaystyle\frac{-3+2N_{c}^{2}}{N_{c}^{2}(N_{c}^{2}-1)(N_{c}^{2}-4)(N_{c}^{2}-9)}\sum_{n=1}^{4!}\delta^{a_{1}a_{2}a_{3}a_{4}}_{P_{n}(d_{1}d_{2}d_{3}d_{4})}$
$\displaystyle\times\bigg{(}\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{2}b_{3}b_{1}b_{4})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{1}b_{3}b_{4}b_{2})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{3}b_{2}b_{4}b_{1})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{2}b_{4}b_{3}b_{1})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{1}b_{4}b_{2}b_{3})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{3}b_{1}b_{2}b_{4})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{4}b_{1}b_{3}b_{2})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{4}b_{2}b_{1}b_{3})}\bigg{)}$
$\displaystyle+$
$\displaystyle\frac{-5N_{c}}{N_{c}^{2}(N_{c}^{2}-1)(N_{c}^{2}-4)(N_{c}^{2}-9)}\sum_{n=1}^{4!}\delta^{a_{1}a_{2}a_{3}a_{4}}_{P_{n}(d_{1}d_{2}d_{3}d_{4})}$
$\displaystyle\times\bigg{(}\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{2}b_{3}b_{4}b_{1})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{4}b_{1}b_{2}b_{3})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{2}b_{4}b_{1}b_{3})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{3}b_{1}b_{4}b_{2})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{3}b_{4}b_{2}b_{1})}+\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{4}b_{3}b_{1}b_{2})}\bigg{)}$

In the large $N_{c}$ limit, it simplifies considerably with the result

$\displaystyle\bigg{<}U^{a_{1}}_{c_{1}}U^{\dagger
b_{1}}_{d_{1}}U^{a_{2}}_{c_{2}}U^{\dagger
b_{2}}_{d_{2}}U^{a_{3}}_{c_{3}}U^{\dagger
b_{3}}_{d_{3}}U^{a_{4}}_{c_{4}}U^{\dagger
b_{4}}_{d_{4}}\bigg{>}_{U}=\frac{1}{N_{c}^{4}}\sum_{n=1}^{4!}\delta^{a_{1}a_{2}a_{3}a_{4}}_{P_{n}(d_{1}d_{2}d_{3}d_{4})}\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{1}b_{2}b_{3}b_{4})}+{\cal
O}\bigg{(}\frac{1}{N_{c}^{5}}\bigg{)}$ (53)

We are using the shorthand notation

$\displaystyle\sum_{n=1}^{4!}\delta^{a_{1}a_{2}a_{3}a_{4}}_{P_{n}(d_{1}d_{2}d_{3}d_{4})}\delta^{c_{1}c_{2}c_{3}c_{4}}_{P_{n}(b_{1}b_{2}b_{3}b_{4})}=\bigg{(}\delta^{a_{1}}_{d_{1}}\delta^{a_{2}}_{d_{2}}\delta^{a_{3}}_{d_{3}}\delta^{a_{4}}_{d_{4}}\delta^{c_{1}}_{b_{1}}\delta^{c_{2}}_{b_{2}}\delta^{c_{3}}_{b_{3}}\delta^{c_{4}}_{b_{4}}+{\rm
perm.}\bigg{)}$ (54)

where the sum is over the 4! $P_{n}$ elements of the permutation group
$S_{4}$, which are explicitly

$\displaystyle P_{1}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle
P_{2}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{1}a_{2}a_{4}a_{3})$ $\displaystyle
P_{3}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{1}a_{3}a_{2}a_{4})$ $\displaystyle
P_{4}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{1}a_{3}a_{4}a_{2})$ $\displaystyle
P_{5}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{1}a_{4}a_{2}a_{3})$ $\displaystyle
P_{6}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle
P_{7}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{1}a_{4}a_{3}a_{2})$ $\displaystyle
P_{8}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{2}a_{1}a_{3}a_{4})$ $\displaystyle
P_{9}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{2}a_{1}a_{4}a_{3})$ $\displaystyle
P_{10}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{2}a_{3}a_{1}a_{4})$ $\displaystyle
P_{11}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{2}a_{3}a_{4}a_{1})$ $\displaystyle
P_{12}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{2}a_{4}a_{1}a_{3})$ $\displaystyle
P_{13}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{2}a_{4}a_{3}a_{1})$ $\displaystyle
P_{14}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{3}a_{1}a_{2}a_{4})$ $\displaystyle
P_{15}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{3}a_{1}a_{4}a_{2})$ $\displaystyle
P_{16}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{3}a_{2}a_{4}a_{1})$ $\displaystyle
P_{17}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{3}a_{4}a_{1}a_{2})$ $\displaystyle
P_{18}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{3}a_{4}a_{2}a_{1})$ $\displaystyle
P_{19}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{4}a_{1}a_{2}a_{3})$ $\displaystyle
P_{20}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{4}a_{1}a_{3}a_{2})$ $\displaystyle
P_{21}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{4}a_{2}a_{1}a_{3})$ $\displaystyle
P_{22}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{4}a_{2}a_{3}a_{1})$ $\displaystyle
P_{23}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{4}a_{3}a_{1}a_{2})$ $\displaystyle
P_{24}(a_{1}a_{2}a_{3}a_{4})$ $\displaystyle=$
$\displaystyle(a_{4}a_{3}a_{2}a_{1})$ (55)

Note that coefficients have different powers of $1/N_{c}$, so when
$N_{c}\rightarrow\infty$ the only one left is the first raw. However at
$N_{c}=3$ they are all identical (by modulus) and have alternating signs, so
they conspire to cancel. At the other hand, coefficients have
$1/(N_{c}^{2}-9)$ factors prohibiting direct usage of this expressions when
$N_{c}=3$ (or smaller), without investigations of those cancellations.
Clearly, the result (E.1) is useful for numerical implementation for
$N_{c}>3$. For $N_{c}\leq p=4$ we now provide alternative and finite
identities.

###  E.2 Creutz formula

Another way to carry the color averaging in (E.1) which is pole free, is by
determinantal reduction. More specifically, we will use Creutz’s identities
Creutz (1978)

$\bigg{<}\prod_{i=1}^{N_{c}}U^{a_{i}}_{c_{i}}\bigg{>}_{U}=\frac{1}{N_{c}!}\epsilon^{a_{1}...a_{N_{c}}}\epsilon_{c_{1}...c_{N_{c}}}$
(56)

and

$\displaystyle U^{\dagger a}_{c}=$ $\displaystyle\frac{1}{(N_{c}-1)!}$
$\displaystyle\times\epsilon^{aa_{1}...a_{N_{c}-1}}\epsilon_{cc_{1}...c_{N_{c}-1}}U^{a_{1}}_{c_{1}}...U^{a_{N_{c}-1}}_{c_{{N_{c}-1}}}$

where $\epsilon^{a_{1}...a_{N_{c}}}$ is the Levi-Cevita tensor of
rank-$N_{c}$, with $\epsilon^{1...N_{c}}=1$.

With this in mind, we can substitute (E.2) for the $p=4$ string of
$U^{\dagger}$ in (E.1), with the result

$\displaystyle\bigg{<}U^{a_{1}}_{c_{1}}U^{\dagger
b_{1}}_{d_{1}}U^{a_{2}}_{c_{2}}U^{\dagger
b_{2}}_{d_{2}}U^{a_{3}}_{c_{3}}U^{\dagger
b_{3}}_{d_{3}}U^{a_{4}}_{c_{4}}U^{\dagger b_{4}}_{d_{4}}\bigg{>}_{U}=$
$\displaystyle\frac{1}{((N_{c}-1)!)^{p=4}}\prod^{p=4}_{m=1}\epsilon^{b_{m}b_{1m}...b_{(N_{c}-1)m}}\epsilon_{d_{m}d_{1m}..._{(N_{c}-1)m}}\bigg{<}\prod_{m=1}^{p=4}U^{a_{m}}_{c_{m}}U^{b_{1m}}_{d_{1m}}...U^{b_{(N_{c}-1)m}}_{d_{(N_{c}-1)m}}\bigg{>}_{U}$
(58)

Each of the product is composed rank-$N_{c}$ Levi-Cevita tensors. The last
unitary averaging in (E.2) can be undone by the determinantal identity (56)

$\displaystyle\bigg{<}U^{a_{1}}_{c_{1}}U^{\dagger
b_{1}}_{d_{1}}U^{a_{2}}_{c_{2}}U^{\dagger
b_{2}}_{d_{2}}U^{a_{3}}_{c_{3}}U^{\dagger
b_{3}}_{d_{3}}U^{a_{4}}_{c_{4}}U^{\dagger
b_{4}}_{d_{4}}\bigg{>}_{U}=\bigg{(}\frac{1}{((N_{c}-1)!)^{p}}\frac{2!...(N_{c}-1)!}{(p+1)!...(p+N_{c}-1)!}\bigg{)}_{p=4}$
$\displaystyle\prod^{p=4}_{m=1}\epsilon^{a_{m}b_{1m}...b_{(N_{c}-1)m}}\epsilon_{c_{m}d_{1m}...d_{(N_{c}-1)m}}\bigg{(}\prod_{m=1}^{p=4}\epsilon^{a_{m}b_{1m}...b_{(N_{c}-1)m}}\epsilon_{c_{m}d_{1m}...d_{(N_{c}-1)m}}+{\rm
perm.}\bigg{)}$ (59)

The combinatorical pre-factor follows from the fact that the group integration
is invariant under the permutation of the U-factors, with a total of
$(pN_{c})!$ permutations. Note that many permutations are identical, with only
$\frac{(pN_{c})!}{(p!(N_{c}!)^{p})}$ independent ones. Clearly all
coefficients in (E.2) are finite for any value of $N_{c}$, in contrast to
(E.1). We note that (E.2) holds for the averaging of even higher products
$(UU^{\dagger})^{p}$, and is free of poles whatever $p,N_{c}$.

##  Appendix F CNZ formula

Finally, another alternative formula for the color averaging can be obtained
by using the graphical color projection rules developed in CNZ paper
Chernyshev  _et al._ (1996). More specifically, the color averaging
$(UU^{\dagger})^{4}$ with $p=4$, can be tied recursively to $p=3$ by reduction

$\displaystyle\bigg{<}U^{a_{1}}_{c_{1}}U^{\dagger
b_{1}}_{d_{1}}U^{a_{2}}_{c_{2}}U^{\dagger
b_{2}}_{d_{2}}U^{a_{3}}_{c_{3}}U^{\dagger
b_{3}}_{d_{3}}U^{a_{4}}_{c_{4}}U^{\dagger b_{4}}_{d_{4}}\bigg{>}_{U}=$
$\displaystyle+\bigg{(}\bigg{[}\frac{1}{N_{c}}\delta^{a_{1}}_{d_{1}}\delta^{b_{1}}_{c_{1}}\,{\bf
1}_{1}\bigg{]}\bigg{<}U^{a_{2}}_{c_{2}}U^{\dagger
b_{2}}_{d_{2}}U^{a_{3}}_{c_{3}}U^{\dagger
b_{3}}_{d_{3}}U^{a_{4}}_{c_{4}}U^{\dagger b_{4}}_{d_{4}}\bigg{>}_{U}+{\rm
perm.}\bigg{)}$
$\displaystyle+\bigg{(}\big{(}[\lambda^{A}_{1}]^{a_{1}}_{d_{1}}[\lambda^{B}_{2}]^{a_{2}}_{d_{2}}[\lambda^{C}_{3}]^{a_{3}}_{d_{3}}[\lambda^{D}_{4}]^{a_{4}}_{d_{4}}\big{)}\,\big{(}[\lambda^{I}_{1}]^{b_{1}}_{c_{1}}[\lambda^{J}_{2}]^{b_{2}}_{c_{2}}[\lambda^{K}_{3}]^{b_{3}}_{c_{3}}[\lambda^{L}_{4}]^{b_{4}}_{c_{4}}\big{)}\,X^{ABCD}\,X^{IJKL}\bigg{)}$
(60)

with the color factors

$\displaystyle X^{ABCD}=$
$\displaystyle\frac{1}{16{(N_{c}^{2}-1)^{2}}}(\delta^{AB}\delta^{CD}+\delta^{AC}\delta^{BD}+\delta^{AD}\delta^{BC})$
(61)
$\displaystyle+\frac{1}{4(N_{c}^{2}-1)}\frac{N_{c}^{2}}{4{(N_{c}^{2}-4)^{2}}}\big{(}d^{ABM}d^{CDM}+d^{ACM}d^{BDM}+d^{ADM}d^{BCM}\big{)}$
$\displaystyle+\frac{1}{4(N_{c}^{2}-1)}\frac{1}{4N_{c}^{2}}\big{(}f^{ABM}f^{CDM}+f^{ACM}f^{BDM}+f^{ADM}f^{BCM}\big{)}$
$\displaystyle+\frac{1}{4(N_{c}^{2}-1)}\frac{1}{4{(N_{c}^{2}-4)}}\big{(}d^{ABM}f^{CDM}+d^{ACM}f^{BDM}+d^{ADM}f^{BCM}\big{)}$
$\displaystyle+\frac{1}{4(N_{c}^{2}-1)}\frac{1}{4{(N_{c}^{2}-4)}}\big{(}f^{ABM}d^{CDM}+f^{ACM}d^{BDM}+f^{ADM}d^{BCM}\big{)}$

and similarly for $X^{IJKL}$. The reduced averaging over $(UU^{\dagger})^{3}$
gives

$\displaystyle\bigg{<}U^{a_{2}}_{c_{2}}U^{\dagger
b_{2}}_{d_{2}}U^{a_{3}}_{c_{3}}U^{\dagger
b_{3}}_{d_{3}}U^{a_{4}}_{c_{4}}U^{\dagger
b_{4}}_{d_{4}}\bigg{>}_{U}=\bigg{[}\frac{1}{N_{c}}\delta^{a_{2}}_{d_{2}}\delta^{b_{2}}_{c_{2}}\,{\bf
1}_{2}\bigg{]}\bigg{[}\frac{1}{N_{c}}\delta^{a_{3}}_{d_{3}}\delta^{b_{3}}_{c_{3}}\,{\bf
1}_{3}\bigg{]}\bigg{[}\frac{1}{N_{c}}\delta^{a_{4}}_{d_{4}}\delta^{b_{4}}_{c_{4}}\,{\bf
1}_{4}\bigg{]}$
$\displaystyle+\bigg{(}\bigg{[}\frac{1}{N_{c}}\delta^{a_{2}}_{d_{2}}\delta^{b_{2}}_{c_{2}}\,{\bf
1}_{2}\bigg{]}\bigg{[}\frac{1}{4(N_{c}^{2}-1)}[\lambda_{3}^{A}]^{a_{3}}_{d_{3}}[\lambda_{4}^{A}]^{a_{4}}_{d_{4}}[\lambda_{3}^{B}]^{b_{3}}_{c_{3}}[\lambda_{4}^{B}]^{b_{4}}_{c_{4}}\bigg{]}+2\,\,{\rm
perm.}\bigg{)}$
$\displaystyle+\frac{N_{c}}{8(N_{c}^{2}-1)(N_{c}^{2}-4)}\bigg{[}d^{ABC}[\lambda_{2}^{A}]^{a_{2}}_{d_{2}}[\lambda_{3}^{B}]^{a_{3}}_{d_{3}}[\lambda_{4}^{C}]^{a_{4}}_{d_{4}}\bigg{]}\bigg{[}d^{IJK}[\lambda_{2}^{I}]^{b_{2}}_{c_{2}}[\lambda_{3}^{J}]^{b_{3}}_{c_{3}}[\lambda_{4}^{K}]^{b_{4}}_{c_{4}}\bigg{]}$
$\displaystyle+\frac{1}{8N_{c}(N_{c}^{2}-1)}\bigg{[}f^{ABC}[\lambda_{2}^{A}]^{a_{2}}_{d_{2}}[\lambda_{3}^{B}]^{a_{3}}_{d_{3}}[\lambda_{4}^{C}]^{a_{4}}_{d_{4}}\bigg{]}\bigg{[}f^{IJK}[\lambda_{2}^{I}]^{b_{2}}_{c_{2}}[\lambda_{3}^{J}]^{b_{3}}_{c_{3}}[\lambda_{4}^{K}]^{b_{4}}_{c_{4}}\bigg{]}$
(62)

## References

  * Shuryak and Zahed (2023a) E. Shuryak and I. Zahed, Phys. Rev. D 107, 034023 (2023a), arXiv:2110.15927 [hep-ph] .
  * Simonov (1966) Y. A. Simonov, Sov. J. Nucl. Phys. 3, 461 (1966).
  * Shuryak and Torres-Rincon (2019) E. Shuryak and J. M. Torres-Rincon, Phys. Rev. C 100, 024903 (2019), arXiv:1805.04444 [hep-ph] .
  * Badalian  _et al._ (1987) A. M. Badalian, B. L. Ioffe, and A. V. Smilga, Nucl. Phys. B 281, 85 (1987).
  * Chernyshev  _et al._ (1996) S. Chernyshev, M. A. Nowak, and I. Zahed, Phys. Rev. D 53, 5176 (1996), arXiv:hep-ph/9510326 .
  * Ferraris  _et al._ (1995) M. Ferraris, M. M. Giannini, M. Pizzo, E. Santopinto, and L. Tiator, Phys. Lett. B 364, 231 (1995).
  * Gandhi  _et al._ (2018) K. Gandhi, Z. Shah, and A. K. Rai, Eur. Phys. J. Plus 133, 512 (2018), arXiv:1811.00251 [hep-ph] .
  * Brambilla  _et al._ (2022) N. Brambilla _et al._ , (2022), arXiv:2203.16583 [hep-ph] .
  * Koma and Koma (2017) Y. Koma and M. Koma, Phys. Rev. D 95, 094513 (2017), arXiv:1703.06247 [hep-lat] .
  * Bicudo  _et al._ (2017) P. Bicudo, M. Cardoso, O. Oliveira, and P. J. Silva, Phys. Rev. D 96, 074508 (2017), arXiv:1702.07789 [hep-lat] .
  * Shuryak (1982) E. V. Shuryak, Nucl. Phys. B 203, 93 (1982).
  * Schäfer and Shuryak (1998) T. Schäfer and E. V. Shuryak, Rev. Mod. Phys. 70, 323 (1998), arXiv:hep-ph/9610451 .
  * Workman  _et al._ (2022) R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Grosse and Martin (1980) H. Grosse and A. Martin, Phys. Rept. 60, 341 (1980).
  * Callan  _et al._ (1978) C. G. Callan, Jr., R. F. Dashen, D. J. Gross, F. Wilczek, and A. Zee, Phys. Rev. D 18, 4684 (1978).
  * Shuryak and Zahed (2023b) E. Shuryak and I. Zahed, Phys. Rev. D 107, 034024 (2023b), arXiv:2111.01775 [hep-ph] .
  * Shuryak and Zahed (2021) E. Shuryak and I. Zahed, (2021), arXiv:2112.15586 [hep-ph] .
  * Guimarães  _et al._ (1981) A. B. Guimarães, H. T. Coelho, and R. Chanda, Phys. Rev. D 24, 1343 (1981).
  * Aaij  _et al._ (2020) R. Aaij _et al._ (LHCb), Sci. Bull. 65, 1983 (2020), arXiv:2006.16957 [hep-ex] .
  * Zhang and Yi (2022) J. Zhang and K. Yi (CMS), PoS ICHEP2022, 775 (2022), arXiv:2212.00504 [hep-ex] .
  * Xu (2023) Y. Xu (ATLAS), Acta Phys. Polon. Supp. 16, 21 (2023), arXiv:2209.12173 [hep-ex] .
  * Karliner and Rosner (2020) M. Karliner and J. L. Rosner, Phys. Rev. D 102, 114039 (2020), arXiv:2009.04429 [hep-ph] .
  * Giron and Lebed (2020) J. F. Giron and R. F. Lebed, Phys. Rev. D 102, 074003 (2020), arXiv:2008.01631 [hep-ph] .
  * Badalian (2023) A. M. Badalian, (2023), arXiv:2305.04585 [hep-ph] .
  * Fabre de la Ripelle and Navarro (1979) M. Fabre de la Ripelle and J. Navarro, Annals of Physics 123, 185 (1979).
  * Creutz (1978) M. Creutz, J. Math. Phys. 19, 2043 (1978).


* [10] Y. Ameur, H. Hedenmalm and N. Makarov, _Fluctuations of eigenvalues of random normal matrices_ , Duke Math. J. 159 (2011), 31–81.
  * [11] S. Belinschi, M. A. Nowak, R. Speicher and W. Tarnowski, _Squared eigenvalue condition numbers and eigenvector correlations from the single ring theorem_ , J. Phys. A 50 (2017), 105204.
  * [12] F. Benaych-Georges and F. Chapon, _Random right eigenvalues of Gaussian quaternionic matrices_ , Random Matrices Theory Appl. 1 (2012), 1150009.
  * [13] F. Benaych-Georges and O. Zeitouni, _Eigenvectors of non normal random matrices_ , Electron. Commun. Probab. 23 (2018), 1–12.
  * [14] P. Bourgade and G. Dubach, _The distribution of overlaps between eigenvectors of Ginibre matrices_ , Probab. Theory Relat. Fields 177 (2020), 397–464.
  * [15] P. Bourgade and H.-T. Yau, _The eigenvector moment flow and local quantum unique ergodicity_ , Comm. Math. Phys. 350 (2017), 231–278.
  * [16] Z. Burda, B.J. Spisak and P. Vivo, Eigenvector statistics of the product of Ginibre matrices, Phys. Rev. E 95 (2017), 022134.
  * [17] S.-S. Byun and C. Charlier, _On the almost-circular symplectic induced Ginibre ensemble_ , Stud. Appl. Math. 150 (2023), 184–217.
  * [18] S.-S. Byun and M. Ebke, _Universal scaling limits of the symplectic elliptic Ginibre ensembles_ , Random Matrices Theory Appl. 12 (2023), 2250047.
  * [19] S.-S. Byun, M. Ebke and S.-M. Seo, _Wronskian structures of planar symplectic ensembles_ , Nonlinearity 36 (2023), 809–844.
  * [20] S.-S. Byun and P. J. Forrester, _Spherical induced ensembles with symplectic symmetry_ , SIGMA Symmetry Integrability Geom. Methods Appl. 19 (2023), 033, 28pp.
  * [21] S.-S. Byun and P. J. Forrester, _Progress on the study of the Ginibre ensembles I: GinUE_ , arXiv:2211.16223.
  * [22] S.-S. Byun and P. J. Forrester, _Progress on the study of the Ginibre ensembles II: GinOE and GinSE_ , arXiv:2301.05022.
  * [23] S.-S. Byun, S.-Y. Lee and M. Yang, _Lemniscate ensembles with spectral singularity_ , arXiv:2107.07221.
  * [24] S.-S. Byun and K. Noda, _Scaling limits of complex and symplectic non-Hermitian Wishart ensembles_ , arXiv:2402.18257.
  * [25] J. T. Chalker and B. Mehlig, _Eigenvector statistics in non-Hermitian random matrix ensembles_ , Phys. Rev. Lett. 81 (1998), 3367–3370.
  * [26] J. T. Chalker and B. Mehlig, _Statistical properties of eigenvectors in non-Hermitian Gaussian random matrix ensembles_ , J. Math. Phys. 41 (2000), 3233–3256.
  * [27] C. Charlier, _Asymptotics of determinants with a rotation-invariant weight and discontinuities along circles_ , Adv. Math. 408 (2022), 108600.
  * [28] C. Cipolloni and D. Schröder, _On the condition number of the shifted real Ginibre ensemble_ , SIAM J. Matrix Anal. Appl. 43 (2022), 1469–1487.
  * [29] C. Cipolloni, L. Erdös, J. Henheik and D. Schröder, _Optimal lower bound on eigenvector overlaps for non-Hermitian random matrices_ , J. Funct. Anal. 287 (2024), 110495.
  * [30] N. Crawford and R. Rosenthal, _Eigenvector correlations in the complex Ginibre ensemble_ , Ann. Appl. Probab. 32 (2022), 2706–2754.
  * [31] M. J. Crumpton, Y. V. Fyodorov and T. R. Würfel, _Mean eigenvector self-overlap in the real and complex elliptic Ginibre ensembles at strong and weak non-Hermiticity_ , arXiv:2402.09296.
  * [32] M. Crumpton, T. R. Würfel, _Spectral density of complex eigenvalues and associated mean eigenvector self-overlaps at the edge of elliptic Ginibre ensembles_ , arXiv:2405.02103.
  * [33] D. S. Dean, P. Le Doussal, S. N. Majumdar and G. Schehr, _Noninteracting fermions in a trap and random matrix theory_ , J. Phys. A 52 (2019), 144006.
  * [34] G. Dubach, On eigenvector statistics in the spherical and truncated unitary ensembles, Electron. J. Probab. 26 (2021), Paper No. 124, 29 pp.
  * [35] G. Dubach, _Symmetries of the quaternionic Ginibre ensemble_ , Random Matrices Theory Appl. 10 (2021), 2150013.
  * [36] G. Dubach, _Explicit formulas concerning eigenvectors of weakly non-unitary matrices_ , Electron. Commun. Probab. 28 (2023), Paper No. 6, 11 pp.
  * [37] L. Erdös and H. C. Ji, _Wegner estimate and upper bound on the eigenvalue condition number of non-Hermitian random matrices_ , Comm. Pure Appl. Math. (to appear), arXiv:2301.04981.
  * [38] S. Esaki, M. Katori and S. Yabuoku, _Eigenvalues, eigenvector-overlaps, and regularized Fuglede-Kadison determinant of the non-Hermitian matrix-valued Brownian motion_ , arXiv:2306.00300.
  * [39] P. J. Forrester, _Log-gases and random matrices_ , Princeton University Press, Princeton, NJ, 2010.
  * [40] P. J. Forrester, _Analogies between random matrix ensembles and the one-component plasma in two-dimensions_ , Nuclear Phys. B 904 (2016), 253–281.
  * [41] Y. V. Fyodorov and B. Mehlig, _Statistics of resonances and nonorthogonal eigenfunctions in a model for single-channel chaotic scattering_ , Phys. Rev. E. 66 (2002), 045202.
  * [42] Y. V. Fyodorov and D. V. Savin, _Statistics of resonance width shifts as a signature of eigenfunction non-orthogonality_ , Phys. Rev. Lett. 108 (2012), 184101.
  * [43] Y. V. Fyodorov, _On statistics of bi-orthogonal eigenvectors in real and complex Ginibre ensembles: combining partial Schur decomposition with supersymmetry_ , Comm. Math. Phys. 363 (2018), 579–603.
  * [44] Y. V. Fyodorov and W. Tarnowski, _Condition numbers for real eigenvalues in the real elliptic Gaussian ensemble_ , Ann. Henri Poincaré 22, (2021), 309–330.
  * [45] Y. V. Fyodorov and H.-J. Sommers, _Random matrices close to hermitian or unitary: overview of methods and results_ , J. Phys. A 36 (2003), 3303–3347.
  * [46] H. Hedenmalm and A. Wennman, _Planar orthogonal polynomials and boundary universality in the random normal matrix model_ , Acta Math. 227 (2021), 309–406.
  * [47] J. Grela and P. Warchol, _Full Dysonian dynamics of the complex Ginibre ensemble_ , J. Phys. A 51 (2018), 425203.
  * [48] J.-B. Gros, U. Kuhl, O. Legrand, F. Mortessagne, E. Richalot and D.V. Savin, _Experimental width shift distribution: A test of nonorthogonality for local and global perturbations_ , Phys. Rev. Lett. 113 (2014), 224101.
  * [49] M. B. Hastings, _Fermionic mapping for eigenvalue correlation functions of weakly non-Hermitian symplectic ensemble_ , Nucl. Phys. B 572 (2000), 535–546.
  * [50] M. E. H. Ismail and D. R. Masson, _Generalized orthogonality and continued fractions_ , J. Approx. Theory 83 (1995), 1–40.
  * [51] R. A. Janik, W. Nörenberg, M. A. Nowak, G. Papp and I. Zahed, _Correlations of eigenvectors for non Hermitian random matrix models_ , Phys. Rev. E. 60 (1999), 2699–2705.
  * [52] E. Kanzieper, _Eigenvalue correlations in non-Hermitean symplectic random matrices_ , J. Phys. A 35 (2002), 6631–6644.
  * [53] J. S. Kim and D. Stanton, _Combinatorics of orthogonal polynomials of type $R_{I}$_, Ramanujan J. 61 (2023), 329–390.
  * [54] A. Knowles and J. Yin, _Eigenvector distribution of Wigner matrices_ , Probab. Theory Related Fields, 155 (2013), 543–582.
  * [55] A. V. Kolesnikov and K. B. Efetov, _Distribution of complex eigenvalues for symplectic ensembles of non-Hermitian matrices_ , Waves in Random Media 9 (1999), 71–82
  * [56] S.-Y. Lee and M. Yang, _Strong asymptotics of planar orthogonal polynomials: Gaussian weight perturbed by finite number of point charges_ , Comm. Pure Appl. Math. 76 (2023), 2888–2956.
  * [57] T. A. Loring, Factorization of matrices of quaternions, Expo. Math. 30 (2012), 250.
  * [58] B. Mehlig and M. Santer, _Universal eigenvector statistics in a quantum scattering ensemble_ , Phys. Rev. E. (2001) 63, 020105.
  * [59] M. L. Mehta, _Random Matrices_ , Academic Press, New York, second edition, 1991.
  * [60] K. Noda, _Determinantal structure of the overlaps of the induced Ginibre unitary ensemble_ , arXiv:2310.15362.
  * [61] K. Noda, _Determinantal structure of the overlaps for induced spherical unitary ensemble_ , arXiv:2312.12690.
  * [62] M. A. Nowak and W. Tarnowski, _Probing non-orthogonality of eigenvectors in non-Hermitian matrix models: diagrammatic approach_ , J. High Energ. Phys. 2018 (2018), 152.
  * [63] F. W. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark (Editors). NIST Handbook of Mathematical Functions. Cambridge University Press, Cambridge, 2010.
  * [64] S. O’Rourke, V. Vu and K. Wang, _Eigenvectors of random matrices: a survey_ , J. Combin. Theory. Ser. A 144 (2016), 361–442.
  * [65] A. Serebryakov and N. Simm, _Schur function expansion in non-Hermitian ensembles and averages of characteristic polynomials_ , arXiv:2310.20686.
  * [66] S. Serfaty, _Lectures on Coulomb and Riesz Gases_ , 2024.
  * [67] T. Shirai and Y. Takahashi, _Random point fields associated with certain Fredholm determinants I: fermion, Poisson and boson point processes_ , J. Funct. Anal. 205 (2003), 414–463.
  * [68] T. Tao and V. Vu, _Random matrices: universal properties of eigenvectors_ , Random Matrices Theory Appl. 1 (2012), 1150001.
  * [69] W. Tarnowski, _Condition numbers for real eigenvalues of real elliptic ensemble: weak non-normality at the edge_ , J. Phys. A 57 (2024), 255204.
  * [70] M. Walters and S. Starr, _A note on mixed matrix moments for the complex Ginibre ensemble_ , J. Math. Phys. 56 (2015), 013301.
  * [71] T. R. Würfel, M. J. Crumpton and Y. V. Fyodorov, _Mean left-right eigenvector self-overlap in the real Ginibre ensemble_ , Random Matrices Theory Appl. (to appear), https://doi.org/10.1142/S2010326324500175, arXiv:2310.04307.
  * [72] S. Yabuoku, _Eigenvalue processes of elliptic Ginibre ensemble and their overlaps_ , Int. J. Math. Ind. 12 (2020), 2050003, 17 pp.


11institutetext: Chungnam National University, Daejeon, South Korea
22institutetext: The Affiliated Institute of ETRI, Daejeon, South Korea
22email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# SSLEM: A Simplifier for MBA Expressions based on Semi-linear MBA Expressions
and Program Synthesis

Seong-Kyun Mok 11 Seoyeon Kang 11 Jeongwoo Kim 11 Eun-Sun Cho 11 Seokwoo Choi
22

###### Abstract

MBA (mixed boolean and arithmetic) expressions are hard to simplify, so used
for malware obfuscation to hinder analysts’ diagnosis. Some MBA simplification
methods with high performance have been developed, but they narrowed the
target to ”linear” MBA expressions, which allows efficient solutions based on
logic/term-rewriting. However such restrictions are not appropriate for
general forms of MBA expressions usually appearing in malware. To overcome
this limitation, we introduce a ”semi-linear” MBA expression, a new class of
MBA expression extended from a linear MBA expression, and propose a new MBA
simplifier called ”SSLEM”, based on a simplification idea of semi-linear MBA
expressions and program synthesis.

###### Keywords:

MBA obfuscation logic-based simplifiers program synthesis-based simplifiers

##  1 Introduction

MBA (Mixed Boolean and Arithmetic) expressions are the expressions mixed with
arithemethic operators like ‘$+$’, ‘$-$’ and $\times$’, and boolean operators
like ‘$\&$’, ‘$|$’ and ‘$\wedge$’. Program rewriting to insert MBA expressions
is a popular way of program obfuscation because such expressions are difficult
to analyze with popular solvers like Z3 and Mathematica.

While MBA obfuscation is a helpful method for software protection, it is also
used for malware obfuscation. Once malware is obfuscated with the MBA-
obfuscation option, it entails a herculean effort of analysts to understand
its behavior. Thus there have been various simplification efforts to handle
this problem.

The methods used by existing works to simplify MBA expressions are classified
into two kinds; First, logic-based simplifiers have been devised to simplify
MBA expressions. Based on theoretical foundation [1][2], this logic-based
approach suffers from the scalability problem [3]. Thus most of the previous
works narrowed the problem scope. For instance, typically logic-based
simplifiers have been dedicated to a specific form of MBA expressions like
linear MBA expressions to make the problem tractable, systematic and fast [2].
However, it does not seem appropriate for obfuscated malware, because trivial
MBA expressions shown in obfuscated malware like $x\times y$, $x\ \&\ (x+y)$
and $x\ \&\ 0x1$ are not linear expressions.

Second, synthesis-based MBA expression simplication [4][5][6] has been
recently introduced. This approach searches a simpler expression of the
original program, if any, by program synthesis. It is applicable to general
complicated forms of MBA expressions, besides specific forms of expressions.
However, the solution is imcomplete and relatively slow [6].

In this poster, we proposed a method and tool to develop an enhanced
simplifier of MBA expressions to help deobfuscate MBA-obfuscated malware,
taking advantage of program synthesis as well as logic-based solutions.

##  2 Previous Works

Previous works recognized nested categories of MBA expressions. “Linear MBA
expressions” is a summation of terms, where a summation means an arithmetic
operation with $+$, a term is a multiplication of a coefficient $a_{i}$ and a
pure boolean expression $e_{i}(x_{1},\dots x_{T})$, $T$ is the number of
variables and $N$ is the number of terms, as depicted in the following fomula
[1][2][3];

$\sum_{i=0}^{N}a_{i}e_{i}(x_{1},\dots x_{T})$ (1)

Examples of linear expressions include $10\times(x\ \&\ y)$, $5\times(x\ |\
y)+6\times(x\ \&\ y)$, $3\times(x\ \&\ y)-20\times(x\ |\ y)-1$ and
$y+376\times(x\ \&\ -1)$.

As mentioned earlier, most traditional simplification of MBA expressions
dedicated to linear MBA expressions. Depending on pattern matching and rule
rewritting, they are not feasible for complicated non-linear expressions due
to the lack of the scalability.

It is proved a valuable proposition that a simplified form of any of the
linear MBA expressions with 1-bit variables also represents a simplified form
of the one with the n-bit variables, and vice versa [2][1]. This means that if
we simplify any linear MBA expression with 1-bit variables, the simplification
result can be also considered as a simplification of the corresonding, same
looking expression with n-bit variables [2]. This makes the simplification
process much more efficient.

For some insight, let us consider an expression a function which maps a
combination of input variables to the output. Regardless representation of an
expression whether it is complicated or simple, the number of unique mappings
will be limited by the size and the number of the input and output variables.
For instance, if we use two 1-bit variables $x$ and $y$ as inputs, the
possible combination of $x$ and $y$ will be 00, 01, 10 and 11. The
simplification of a given linear MBA expression would be a searching problem.
Among all the possible output-bit-combinations for each input-bit-combination,
a simplifier should find a simpler expression corresponding to the given
expression. The property proven by the previous works dramatically reduced the
searching space from $2^{2^{2^{16}}}$ (that is, $16^{16}=1.84e+19$) to
$2^{2^{2^{1}}}$ (that is, $8$) .

Note that the proposition proved by previous methods only for linear MBA
expressions. However, solution dedicated to linear MBA expressions cannot
handle a lot of general MBA expressions-only $2^{2^{2^{1}}}$ expressions can
be handled. Thus, they adopted heuristics to spot any linear subexpressions in
a given complicated expression. Otherwise $2^{2^{2^{1}6}}-2^{2^{2^{1}}}$ MBA
expressions cannot be handled, some of which might appear in obfuscated
malware.

##  3 SSLEM: The Proposed MBA Expression Simplifier

In this section, we introduce SSLEM, a new MBA expression simplification
solution To overcome limitations of existing works, SSLEM makes use of both
kinds of existing tools (the logic-based simplifier and the program synthesis-
based simplifier), we envision that the proposed mehod will enhance
performance some corner cases.

As mentioned earlier, a logic-based simplifier is fast and effiecient but
simplifying only limited category of expressions (linear MBA expressions).
Aiming to overcome this problem, the proposed SSLEM proposes to preprocess the
general MBA expressions to linear MBA expressions, and input the result to a
logic-based simplifier. However, transformation to linear expressions would
fail in most cases, because linear MBA expressions are only a small subset of
the whole MBA expressions, as we compared the number of unique semantics in
the previous section.

Thus, SSLEM uses an intermediate form of MBA expressions called “semi-linear”
MBA expressions, and transforms a general MBA expressions to a semi-linear MBA
expression, first, by the program synthesis-based simplifier. Then the
resulting semi-linear MBA expression is transformed into a combined expression
which is made of several linear expressions. Since each of the resulting
linear MBA expressions can be simplified by the logic-based simplifier, we can
finally earn simplification of a non-linear expression with a logic-based
simplier.

Figure 1:  Overview of the proposed method and existing method

Overall process of the proposed SSLEM is composed of three steps, as depicted
in Fig. 1:

  1. 1.

PLASynth-SL finds preliminary candidates for the logic-based simplifier
working at Step III. It is a search-based synthesizer to find a semi-linear
MBA expression corresponding to a general non-linear MBA expressions, keeping
the semantics the same.

  2. 2.

LEM transforms a semi-linear MBA expression to a combination of linear MBA
expressions

  3. 3.

The linear subexpressions of the semi-linear MBA expression resulting from
LEM, are in turns fed to the logic-based simplifier. This step might leverage
an existing powerful linear MBA expression simplifier like MBA Blast [2]

###  3.1 Semi-linear MBA expressions and PLASynth-SL

Let us start with a “heterogeneous-bit-valued constant”, a constant that has
different bit values in it (eg. 0xFE, 0x1, 0x003F), Note that a linear
expression only a “homogeneous-bit-valued constant” in its boolean
subexpressions, which includes 0xFF, -1 and 0x0, each of which consist of the
same valued bits (all 1’s or all 0’s). Semi-linear MBA expressions are in a
slightly extended form of linear expressions, defined as follows;

###### Definition 1

A semi-linear MBA expression is an MBA Expression that could be a linear MBA
expression, if without “heterogeneous-bit-valued constants” in the boolean
subexpressions. Semi-linear expressions inlucde ”$30\times(0xFE\ |\ y)$”,
“$21\times x\ \&\ 0x01)$”

PLASynth-SL is equipped with optimization and timeout management [7] for
better performance. In a preliminary experiment, out of the 500 MBA expression
samples obfuscated by Tigress [8] with the MBA-obfuscation option, PLASynth-SL
earned 68 semi-linear MBA expressions, while gained 46 linear MBA expressions.

###  3.2 LEM (Linear MBA Expression Generator from Semi-linear MBA
Expressions)

Let us assume that a given 16-bit semi-linear MBA expression $f$ is “$x\wedge
0x003F$”, Here, “0x003F” is the heterogenous-bit-valued constant, without it,
$f$ could be a linear MBA expression. As depicted in Fig. 2, this constant is
divided into 0000 0000 00 and 11 1111.

The MBA expression $f$ is considered as a combination of $f_{1}$ and $f_{2}$,
where $x$ = [$X_{1}^{10},X_{0}^{6}$] in the multi-granularity representation
[9] $f_{1}$ = $X_{1}^{10}\wedge 0x00000$, and $f_{2}$ = $X_{0}^{6}\wedge
0x3F$. At this time, each of $f_{1}$ and $f_{2}$ is a linear MBA expression.
This means that the simplified expression of $f_{1}$, say $F_{1}$, is
syntatically same to $F_{1}^{{}^{\prime}}$, where $F_{1}^{{}^{\prime}}$ is the
simplified expression of $f_{1}$ when interpreted as a 1-bit expression
[1][2][3]. (Similarly, $F_{2}$, is syntatically same to
$F_{2}^{{}^{\prime}}$.) as a 1-bit expression. Thus, the final result of the
simplification from $f$ is formulated as $(f_{1}<<6)+f_{2}$.

Figure 2:  A bit representation of a heterogeneous-bit-valued contstant and
the variable

##  4 Discussions and Conclusions

We propose a novel MBA simplifier, called “SSLEM,” based on a new concept of
“semi-linear” MBA expression transformation, as well as existing logic based
simplifiers that restrict their target into linear MBA expressions. We assume
that more semi-linear MBA expressions exist in real-world obfuscated malware,
than linear MBA expressions do. The shorter contiguous sequences of the same
bits in a -bit-valued constant complicates the result. For instance, 0101 0101
0101 01 0101 will require 16 separated linear expressions and 15 shifts in the
simplification result, which might not be simpler than the original
expression. However, our method will be effective with the small number of
changes in the bits of constants like 0x01 and 0xFE. By taking the advantages
of program synthesis and logically proved solutions, we envision that SSELM
will show better performance. Since some expressions are likely to be
simplified by program synthesis, while some others are simplified by logic-
based rewriting methods, we expect that SSLEM would be able to handle the
corner cases that the existing solutions could not hanldle.

## References

  * [1] N. Eyrolles, _Obfuscation with Mixed Boolean-Arithmetic Expressions : reconstruction, analysis and simplification tools_ , Cryptography and Security [cs.CR]. Universite Paris-Saclay, 2017.
  * [2] B. Liu et. al, _MBA-Blast: Unveiling and Simplifying Mixed Boolean-Arithmetic Obfuscation_ , Usenix Security 2021
  * [3] Y. Zhou et. al, _Information Hiding in Software with Mixed Boolean-Arithmetic Transforms_ , WISA 2007
  * [4] T. Blazytco et al, _Syntia: Synthesizing the semantics of obfuscated code_ , USENIX Security Security, 2017
  * [5] R. David et., _QSynth-A Program Synthesis based Approach for Binary Code Deobfuscation_ , BAR 2020 Workshop, NDSS 2020
  * [6] R. David et., _Search-Based Local Black-Box Deobfuscation: Understand, Improve and Mitigate_ , ACM CCS 2021
  * [7] _The Tigress C Diversifier/Obfuscator_ https://tigress.wtf/
  * [8] D. Xu, et. al, _Boosting SMT Solver Performance on Mixed-Bitwise-Arithmetic Expressions_ , ACM PLDI 2021
  * [9] D. Xu, et. al, _VMHunt: A Verifiable Approach to Partially-Virtualized Binary Code Simplification_ ACM CCS 2019