# Fully-heavy baryons $QQQ$ in vacuum and hot QCD medium

Jiaxing Zhaoa<EMAIL_ADDRESS>Shuzhe Shib<EMAIL_ADDRESS>aSUBATECH, Université de Nantes, IMT Atlantique, IN2P3/CNRS, 4 rue Alfred
Kastler, 44307 Nantes cedex 3, France
bPhysics Department, Tsinghua University, Beijing 100084, China

###### Abstract

We study the properties of fully-heavy baryons in the vacuum and the hot QCD
medium, which is created in relativistic heavy-ion collisions. Masses and wave
functions of $\Omega_{ccc}$, $\Omega_{ccb}$, $\Omega_{bbc}$, and
$\Omega_{bbb}$ up to the second radial excited states are obtained by solving
the three-body Schrödinger equation with Hyperspherical Harmonics method. With
parameters completely fixed by fitting quarkonium boundstates in vacuum, we
predicted the masses for $1S$, $2S$, and $3S$ states of fully-heavy baryons.
We also computed the temperature dependence of baryon masses and the thermal
widths in a hot QCD medium. These properties are important to precise study of
fully-heavy baryon production in heavy ion collisions.

##  I introduction

Taking into account the fact that charm and bottom quarks are very heavy and
their moving velocity is small, there exists a hierarchy of scales in the
study of heavy quarks: $m\gg mv\gg mv^{2}$ [1, 2]. Integrating out the degrees
of freedom with momenta larger than $m$ and $mv$ successively in the QCD
Lagrangian, one can derive its nonrelativistic versions NRQCD and pNRQCD [2].
Furthermore, if neglecting the color transition between the color-singlet and
color-octet states, the pNRQCD becomes a potential model [2]. In this case,
one can employ the Schrödinger equation to study the properties of hadrons
consisting of only heavy quarks. It turns out that the mass spectra of
quarkonium are explained very well based on the Schrödinger equation with
Cornell potential, see the review paper [3]. The Schrödinger equation has also
been extended to three-body case to predict the masses of fully-heavy baryons
[4, 5, 6, 7, 8, 9].

Searching for the QCD phase transition is one of the physics goals of
relativistic heavy-ion collisions. The phase transition happens around $160$
MeV predicted by the lattice QCD and effective models [10, 11, 12, 13]. The
high temperature phase is the so called quark-gluon plasma (QGP) and has been
confirmed in experiments [14]. This hot QCD medium changes quarkonium
production rate, compared to the vacuum case which happens in $pp$ or
$e^{+}e^{-}$ collisions. Quarkonium suppression was proposed as a smoking gun
of the production of the QGP [15]. The evolution and production of quarkonium
in QGP medium depends closely on their finite-temperature properties, such as
binding energy and width [16, 17, 18, 19, 20, 21]. Analogously, to study the
yield of fully-heavy baryons $QQQ$ in the QGP, their finite-temperature
properties are demanded. Many theoretical studies show the $\Xi_{cc}$ and
$\Omega_{ccc}$ yields per binary nucleon-nucleon collision in heavy-ion
collisions at RHIC and LHC will be largely enhanced in comparison with
nucleon-nucleon collisions in vacuum [22, 23, 24, 25]. Searching for these
fully-heavy baryons in relativistic heavy ion collisions at LHC energy
attracted many attentions and has been listed as one of goal of next-
generation LHC heavy-ion experiment [26]. Therefore, it is important to figure
out their properties in the hot QCD medium, especially with the help of the
newly obtained heavy quark finite-temperature potential [27, 28]. In this
work, we employ the three-body Schrödinger equation to study the properties of
fully-heavy baryon states $\Omega_{ccc}$, $\Omega_{ccb}$, $\Omega_{bbc}$, and
$\Omega_{bbb}$ in both the vacuum and finite temperatures.

The structure of this paper is as follows. In Sec. II we present the framework
of solving the three-body Schrödinger equation. The baryon properties,
including mass and size, in the vacuum and a hot medium, are investigated in
Secs. III and IV, respectively. A summary is given in Sec. V.

##  II Theoretic framework

###  II.1 Coordinate Transformation

For a system of three quarks with the mass $\mu_{i}(i=1,2,3)$, the wave
function $\Psi({\bf r}_{1},{\bf r}_{2},{\bf r}_{3})$ and the energy $E$
satisfy by the Schrödinger equation

$\left(\sum_{i=1}^{3}{\hat{\bf p}^{2}_{i}\over
2\mu_{i}}+\sum_{i<j}V_{ij}(|{\bf r}_{ij}|)\right)\,\Psi=E\,\Psi,$ (1)

under the boundary condition that the wave function vanishs when coordinates
approaches to infinity. We have neglected the direct three-body potentials and
assumed that the interaction potential is the summation of the two-body
interactions. Taking into account one-gluon-exchange interaction, the two-body
potential can be effectively expressed as [29, 30],

$V_{ij}(|{\bf r}_{ij}|)=-{\lambda_{i}^{a}\cdot\lambda_{j}^{a}\over
4}\left(V^{c}_{ij}(|{\bf r}_{ij}|)+V^{ss}_{ij}(|{\bf r}_{ij}|){\bf
s}_{i}\cdot{\bf s}_{j}\right),$ (2)

where $\lambda_{i}^{a}\ (a=1,...,8)$ are the SU(3) Gell-Mann matrices, the
factor $1/4$ is from the normalization, $V^{c}_{ij}$ is the spin independent
interaction, $V^{ss}_{ij}$ is the strength of the spin-spin interaction, and
$|{\bf r}_{ij}|=|{\bf r}_{i}-{\bf r}_{j}|$ is the distance between two quarks
labeled by $i$ and $j$. We employ the Cornell potential to describe the spin-
independent central interaction $V_{ij}^{c}$ between two quarks and the
lattice result [30] for the spin-spin coupling,

$\displaystyle V^{c}_{ij}(|{\bf r}_{ij}|)$ $\displaystyle=$
$\displaystyle-{\alpha\over|{\bf r}_{ij}|}+\sigma|{\bf r}_{ij}|,$
$\displaystyle V^{ss}_{ij}(|{\bf r}_{ij}|)$ $\displaystyle=$
$\displaystyle\beta e^{-\gamma|{\bf r}_{ij}|}.$ (3)

The parameters in the potential like $\alpha$, $\sigma$, $\beta$, and $\gamma$
are given by lattice QCD with some uncertainties, which can be fixed by
fitting the experimental data of charmonium and bottomonium masses. We will
show this later.

In order to solve the three-body Schrödinger equation, we first introduce the
Jacobi coordinates,

$\displaystyle{\bf R}$ $\displaystyle=$ $\displaystyle{1\over\mu}(\mu_{1}{\bf
r}_{1}+\mu_{2}{\bf r}_{2}+\mu_{3}{\bf r}_{3}),$ $\displaystyle{\bf x}_{1}$
$\displaystyle=$
$\displaystyle{\sqrt{(\mu_{1}+\mu_{2})\mu_{3}}\over\mu}\left({\bf
r}_{3}-{\mu_{1}{\bf r}_{1}+\mu_{2}{\bf r_{2}}\over\mu_{1}+\mu_{2}}\right),$
$\displaystyle{\bf x}_{2}$ $\displaystyle=$
$\displaystyle\sqrt{{\mu_{1}\mu_{2}\over(\mu_{1}+\mu_{2})\mu}}\left({\bf
r}_{2}-{\bf r}_{1}\right),$ (4)

where $\mu\equiv\sum_{i=1}^{3}\mu_{i}$ is the total mass. With such
coordinates, the kinetic energy becomes

$\sum_{i=1}^{3}{\hat{\bf p}^{2}_{i}\over 2\mu_{i}}=\frac{\hat{\bf
P}^{2}+\hat{\bf q}_{1}^{2}+\hat{\bf q}_{2}^{2}}{2\mu},$ (5)

where $\hat{\bf q}_{1}$ and $\hat{\bf q}_{2}$ are the relative momenta that
conjugate to ${\bf x}_{1}$ and ${\bf x}_{2}$, respectively.

Since the potential depends only on the relative coordinates ${\bf x}_{i}$ and
the total momentum is conserved, one can factorize the three-body motion into
a center-of-mass motion and a relative motion, $\Psi({\bf r}_{1},{\bf
r}_{2},{\bf r}_{3})=e^{i{\bf P\cdot R}}\Phi({\bf x}_{1},{\bf x}_{2})$. The
bound state properties only relate to the relative motion of the system, and
we just need to deal with the six-dimensional wave equation. We then express
the relative coordinates ${\bf x}_{1}$ and ${\bf x}_{2}$ in the hyperspherical
frame [31]: hyperradius $\rho=\sqrt{x_{1}^{2}+x_{2}^{2}}$ and hyperangles
$\Omega=\\{\alpha_{2},\theta_{1},\phi_{1},\theta_{2},\phi_{2}\\}$, where the
angle $\alpha_{2}\equiv\arcsin(x_{2}/\rho)$ is defined within the range
$[0,\pi/2]$, and $\\{x_{i},\theta_{i},\phi_{i}\\}$ are the spherical
coordinates corresponding to ${\bf x}_{i}$. With the hyperspherical
coordinates, the Schrödinger equation governing the relative wave function
$\Phi(\rho,\Omega)$ can be written as

$\displaystyle\left[{1\over 2\mu}\left(-{d^{2}\over
d\rho^{2}}-{5\over\rho}{d\over d\rho}+{\widehat{\bf
K}^{2}\over\rho^{2}}\right)+V(\rho,\Omega)\right]\Phi=E_{r}\Phi,$ (6)

where the corresponding energy eigenvalue $E_{r}=E-\frac{{\bf P}^{2}}{2M}$ and
$\widehat{\bf K}$ is the hyperangular momentum operator. The hyperspherical
harmonic(HH) functions ${\mathcal{Y}}_{\kappa}(\Omega)$ are the eigenstates of
$\widehat{\bf K}^{2}$,

$\widehat{\bf
K}^{2}{\mathcal{Y}}_{\kappa}(\Omega)=K(K+4){\mathcal{Y}}_{\kappa}(\Omega),$
(7)

where $K$, referred to as the grand-orbital momentum, is the quantum number
describing the magnitude of the hyper-angular momentum. There are eight
operators that commute with the kinetic energy term in the Hamiltonian and
with each other. In addition to $K$, the other conserved quantum numbers are
total angular momentum($L$), total magnetic quantum number($M$), angular
momentum of corresponding to each Jacobi coordinate($l_{1}$ and $l_{2}$), and
$n$. Here $l_{1}$($l_{2}$) is the orbital angular momentum quantum number of
subsystem 1 and 2 (12 diquark and 3). $n$ is the relative radial quantum
number. These quantum numbers satisfy,

$\displaystyle K$ $\displaystyle=$ $\displaystyle 2n+L,$ $\displaystyle L$
$\displaystyle=$ $\displaystyle l_{1}+l_{2},$ $\displaystyle M$
$\displaystyle=$ $\displaystyle m_{1}+m_{2}.$ (8)

We have introduced the shorthand that $\kappa\equiv\\{K,L,M,n,l_{1},l_{2}\\}$.
In a three-body system, the HH function can be expressed as,

$\displaystyle\begin{split}{\mathcal{Y}}_{\kappa}(\Omega)=&\sum_{m_{1},m_{2}}\langle
l_{1}m_{1}l_{2}m_{2}|LM\rangle
Y_{l_{1}}^{m_{1}}(\theta_{1},\phi_{1})Y_{l_{2}}^{m_{2}}(\theta_{2},\phi_{2})\\\
&\times\left[{\mathcal{N}}\,\sin^{l_{2}}\alpha_{2}\cos^{l_{1}}\alpha_{2}\,P_{n}^{l_{2}+\frac{1}{2},l_{1}+\frac{1}{2}}(\cos
2\alpha_{2})\right],\end{split}$ (9)

with $Y_{l}^{m}$ being the spherical harmonics, $P_{n}^{l}$ the associated
Legendre polynomials, and

$\displaystyle{\mathcal{N}}=\sqrt{(2K+4)n!\Gamma(n+l_{1}+l_{2}+2)\over\Gamma(n+l_{2}+{3\over
2})\Gamma(n+l_{1}+{3\over 2})}.$ (10)

See e.g., [32, 33, 34] for properties of the HH functions.

###  II.2 Spatial, Color, and Spin Wavefunctions

As shown in (2), the potential $V(\rho,\Omega)$ depends on the color and spin
degrees of freedom. We start to construct the color and spin wave-function,
based on the symmetry properties. Identical fermions should fulfill the Pauli
principle, the total wave-function should be antisymmetric under exchange. For
heavy quarks, the flavor wave function is trivial. So, the total wave function
can be expressed as,

$\displaystyle\Psi=\psi_{\mathrm{space}}\,\phi_{\mathrm{color}}\,\chi_{\mathrm{spin}},$
(11)

where the color wave-function
$\phi_{\mathrm{color}}=(QQ)_{\bar{3}_{c}}Q_{3_{c}}$ is antisymmetry. For the
spin space, one has,

$\displaystyle 2\otimes 2\otimes 2=(3\oplus 1)\otimes 2=4\oplus 2\oplus 2.$
(12)

As a result, the othonormal basis is listed as follows,

$\displaystyle\chi_{+\frac{3}{2}}^{s}$ $\displaystyle=$
$\displaystyle\left|\frac{3}{2},+\frac{3}{2}\right>_{\mathrm{S}}=\left|\uparrow\uparrow\uparrow\right>$
$\displaystyle\chi_{+\frac{1}{2}}^{s}$ $\displaystyle=$
$\displaystyle\left|\frac{3}{2},+\frac{1}{2}\right>_{\mathrm{S}}=3^{-1/2}\Big{[}\left|\uparrow\uparrow\downarrow\right>+\left|\uparrow\downarrow\uparrow\right>+\left|\downarrow\uparrow\uparrow\right>\Big{]}$
$\displaystyle\chi_{-\frac{1}{2}}^{s}$ $\displaystyle=$
$\displaystyle\left|\frac{3}{2},-\frac{1}{2}\right>_{\mathrm{S}}=3^{-1/2}\Big{[}\left|\downarrow\downarrow\uparrow\right>+\left|\downarrow\uparrow\downarrow\right>+\left|\uparrow\downarrow\downarrow\right>\Big{]}$
$\displaystyle\chi_{-\frac{3}{2}}^{s}$ $\displaystyle=$
$\displaystyle\left|\frac{3}{2},-\frac{3}{2}\right>_{\mathrm{S}}=\left|\downarrow\downarrow\downarrow\right>$
$\displaystyle\chi_{+\frac{1}{2}}^{ms}$ $\displaystyle=$
$\displaystyle\left|\frac{1}{2},+\frac{1}{2}\right>_{\mathrm{MS}}=6^{-1/2}\Big{[}\left|\uparrow\downarrow\uparrow\right>+\left|\downarrow\uparrow\uparrow\right>-2\left|\uparrow\uparrow\downarrow\right>\Big{]}$
$\displaystyle\chi_{-\frac{1}{2}}^{ms}$ $\displaystyle=$
$\displaystyle\left|\frac{1}{2},-\frac{1}{2}\right>_{\mathrm{MS}}=6^{-1/2}\Big{[}2\left|\downarrow\downarrow\uparrow\right>-\left|\downarrow\uparrow\downarrow\right>-\left|\uparrow\downarrow\downarrow\right>\Big{]}$
$\displaystyle\chi_{+\frac{1}{2}}^{ma}$ $\displaystyle=$
$\displaystyle\left|\frac{1}{2},+\frac{1}{2}\right>_{\mathrm{MA}}=2^{-1/2}\Big{[}\left|\uparrow\downarrow\uparrow\right>-\left|\downarrow\uparrow\uparrow\right>\Big{]}$
$\displaystyle\chi_{-\frac{1}{2}}^{ma}$ $\displaystyle=$
$\displaystyle\left|\frac{1}{2},-\frac{1}{2}\right>_{\mathrm{MA}}=2^{-1/2}\Big{[}\left|\uparrow\downarrow\downarrow\right>-\left|\downarrow\uparrow\downarrow\right>\Big{]}$
(13)

Here the $ms$ or $ma$ represent exchange antisymmetry or symmetry between 1
and 2 particle. It is easy to get,

$\displaystyle{\bf s}_{i}\cdot{\bf s}_{j}|\chi^{s}\rangle$ $\displaystyle=$
$\displaystyle{1\over 4}|\chi^{s}\rangle,$ $\displaystyle{\bf s}_{1}\cdot{\bf
s}_{2}|\chi^{ms}\rangle$ $\displaystyle=$ $\displaystyle{1\over
4}|\chi^{ms}\rangle,$ $\displaystyle{\bf s}_{1}\cdot{\bf
s}_{3}|\chi^{ms}\rangle$ $\displaystyle=$ $\displaystyle-{1\over
2}|\chi^{ms}\rangle-{\sqrt{3}\over 4}|\chi^{ma}\rangle,$ $\displaystyle{\bf
s}_{2}\cdot{\bf s}_{3}|\chi^{ms}\rangle$ $\displaystyle=$
$\displaystyle-{1\over 2}|\chi^{ms}\rangle+{\sqrt{3}\over
4}|\chi^{ma}\rangle,$ $\displaystyle{\bf s}_{1}\cdot{\bf
s}_{2}|\chi^{ma}\rangle$ $\displaystyle=$ $\displaystyle-{3\over
4}|\chi^{ma}\rangle,$ $\displaystyle{\bf s}_{1}\cdot{\bf
s}_{3}|\chi^{ma}\rangle$ $\displaystyle=$ $\displaystyle-{\sqrt{3}\over
4}|\chi^{ms}\rangle,$ $\displaystyle{\bf s}_{2}\cdot{\bf
s}_{3}|\chi^{ma}\rangle$ $\displaystyle=$ $\displaystyle{\sqrt{3}\over
4}|\chi^{ms}\rangle,$ (14)

where the subscript has been omitted.

Noting that the ground state spacial wavefunction is symmetric, we focus on
$S$-wave states (with $L=0$) in this work. For $\Omega_{ccc}$ and
$\Omega_{bbb}$, the spin wavefunction must be symmetric. Therefore, only the
$|\phi_{c}\chi^{s}\rangle$ state with $J^{P}={3\over 2}^{+}$ will be
considered. The potential elements in color-spin space can be expressed as,

$\displaystyle\begin{split}&\langle\phi_{c}\chi^{s}|\sum_{i<j}V_{ij}|\phi_{c}\chi^{s}\rangle\\\
=\;&\frac{V^{c}_{12}+V^{c}_{13}+V^{c}_{23}}{2}+\frac{V^{ss}_{12}+V^{ss}_{13}+V^{ss}_{23}}{8}.\end{split}$
(15)

For the $\Omega_{ccb}$ and $\Omega_{bbc}$, the spin wavefunction should be
symmetric when exchanging the first two particles. The ground states can be
$|\phi_{c}\chi^{s}\rangle$ with $J^{P}={3\over 2}^{+}$or
$|\phi_{c}\chi^{ms}\rangle$ with $J^{P}={1\over 2}^{+}$. The potential
elements of the former are given in Eq. (15), whereas those of the latter can
be written as,

$\displaystyle\begin{split}&\langle\phi_{c}\chi^{ms}|\sum_{i<j}V_{ij}|\phi_{c}\chi^{ms}\rangle\\\
=\;&\frac{V^{c}_{12}+V^{c}_{13}+V^{c}_{23}}{2}+\frac{V^{ss}_{12}}{8}-\frac{V^{ss}_{13}+V^{ss}_{23}}{4}.\end{split}$
(16)

###  II.3 Hyperspherical Harmonics Expansion

The potential $V(\rho,\Omega)$ depends not only on the hyperradius but also on
the eight hyperangles. The Schrödinger equation (6) cannot be further
factorized into a radial part and an angular part. Instead, one expands the
wave function in terms of the hyper-spherical harmonic(HH) functions
${\mathcal{Y}}_{\kappa}(\Omega)$. The total relative wave function can be
expanded as

$\displaystyle\Phi(\rho,\Omega)=\sum_{\kappa}R_{\kappa}(\rho){\mathcal{Y}}_{\kappa}(\Omega)|\phi_{c}\chi^{s}\rangle,$
(17)

for the ${3\over 2}^{+}$ states, and

$\displaystyle\Phi(\rho,\Omega)=\sum_{\kappa}R_{\kappa}(\rho){\mathcal{Y}}_{\kappa}(\Omega)|\phi_{c}\chi^{ms}\rangle\,,$
(18)

for the ${1\over 2}^{+}$ states. Taking such an expansion and employing the
reduced radial wavefunction as
$u_{\kappa}(\rho)\equiv\rho^{5/2}R_{\kappa}(\rho)$, the multi-variable
Schröedinger equation becomes a set of coupled differential equations of
single variable,

$\displaystyle\left[{1\over 2\mu}{d^{2}\over d\rho^{2}}-{4K(K+4)+15\over
8\mu\,\rho^{2}}+E_{r}\right]u_{\kappa}=\sum_{\kappa^{\prime}}V_{\kappa\kappa^{\prime}}u_{\kappa^{\prime}},$
(19)

where $V_{\kappa\kappa^{\prime}}$ is the potential matrix,

$\displaystyle\begin{split}V_{\kappa\kappa^{\prime}}=\;&\int{\mathcal{Y}}_{\kappa}^{*}(\Omega)V(\rho,\Omega){\mathcal{Y}}_{\kappa^{\prime}}(\Omega)d\Omega\\\
=\;&\sum_{i<j}\int
V_{ij}(|\mathbf{r}_{ij}|){\mathcal{Y}}_{\kappa}^{*}(\Omega){\mathcal{Y}}_{\kappa^{\prime}}(\Omega)\mathrm{d}\Omega\,,\end{split}$
(20)

with the volume element being

$d\Omega=\cos^{2}\alpha_{2}\sin^{2}\alpha_{2}\sin\theta_{1}\sin\theta_{2}d\alpha_{2}d\theta_{1}d\phi_{1}d\theta_{2}d\phi_{2}.$
(21)

###  II.4 Computation of the Potential Matrix

Computing the potential matrix is nontrivial. In the most general form, (20)
is a five-dimensional integral, which is computationally expensive. However,
taking the assumption that the total interaction potential is the summation of
two-body interaction $V_{ij}(|{\bf r}_{ij}|)$, one can reduce (20) into a one-
dimensional integral by performing particle permutation. Let us first focus on
$V_{12}$, which depends only on ${\bf x}_{2}$, see Eq. (4). The integral over
$\alpha_{2}$ is the only non-trivial one out of the five hyperangles, whereas
the remainders can be computed analytically using the orthogonal relation of
the HH functions. To compute other sectors of the potential matrix, the trick
is to consider another form of Jacobi coordinate, which differs from the
previous one by a particle permutation and puts ${\bf r}_{2}-{\bf r}_{3}$ or
${\bf r}_{3}-{\bf r}_{1}$ in the newly defined ${\bf x}_{2}$. We label the
covention in Eq. (4) as ${\bf x}^{(12)}_{i}$, and the other two conventions as
${\bf x}^{(23)}_{i}$ and ${\bf x}^{(31)}_{i}$. Explicitly, they are defined as

$\displaystyle{\bf x}^{(23)}_{1}$ $\displaystyle=$
$\displaystyle{\sqrt{(\mu_{2}+\mu_{3})\mu_{1}}\over\mu}\left({\bf
r}_{1}-{\mu_{2}{\bf r}_{2}+\mu_{3}{\bf r_{3}}\over\mu_{2}+\mu_{3}}\right),$
$\displaystyle{\bf x}^{(23)}_{2}$ $\displaystyle=$
$\displaystyle\sqrt{{\mu_{2}\mu_{3}\over(\mu_{2}+\mu_{3})\mu}}\left({\bf
r}_{3}-{\bf r}_{2}\right),$ (22)

and

$\displaystyle{\bf x}^{(31)}_{1}$ $\displaystyle=$
$\displaystyle{\sqrt{(\mu_{3}+\mu_{1})\mu_{2}}\over\mu}\left({\bf
r}_{2}-{\mu_{3}{\bf r}_{3}+\mu_{1}{\bf r_{1}}\over\mu_{3}+\mu_{1}}\right),$
$\displaystyle{\bf x}^{(31)}_{2}$ $\displaystyle=$
$\displaystyle\sqrt{{\mu_{3}\mu_{1}\over(\mu_{3}+\mu_{1})\mu}}\left({\bf
r}_{1}-{\bf r}_{3}\right),$ (23)

which can be related to ${\bf x}^{(12)}$ by

$\displaystyle\left(\begin{array}[]{c}{\bf x}_{1}^{(12)}\\\ {\bf
x}_{2}^{(12)}\\\ \end{array}\right)=\left(\begin{array}[]{cc}a_{11}&a_{12}\\\
a_{21}&a_{22}\\\ \end{array}\right)\cdot\left(\begin{array}[]{c}{\bf
x}_{1}^{(ij)}\\\ {\bf x}_{2}^{(ij)}\\\ \end{array}\right),$ (30)

where $a_{ij}$ is element of the rotation matrix.

It has been shown that under the particle permutation, different HH’s are
related by an unitary transformation,

$\displaystyle{\mathcal{Y}}_{\kappa}(\Omega)=\sum_{\kappa^{\prime}}R^{(ij)}_{\kappa\kappa^{\prime}}{\mathcal{Y}}_{\kappa^{\prime}}(\Omega^{(ij)}),$
(31)

and the coefficients $R^{(ij)}_{\kappa\kappa^{\prime}}$ are refereed to as the
Raynal–Revai coefficients [35]. The Raynal–Revai coefficient is non vanishing
only for the HH’s with the same grand-orbital momentum($K$), total angular
momentum($L$), and magnetic quantum number($M$). Details of the Raynal–Revai
coefficients, especially its explicit form for $K=0,2,4,6$ are shown in
Appendix A.

With these preparation, we are now ready to compute the potential matrix for
an arbitrary $(i,j)$ pair:

$\displaystyle V_{\kappa\kappa^{\prime}}^{(ij)}(\rho)$ $\displaystyle\equiv$
$\displaystyle\int
V^{(ij)}(|\mathbf{r}_{ij}|){\mathcal{Y}}_{\kappa}^{*}(\Omega){\mathcal{Y}}_{\kappa^{\prime}}(\Omega)\mathrm{d}\Omega$
(32) $\displaystyle=$ $\displaystyle\sum_{kk^{\prime}}(R^{(ij)}_{\kappa
k})^{*}R^{(ij)}_{\kappa^{\prime}k^{\prime}}$ $\displaystyle\times$
$\displaystyle\int
V^{(ij)}(|\mathbf{r}_{ij}|){\mathcal{Y}}_{k}^{*}(\Omega^{(ij)}){\mathcal{Y}}_{k^{\prime}}(\Omega^{(ij)})\mathrm{d}\Omega^{(ij)}$
$\displaystyle=$ $\displaystyle\sum_{kk^{\prime}}(R^{(ij)}_{\kappa
k})^{*}R^{(ij)}_{\kappa^{\prime}k^{\prime}}V_{kk^{\prime}}^{(12)}(\rho).$

##  III Baryons in vacuum

Table 1: Potential model parameters

$m_{b}$ | $m_{c}$ | $\alpha$ | $\sigma$
---|---|---|---
4.7 GeV | 1.29 GeV | 0.308 | 0.15 GeV2
$\gamma$ | $\beta_{b\bar{b}}$ | $\beta_{c\bar{c}}$ |  $\beta_{b\bar{c}}$=$\beta_{c\bar{b}}$
1.982 GeV | 0.239 GeV | 1.545 GeV | 0.525 GeV

Table 2: The experimental [36] and calculated quarkonium masses

State |  $M_{\rm exp}$(GeV) |  $M_{\rm theo}$(GeV)
---|---|---
$\eta_{c}$ | 2.981 | 2.968
$J/\psi$ | 3.097 | 3.102
$h_{c}(1P)$ | 3.525 | 3.480
$\chi_{c}(1P)$ | 3.556 | 3.500
$\eta_{c}(2S)$ | 3.639 | 3.654
$\psi(2S)$ | 3.696 | 3.720
$\eta_{b}$ | 9.398 | 9.397
$\Upsilon(1S)$ | 9.460 | 9.459
$h_{b}(1P)$ | 9.898 | 9.845
$\chi_{b}(1P)$ | 9.912 | 9.860
$\eta_{b}(2S)$ | 9.999 | 9.957
$\Upsilon(2S)$ | 10.023 | 9.977
$B_{c}(^{1}S_{0})$ | 6.275 | 6.282
$B_{c}(^{3}S_{1})$ | - | 6.347
$B_{c}(^{1}P_{1})$ | - | 6.726
$B_{c}(^{3}P_{0,1,2})$ | - | 6.738
$B_{c}(2^{1}S_{0})$ | 6.871 | 6.886
$B_{c}(2^{3}S_{1})$ | - | 6.915

Table 3: The calculated fully heavy flavor baryon mass $M_{\rm theo}$ and the root-mean-squared radius $r_{\text{rms}}$ for the ground $1S$ and radial-excited states, $2S$ and $3S$. | $J^{P}={3\over 2}^{+}$ | $J^{P}={1\over 2}^{+}$
---|---|---
| $\Omega_{ccc}$ | $\Omega_{ccb}$ | $\Omega_{bbc}$ | $\Omega_{bbb}$ | $\Omega_{ccb}$ | $\Omega_{bbc}$
State | $1S$ | $2S$ | $3S$ | $1S$ | $2S$ | $3S$ | $1S$ | $2S$ | $3S$ | $1S$ | $2S$ | $3S$ | $1S$ | $2S$ | $3S$ | $1S$ | $2S$ | $3S$
$M_{\rm theo}$(GeV) | 4.80 | 5.31 | 5.74 | 8.19 | 8.68 | 9.08 | 10.96 | 11.33 | 11.63 | 14.36 | 14.77 | 15.09 | 8.17 | 8.66 | 9.07 | 10.87 | 11.28 | 11.59
$r_{\text{rms}}$(fm) | 0.30 | 0.45 | 0.57 | 0.22 | 0.33 | 0.43 | 0.19 | 0.31 | 0.41 | 0.16 | 0.26 | 0.34 | 0.22 | 0.33 | 0.42 | 0.18 | 0.30 | 0.40
Figure 1: (Left) The vacuum radial probabilities for $1S$(top), $2S$(mid), and
$3S$(bottom) states of $\Omega_{ccc}$(solid) and $\Omega_{bbb}$(dotted) with
quantum number $J^{P}=(3/2)^{+}$. The same for $\Omega_{ccb}$(middle) and
$\Omega_{bbc}$(right) particles with $J^{P}=(3/2)^{+}$(solid) and
$J^{P}=(1/2)^{+}$(dotted). Lines with different colors correspond to
$K=0$(red), $2$(gold), $4$(blue), and $6$(purple). See texts for detailed
explanation of the states.

In this section, we start with computing the heavy flavor baryons $QQQ$ bound
states in a vacuum. As mentioned before, we first fix the parameters in
potential and the quark masses ($m_{c}$ and $m_{b}$) with the knew charmonium,
bottomonium, and also $B_{c}$ mesons data. Their masses can be calculated via
the two-body Schrödinger equation with the potential

$V_{Q\bar{Q}}={4\over 3}\left(V^{c}_{ij}(r)+V^{ss}_{ij}(r){\bf s}_{i}\cdot{\bf
s}_{j}\right),$ (33)

where $4/3$ is the color factor for color-singlet states $Q\bar{Q}$. With the
model parameters presented in Table 2, we obtain the quarkonium masses and
$B_{c}$ mesons shown in Table 2. Here the parameters and quarkonium masses
have been used and calculated in previous work [37]. One can see that, the
quarkonia mass can be well described.

With the known parameters, we then solve the three-body Schrödinger equations
(19) for fully-heavy baryons $QQQ$. The baryon mass comes from the summation
of the constituent masses $\mu=\sum_{i=1}^{3}\mu_{i}$ and the binding energy
$E_{r}$ which is determined by the radial equations,

$M_{B}=\mu+E_{r}\,,$ (34)

and the root-mean-squared radius is defined as

$\displaystyle
r_{\text{rms}}^{2}=\int\sum_{\kappa}|R_{\kappa}(\rho)|^{2}\rho^{7}d\rho.$ (35)

One may find $r_{\text{rms}}^{2}=\big{\langle}{1\over 3}\sum_{i=1}^{3}({\bf
r}_{i}-{\bf X})^{2}\big{\rangle}$ when all the three quarks have the same
mass. The normalization
$\int\sum_{\kappa}|R_{\kappa}(\rho)|^{2}\rho^{5}d\rho=1$ for the radial
functions $R_{\kappa}(\rho)$.

One can only include a finite number of hyperspherical harmonics in a
numerical calculation, and our truncation is made according to the symmetry
properties of the system. Since we focus in this work on the $S$-wave baryon
states, the relevant hyperspherical harmonics are those corresponding to
vanishing total orbital angular momentum $L$ and magnetic quantum number $M$,
i.e. $L=M=0$. We choose all such hyperspherical harmonic functions with
hyperangular quantum number $K\leq 6$. This leads to coupled differential
equations which are numerically solved by using the inverse power method [38].
The main advantage of taking the inverse power method is its high precision
for both ground and excited states.

We show the baryon wavefunctions for the ground state $1S$, first radial
excited state $2S$, and the second radial excited state $3S$ states of
$\Omega_{ccc}$, $\Omega_{ccb}$, $\Omega_{bbc}$, and $\Omega_{bbb}$ in Fig. 1
and present the mass and root-mean-squared radius in Table 3. There are
respectively one, two, three, and four states with quantum number $K=0$, $2$,
$4$, and $6$. They are respectively represented by red, gold, blue, and purple
curves in Fig. 1. The clear hierarchy of their magnitudes shows the
convergence of the hyperspherical harmonic expansion. For the equal-mass
system, i.e., $\Omega_{ccc}$ and $\Omega_{bbb}$, we notice the vanishment of
the $K=2$ states, which always correspond to $n=1$. When interchanging the
coordinates of two quarks, $S$ states of the equal-mass system is always
symmetric whereas $n=1$ states are anti-symmetric. Thus, $\Omega_{ccc}$ and
$\Omega_{bbb}$ do not contain $K=2$ component. For the unequal-mass baryons,
i.e., $\Omega_{ccb}$ and $\Omega_{bbc}$, the $K=2$ states become non-
vanishing, and the contributions of other higher order hyperspherical harmonic
functions are more sizable. Comparing the spin half states of $\Omega_{ccb}$
and $\Omega_{bbc}$ are respectively lighter than the corresponding three-
halves states, we find the spin half states are, respectively, lighter than
the corresponding three-halves state. The modest differences in the mass
difference and the wave functions indicate the spin-spin interaction is a
higher-order effect compared to the central interaction. We also check the
convergence of hyperspherical harmonics expansion by observing same results,
up to desired accuracy, when keeping states with $K\leq 4$ only. The ground
state masses of these fully heavy baryons are consistent with previous studies
based on lattice QCD [39, 40, 41, 42], the potential model [4, 5, 6, 7, 8, 9],
and other effective theories [43, 44, 45, 46].

##  IV Baryons in hot medium

Figure 2: Distance dependence of real(upper) and temperature-scaled
imaginary(lower) potentials taking potential I [47, 27](colored dashed) and II
[28](black solid). Results for Potential II are temperature independent.

As mentioned before, a hot QCD medium—quark-gluon plasma(QGP) is created in
relativistic heavy ion collisions. The typical temperature of the QGP is
hundreds of MeV, estimated by the spectrum of the direct photon [48], which is
much larger than the binding energy of most hadrons. Only the tightly bound
states of heavy quarks, such as $J/\psi$ and $\Upsilon$, can survive in the
QGP but with a large thermal width, as shown in lattice results [21, 49, 50].
In this section, we come to study the finite-temperature properties of fully
heavy baryon states.

Figure 3: (Left to right) Masses(upper) and thermal widths(lower) for
$1S$(red), $2S$(green), and $3S$(blue) states of $\Omega_{ccc}$,
$\Omega_{ccb}$, $\Omega_{bbc}$, and $\Omega_{bbb}$ as functions of temperature
$T$. The solid lines are with potential II [28], while the dashed lines are
with potential I [47, 27].

The finite temperature properties of quarkonium states are encoded in the
finite temperature potential between heavy quark $Q$ and antiquark $\bar{Q}$.
For the baryons, there are no study on their finite-temperature potential in
either weak or strong coupling regions. In the weak coupling limit, the HTL
study shows the heavy quark potential is proportional to the color factor
[51]. So, for the color wave-function $(QQ)_{\bar{3}_{c}}Q_{3_{c}}$ of the
baryons, both the real and imaginary potential still satisfy the same relation
as shown in Eq. (2). Besides, we neglect the influence of hot QCD medium to
the heavy quark spin-spin interaction.

The heavy-quark interaction potential is screened by other color objects in
the QGP, and consequently the long-range interaction (3) is strongly
suppressed when the temperature is high enough, as shown in the hard-thermal-
loop (HTL) perturbative calculation [51]. Besides, the potential develops an
imaginary part which originates from the Landau-damping [51]. While in the
strong coupling region, the heavy quark potential can be extracted from the
Wilson loop in lattice QCD [47, 27, 28]. The result also shows the heavy quark
potential is complex-valued. The real part potential shows the screening
effect, while the imaginary part reflect the decay of quarkonium under
scattering with thermal partons. However, the value of the potential depends
on the extraction strategies, the detail as shown in [28]. To account for the
theoretical uncertainty, we choose two qualitatively different schemes of
finite-temperature potentials, both of which are given by recent lattice QCD
calculation, in this study.

The first potential scheme is the from Ref. [27], which has a obvious color-
screened real potential, named potential (I). Both the real and imaginary
potentials can be fitted by a functional form based on the Gauss’ law approach
[21],

$\displaystyle{\rm Re}V(T,r)$ $\displaystyle=$
$\displaystyle-\alpha\left(m_{D}+{e^{-m_{D}r}\over r}\right)$ $\displaystyle+$
$\displaystyle{\sigma\over m_{D}}\left(2-(2+m_{D}r)e^{-m_{D}r}\right),$
$\displaystyle{\rm Im}V(T,r)$ $\displaystyle=$ $\displaystyle-\alpha
T\phi(m_{D}r)$ $\displaystyle-$ $\displaystyle{\sqrt{\pi}\over 4}m_{D}T\sigma
r^{3}G_{2,4}^{2,2}\left({}_{{1\over 2},{1\over 2},-{3\over 2},-1}^{-{1\over
2},-{1\over 2}}\Big{|}{1\over 4}m_{D}^{2}r^{2}\right)$

with

$\displaystyle\phi(x)=2\int_{0}^{\infty}dz{z\over(z^{2}+1)^{2}}\left(1-{\sin(xz)\over
xz}\right),$ (37)

where $G$ is the Meijer-G function. When going into a high-temperature region,
${\rm Im}V(r)$ is consistent with the result from the pure HTL theory [51].
Both ${\rm Re}V$ and ${\rm Im}V$ only depend on a single temperature dependent
parameter, the Debye mass $m_{D}(T)$, which can be obtained by fitting the
lattice data [47, 27]. $\alpha$ and $\beta$ in the potential is the same as
the vacuum case.

The second potential scheme is from Ref. [28], which is the result of Lattice
QCD with dynamical fermions and extracted by assuming Gaussian spectral
function. We call it potential (II). In this case, the real part potential
exhibits negligible screening effect even up to $\sim 700$ MeV, and it can be
parameterized via Eq. (3). Meanwhile, the imaginary part is found to follow a
simple form that $V_{I}/T$ is a single-variable function of $r\,T$. For the
relevant region, it value is greater that the potential (I). Potential (II) is
quantitatively consistent with the deep-learning extraction of heavy quark
potential [52] from lattice QCD results of masses and widths [49, 50]. The
temperature-dependent real and imaginary parts of the potential (I) and (II)
are shown by colored and black curves, respetively, in Fig. 2.

With the complex potential, we solve the coupled radial equations (19) and
both the energy eigenvalues and wavefunctions are also complex-valued. The
real part of the energy eigenvalue gives the baryon mass: $M_{B}(T)=\mu+{\rm
Re}[E_{r}(T)]$, while the imaginary part gives the thermal width of the
baryons $\Gamma(T)=-{\rm Im}[E_{r}(T)]$. The results for fully heavy baryons,
$\Omega_{ccc}$, $\Omega_{ccb}$, $\Omega_{bbc}$, and $\Omega_{bbb}$, are shown
in Fig. 3. First, we can see the large difference in masses and widths between
two potential schemes. With potential (I), masses decrease with the
temperature while the thermal widths are generated with temperature increases.
With potential (II), we observe weaker temperature-dependence of masses,
owning to fact that the mass is mostly controlled by the real part potential.
The thermal widths are increasing with temperature and are quantitatively
larger than the those using potential (I). For the $\Omega_{ccc}$ and
$\Omega_{ccb}$ the second radial excited states $3S$ disappear when the
temperature high than $0.15$ GeV as presented when using the potential (I).
Second, we can see that replace one charm to the bottom quark, the thermal
decay width of $\Omega_{ccb}$ is very close to $\Omega_{ccc}$ (also for
$\Omega_{bbc}$ and $\Omega_{bbb}$). Comparing to the widths of quarkonium with
potential (I) computed in Ref. [21], we find the thermal width of
$\Omega_{ccc}$ is obviously larger than $J/\psi$, while $\Omega_{bbb}$ is
almost same as the $\Upsilon$. The results indicate most of the fully-heavy
baryons can also survive in a hot QCD medium for a long time in both cases.
Meanwhile, with potential (I), there exists an obvious dissociation
temperature $T_{d}$. When the medium temperature is higher than $T_{d}$, the
fully-heavy baryon will disappear immediately. With the potential (II), the
fully-heavy baryon can still survive at high temperatures as long as the
temperature of the medium drops down fast. These two different behaviors will
be reflected in the transport and final production of fully-heavy baryons in
relativistic heavy ion collisions. Precisely studying the yield of these
baryons in the QGP should go beyond the hadronization at phase transition
hypersurface as in previous studies [22, 23, 24, 25].

##  V Summary

In this work, we study the properties of fully-heavy baryons in the vacuum and
the finite temperature, which is created in high-energy nuclear collisions. We
employ the Hyperspherical Harmonics expansion and solve the three-body
Schrödinger equation. We obtain the masses and wave functions of
$\Omega_{ccc}$, $\Omega_{ccb}$, $\Omega_{bbc}$, and $\Omega_{bbb}$, up to
second radial excited states, for both zero and finite temperatures. In the
vacuum, the ground state masses of predicted are consistent with other
studies. In a hot medium, the temperature dependence of the baryon masses and
the thermal decay widths are calculated. These properties are important to
their productions in heavy ion collisions, as they are expected to affect the
yield and momentum distributions of heavy-flavor hadrons in high-energy
nuclear collisions.

Acknowledgement: The work is supported by the European Union’s Horizon 2020
research and the innovation program under grant agreement No. 824093
(STRONG-2020) (J.Z.) and Tsinghua University under grant No. 53330500923
(S.S.).

##  Appendix A Raynal–Revai coefficient

The Raynal–Revai coefficient is non vanishing only for the HH’s with the same
grand-orbital momentum($K$), total angular momentum($L$), and magnetic quantum
number($M$). A general expression for the Raynal–Revai coefficients is
obtained in the form of a multiple sums over the powers of the coefficients
$a_{ij}$,

$\displaystyle R^{(ij)}_{\kappa\kappa^{\prime}}$ $\displaystyle\equiv$
$\displaystyle\langle l_{i}^{\prime},l_{j}^{\prime}|l_{i},l_{j}\rangle_{K,L}$
$\displaystyle=$
$\displaystyle(-1)^{n+n^{\prime}}(C_{l_{i},l_{j}}^{n}C_{l_{i}^{\prime},l_{j}^{\prime}}^{n^{\prime}})^{-1/2}\sum_{l_{1},l_{2},l_{3},l_{4}}i^{l_{2}-l_{1}+l_{j}-l_{j}^{\prime}}$
$\displaystyle\times$ $\displaystyle
f(l_{1},l_{3};l_{i}^{\prime})f(l_{2},l_{3};l_{i})f(l_{2},l_{4};l_{j}^{\prime})f(l_{1},l_{4};l_{j})$
$\displaystyle\times$ $\displaystyle{\rm sgn}(a_{12})^{l_{1}}{\rm
sgn}(a_{21})^{l_{2}}{\rm sgn}(a_{11})^{l_{3}}{\rm sgn}(a_{22})^{l_{4}}$
$\displaystyle\times$
$\displaystyle\left(\begin{array}[]{ccc}l_{3}&l_{1}&l_{i}^{\prime}\\\
l_{2}&l_{4}&l_{j}^{\prime}\\\
l_{x}&l_{y}&L\end{array}\right)\sum_{\nu_{1},\nu_{2}}(-1)^{\nu_{1}}C_{l_{1}l_{2}}^{\nu_{1}}C_{l_{3}l_{4}}^{\nu_{2}}$
$\displaystyle\times$
$\displaystyle|a_{12}|^{2\nu_{1}+l_{1}+l_{2}}|a_{11}|^{2\nu_{2}+l_{3}+l_{4}},$

where

$\displaystyle C_{jl}^{n}$ $\displaystyle={(2n+j+l+1)!\over
n!(n+j+l+1)!(2(n+j)+1)!!(2(n+l)+1)!!},$

and

$\displaystyle f(a,b;c)=\sqrt{(2a+1)(2b+1)}\langle a0b0|c0\rangle.$ (44)

The notation in the three-by-three bracket is the 9j Clebsch–Gordan
coefficient. The summation is restricted by,

$\displaystyle K$ $\displaystyle=$ $\displaystyle
2n+l_{i}+l_{j}=2n^{\prime}+l_{i}^{\prime}+l_{j}^{\prime}$ (45)
$\displaystyle=$ $\displaystyle 2(\nu_{1}+\nu_{2})+l_{1}+l_{2}+l_{3}+l_{4}.$

For the quantum number of $K=0$ and $L=0$, the Raynal–Revai coefficient
$R^{(ij)}=1$. While for the quantum number of $K=2$ and $L=0$, the
Raynal–Revai coefficient is a two-by-two matrix and can be expressed as (here
we already take the equalities $a_{11}=a_{22}$ and $a_{12}=-a_{21}$, which is
given by the rotation matrix.),

$\displaystyle R^{(ij)}_{2\times
2}=\left(\begin{array}[]{cc}a_{11}^{2}-a_{12}^{2}&-2a_{11}a_{12}\\\
2a_{11}a_{12}&a_{11}^{2}-a_{12}^{2}\end{array}\right).$ (48)

For the quantum number of $K=4$ and $L=0$,

$\displaystyle R^{(ij)}_{3\times
3}=\left(\begin{array}[]{ccc}c_{11}&c_{12}&c_{13}\\\ c_{21}&c_{22}&c_{23}\\\
c_{31}&c_{32}&c_{33}\end{array}\right),$ (52)

with

$\displaystyle c_{11}$ $\displaystyle=$ $\displaystyle-{10\over
3}a_{11}^{2}a_{12}^{2}+a_{12}^{4}+a_{11}^{4},$ $\displaystyle c_{12}$
$\displaystyle=$ $\displaystyle 4\sqrt{2\over 3}a_{11}a_{12}^{3}-4\sqrt{2\over
3}a_{11}^{3}a_{12},$ $\displaystyle c_{13}$ $\displaystyle=$
$\displaystyle{8\sqrt{2}\over 3}a_{11}^{2}a_{12}^{2},$ $\displaystyle c_{22}$
$\displaystyle=$ $\displaystyle a_{12}^{4}+a_{11}^{4}-6a_{11}^{2}a_{12}^{2},$
$\displaystyle c_{23}$ $\displaystyle=$
$\displaystyle{4\over\sqrt{3}}a_{11}a_{12}^{3}-{4\over\sqrt{3}}a_{11}^{3}a_{12},$
$\displaystyle c_{33}$ $\displaystyle=$ $\displaystyle-{2\over
3}a_{11}^{2}a_{12}^{2}+a_{12}^{4}+a_{11}^{4}.$ (53)

and $c_{21}=-c_{12}$, $c_{32}=-c_{23}$, and $c_{31}=c_{13}$.

For the quantum number of $K=6$ and $L=0$,

$\displaystyle R^{(ij)}_{4\times
4}=\left(\begin{array}[]{cccc}c_{11}&c_{12}&c_{13}&c_{14}\\\
c_{21}&c_{22}&c_{23}&c_{24}\\\ c_{31}&c_{32}&c_{33}&c_{34}\\\
c_{41}&c_{42}&c_{43}&c_{44}\end{array}\right),$ (58)

with

$\displaystyle c_{11}$ $\displaystyle=$ $\displaystyle
a_{11}^{6}+7a_{11}^{2}a{12}^{4}-7a_{11}^{4}a_{12}^{2}-a_{12}^{6}$
$\displaystyle c_{12}$ $\displaystyle=$
$\displaystyle-2\sqrt{5}a_{11}a_{12}^{5}-2\sqrt{5}a_{11}^{5}a_{12}+{28\sqrt{5}\over
5}a_{11}^{3}a_{12}^{3},$ $\displaystyle c_{13}$ $\displaystyle=$
$\displaystyle-8a_{11}^{2}a_{12}^{4}+8a_{11}^{4}a_{12}^{2},$ $\displaystyle
c_{14}$ $\displaystyle=$ $\displaystyle-{16\sqrt{5}\over
5}a_{11}^{3}a_{12}^{3},$ $\displaystyle c_{22}$ $\displaystyle=$
$\displaystyle a_{11}^{6}+{67\over 5}a_{11}^{2}a_{12}^{4}-{67\over
5}a_{11}^{4}a_{12}^{2}-a_{12}^{6},$ $\displaystyle c_{23}$ $\displaystyle=$
$\displaystyle{32\over\sqrt{5}}a_{11}^{3}a_{12}^{3}-{8\over\sqrt{5}}a_{11}a_{12}^{5}-{8\over\sqrt{5}}a_{11}^{5}a_{12},$
$\displaystyle c_{24}$ $\displaystyle=$ $\displaystyle-{24\over
5}a_{11}^{2}a_{12}^{4}+{24\over 5}a_{11}^{4}a_{12}^{2},$ $\displaystyle
c_{33}$ $\displaystyle=$ $\displaystyle
a_{11}^{6}-7a_{11}^{4}a_{12}^{2}+7a_{11}^{2}a_{12}^{4}-a_{12}^{6},$
$\displaystyle c_{34}$ $\displaystyle=$
$\displaystyle-{6\over\sqrt{5}}a_{11}^{5}a_{12}+{4\over\sqrt{5}}a_{11}^{3}a_{12}^{3}-{6\over\sqrt{5}}a_{11}a_{12}^{5},$
$\displaystyle c_{44}$ $\displaystyle=$ $\displaystyle a_{11}^{6}-{3\over
5}a_{11}^{4}a_{12}^{2}+{3\over 5}a_{11}^{2}a_{12}^{4}-a_{12}^{6},$ (59)

and $c_{21}=-c_{12}$, $c_{31}=c_{13}$, $c_{41}=-c_{14}$, $c_{32}=-c_{23}$,
$c_{42}=c_{24}$, and $c_{43}=-c_{34}$.

## References

  * Caswell and Lepage [1986] W. E. Caswell and G. P. Lepage, Phys. Lett. B 167, 437 (1986).
  * Brambilla  _et al._ [2000] N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Nucl. Phys. B 566, 275 (2000), arXiv:hep-ph/9907240 .
  * Zhao  _et al._ [2020a] J. Zhao, K. Zhou, S. Chen, and P. Zhuang, Prog. Part. Nucl. Phys. 114, 103801 (2020a), arXiv:2005.08277 [nucl-th] .
  * Silvestre-Brac [1996] B. Silvestre-Brac, Few Body Syst. 20, 1 (1996).
  * Martynenko [2008] A. P. Martynenko, Phys. Lett. B 663, 317 (2008), arXiv:0708.2033 [hep-ph] .
  * Vijande  _et al._ [2015] J. Vijande, A. Valcarce, and H. Garcilazo, Phys. Rev. D 91, 054011 (2015), arXiv:1507.03735 [hep-ph] .
  * Liu  _et al._ [2020] M.-S. Liu, Q.-F. Lü, and X.-H. Zhong, Phys. Rev. D 101, 074031 (2020), arXiv:1912.11805 [hep-ph] .
  * Yang  _et al._ [2020] G. Yang, J. Ping, P. G. Ortega, and J. Segovia, Chin. Phys. C 44, 023102 (2020), arXiv:1904.10166 [hep-ph] .
  * Shi  _et al._ [2020] S. Shi, J. Zhao, and P. Zhuang, Chin. Phys. C 44, 084101 (2020), arXiv:1905.10627 [nucl-th] .
  * Bazavov  _et al._ [2012] A. Bazavov _et al._ , Phys. Rev. D 85, 054503 (2012), arXiv:1111.1710 [hep-lat] .
  * Fodor and Katz [2004] Z. Fodor and S. D. Katz, JHEP 04, 050 (2004), arXiv:hep-lat/0402006 .
  * Xin  _et al._ [2014] X.-y. Xin, S.-x. Qin, and Y.-x. Liu, Phys. Rev. D 90, 076006 (2014), arXiv:2109.09935 [hep-ph] .
  * Li  _et al._ [2019] Z. Li, K. Xu, X. Wang, and M. Huang, Eur. Phys. J. C 79, 245 (2019), arXiv:1801.09215 [hep-ph] .
  * Akiba  _et al._ [2015] Y. Akiba _et al._ , (2015), arXiv:1502.02730 [nucl-ex] .
  * Matsui and Satz [1986] T. Matsui and H. Satz, Phys. Lett. B 178, 416 (1986).
  * Blaizot  _et al._ [2016] J.-P. Blaizot, D. De Boni, P. Faccioli, and G. Garberoglio, Nucl. Phys. A 946, 49 (2016), arXiv:1503.03857 [nucl-th] .
  * Yao and Mehen [2021] X. Yao and T. Mehen, JHEP 02, 062 (2021), arXiv:2009.02408 [hep-ph] .
  * Brambilla  _et al._ [2021] N. Brambilla, M. A. Escobedo, M. Strickland, A. Vairo, P. Vander Griend, and J. H. Weber, JHEP 05, 136 (2021), arXiv:2012.01240 [hep-ph] .
  * Miura  _et al._ [2022] T. Miura, Y. Akamatsu, M. Asakawa, and Y. Kaida, Phys. Rev. D 106, 074001 (2022), arXiv:2205.15551 [nucl-th] .
  * Villar  _et al._ [2023] D. Y. A. Villar, J. Zhao, J. Aichelin, and P. B. Gossiaux, Phys. Rev. C 107, 054913 (2023), arXiv:2206.01308 [nucl-th] .
  * Lafferty and Rothkopf [2020] D. Lafferty and A. Rothkopf, Phys. Rev. D 101, 056010 (2020), arXiv:1906.00035 [hep-ph] .
  * He  _et al._ [2015] H. He, Y. Liu, and P. Zhuang, Phys. Lett. B 746, 59 (2015), arXiv:1409.1009 [hep-ph] .
  * Zhao  _et al._ [2017] J. Zhao, H. He, and P. Zhuang, Phys. Lett. B 771, 349 (2017), arXiv:1603.04524 [nucl-th] .
  * Cho and Lee [2020] S. Cho and S. H. Lee, Phys. Rev. C 101, 024902 (2020), arXiv:1907.12786 [nucl-th] .
  * Andronic  _et al._ [2021] A. Andronic, P. Braun-Munzinger, M. K. Köhler, A. Mazeliauskas, K. Redlich, J. Stachel, and V. Vislavicius, JHEP 07, 035 (2021), arXiv:2104.12754 [hep-ph] .
  * Adamová  _et al._ [2019] D. Adamová _et al._ , (2019), arXiv:1902.01211 [physics.ins-det] .
  * Burnier  _et al._ [2015a] Y. Burnier, O. Kaczmarek, and A. Rothkopf, JHEP 12, 101 (2015a), arXiv:1509.07366 [hep-ph] .
  * Bala  _et al._ [2022] D. Bala, O. Kaczmarek, R. Larsen, S. Mukherjee, G. Parkar, P. Petreczky, A. Rothkopf, and J. H. Weber (HotQCD), Phys. Rev. D 105, 054513 (2022), arXiv:2110.11659 [hep-lat] .
  * Wong  _et al._ [2002] C.-Y. Wong, E. S. Swanson, and T. Barnes, Phys. Rev. C 65, 014903 (2002), [Erratum: Phys.Rev.C 66, 029901 (2002)], arXiv:nucl-th/0106067 .
  * Kawanai and Sasaki [2012] T. Kawanai and S. Sasaki, Phys. Rev. D 85, 091503 (2012), arXiv:1110.0888 [hep-lat] .
  * Krivec [1998] R. Krivec, Few-Body Systems 25, 199 (1998).
  * Barnea  _et al._ [2000] N. Barnea, W. Leidemann, and G. Orlandini, Phys. Rev. C 61, 054001 (2000), arXiv:nucl-th/9910062 .
  * Barnea  _et al._ [2006] N. Barnea, J. Vijande, and A. Valcarce, Phys. Rev. D 73, 054004 (2006), arXiv:hep-ph/0604010 .
  * Marcucci  _et al._ [2020] L. E. Marcucci, J. Dohet-Eraly, L. Girlanda, A. Gnech, A. Kievsky, and M. Viviani, Front. in Phys. 8, 69 (2020), arXiv:1912.09751 [nucl-th] .
  * Raynal and Revai [1970] J. Raynal and J. Revai, Nuovo Cim. A 68, 612 (1970).
  * et al [2020] Z. et al, Progress of Theoretical and Experimental Physics 2020, 083C01 (2020), https://academic.oup.com/ptep/article-pdf/2020/8/083C01/34673722/ptaa104.pdf .
  * Zhao  _et al._ [2020b] J. Zhao, S. Shi, and P. Zhuang, Phys. Rev. D 102, 114001 (2020b), arXiv:2009.10319 [hep-ph] .
  * Crater [1994] H. Crater, Journal of Computational Physics 115, 470 (1994).
  * Padmanath  _et al._ [2014] M. Padmanath, R. G. Edwards, N. Mathur, and M. Peardon, Phys. Rev. D 90, 074504 (2014), arXiv:1307.7022 [hep-lat] .
  * Namekawa  _et al._ [2013] Y. Namekawa _et al._ (PACS-CS), Phys. Rev. D 87, 094512 (2013), arXiv:1301.4743 [hep-lat] .
  * Meinel [2012] S. Meinel, Phys. Rev. D 85, 114510 (2012), arXiv:1202.1312 [hep-lat] .
  * Briceno  _et al._ [2012] R. A. Briceno, H.-W. Lin, and D. R. Bolton, Phys. Rev. D 86, 094504 (2012), arXiv:1207.3536 [hep-lat] .
  * Hasenfratz  _et al._ [1980] P. Hasenfratz, R. R. Horgan, J. Kuti, and J. M. Richard, Phys. Lett. B 94, 401 (1980).
  * Zhang and Huang [2009] J.-R. Zhang and M.-Q. Huang, Phys. Lett. B 674, 28 (2009), arXiv:0902.3297 [hep-ph] .
  * Wang [2012] Z.-G. Wang, Commun. Theor. Phys. 58, 723 (2012), arXiv:1112.2274 [hep-ph] .
  * Migura  _et al._ [2006] S. Migura, D. Merten, B. Metsch, and H.-R. Petry, Eur. Phys. J. A 28, 41 (2006), arXiv:hep-ph/0602153 .
  * Burnier  _et al._ [2015b] Y. Burnier, O. Kaczmarek, and A. Rothkopf, Phys. Rev. Lett. 114, 082001 (2015b), arXiv:1410.2546 [hep-lat] .
  * ALI [2022] (2022), arXiv:2211.04384 [nucl-ex] .
  * Larsen  _et al._ [2019] R. Larsen, S. Meinel, S. Mukherjee, and P. Petreczky, Phys. Rev. D 100, 074506 (2019), arXiv:1908.08437 [hep-lat] .
  * Larsen  _et al._ [2020] R. Larsen, S. Meinel, S. Mukherjee, and P. Petreczky, Phys. Lett. B 800, 135119 (2020), arXiv:1910.07374 [hep-lat] .
  * Laine  _et al._ [2007] M. Laine, O. Philipsen, P. Romatschke, and M. Tassler, JHEP 03, 054 (2007), arXiv:hep-ph/0611300 .
  * Shi  _et al._ [2022] S. Shi, K. Zhou, J. Zhao, S. Mukherjee, and P. Zhuang, Phys. Rev. D 105, 014017 (2022), arXiv:2105.07862 [hep-ph] .