Production of $\Omega NN$ and $\Omega\Omega N$ in ultra-relativistic heavy-ion collisions

Cluster formation in heavy-ion collision can be realized by the coalescence model  Yan et al. (2006); Schaffner-Bielich et al. (2000); Chen et al. (2003); Zhang et al. (2010); Sun and Chen (2017) or other methods like kinetic approaches Sun et al. (2021a); Staudenmaier et al. (2021); He et al. (2020). In this work, a coalescence model constructed by the particle emission distribution and the Wigner density distribution is used to calculate few-body system production in heavy-ion collisions. The multiplicity of three-constituent-cluster is given by,

$\displaystyle\begin{aligned} N_{3b}=g_{3}\int&\prod_{i=1}^{3}\left(d^{4}x_{i}S_{i}\left(x_{i},p_{i}\right)\frac{d^{3}p_{i}}{E_{i}}\right)\times\\ &\ \rho_{3}^{W}\left(x_{1},x_{2},x_{3};p_{1},p_{2},p_{3}\right),\end{aligned}$

(1)

where $\rho^{W}_{3}(x_{1},x_{2},x_{3};p_{1},p_{2},p_{3})$ is the Wigner density function which describes the coalescence probability, and $g_{3}=(2S+1)/\left((2s_{1}+1)(2s_{2}+1)(2s_{3}+1)\right)$ is the coalescence statistical factor Polleri et al. (1999); Zhao et al. (2018); Zhu et al. (2015); Sun et al. (2019, 2021b), $S$ is the total spin for the three-body system and $s_{i}$ is the spin for each constituent particle. Table 1 lists the $g_{3}$ used in this paper for each $\Omega$-hypernucleus ($A$ = 3) and triton.

In this work the particle emission distribution, $S_{i}\left(x_{i},p_{i}\right)$, is given by the blast-wave model Retière and Lisa (2004); Sun and Chen (2017); Zhang et al. (2014); Zhang and Ma (2020) which can describe the particle phase-space distribution in heavy-ion collisions. It assumes that in the rest frame the distribution of momenta is described by either a Bose or Fermi distribution of single particle and then the distribution is boosted into the center-of-mass frame of the total number of particles to describe the probability of finding a particle Cooper and Frye (1974). In heavy-ion collisions, the freeze-out time is considered following a Gaussian distribution Retière and Lisa (2004); Sun and Chen (2017); Zhang et al. (2014); Zhang and Ma (2020); Waqas et al. (2020). The blast-wave model is formalized as,

	$\displaystyle S(x,p)d^{4}x=M_{T}\cosh\left(\eta_{s}-y_{p}\right)$	$\displaystyle f(x,p)\times$		(2)
		$\displaystyle J(\tau)\tau d\tau d\eta_{s}rdrd\varphi_{s},$

where $M_{T}$ and $y_{p}$ are the transverse mass and the rapidity of a single particle, $r$ and $\varphi_{s}$ are the radius and azimuthal angle of coordinate space, $\tau$ and $\eta_{s}$ are proper time and space pseudorapidity. $J(\tau)=\frac{1}{\Delta\tau\sqrt{2\pi}}\exp\left[-\frac{\left(\tau-\tau_{0}\right)^{2}}{2(\Delta\tau)^{2}}\right]$ is the Gaussian distribution of freeze-out proper time, where $\tau_{0}$ and $\Delta\tau$ are the mean value and dispersion of this distribution. $f(x,p)=\frac{2s+1}{(2\pi)^{3}}\left[\exp\left(p^{\mu}u_{\mu}/T_{kin}\right)\pm 1\right]^{-1}$ is the Fermi or Bose distribution of a single particle boosted into the center-of-mass frame, where $s$ is the spin of the particle, $u_{\mu}$ is the four-velocity of a fluid element in the fireball of the particle source and $T_{kin}$ is the freeze-out temperature. The Lorentz invariant can be expressed as,

	$\displaystyle p^{\mu}u_{\mu}=M_{T}\cosh\rho_{\perp}$	$\displaystyle\cosh(\eta_{s}-y_{p})-$		(3)
		$\displaystyle p_{T}\sinh\rho_{\perp}\cos\left(\varphi_{p}-\varphi_{s}\right),$

where $\varphi_{p}$ is the azimuthal angle in momentum space and $\rho_{\perp}$ is the transverse rapidity of fireball with a transverse radius $R_{0}$, defined as $\rho_{\perp}=v\rho_{0_{\perp}}\frac{r}{R_{0}}$. If the parameters ($\tau_{0}$, $\Delta\tau$, $\rho_{0_{\perp}}$, $R_{0}$ and $T_{kin}$) are fixed, the transverse momentum distribution is given as Zhang and Ma (2020):

$$\frac{dN}{2\pi p_{T}dp_{T}dy_{p}}=\int S(x,p)d^{4}x.$$

(4)

In order to obtain the Wigner function, bound state wave functions of $\Omega NN$ and $\Omega\Omega N$ need to be calculated. The non-relativistic Schrödinger equations of $\Omega NN$ and $\Omega\Omega N$’s bound state can be written as,

$$\hat{H}\psi(\bm{x}_{1},\bm{x}_{2},\bm{x}_{3})=E_{b}\psi(\bm{x}_{1},\bm{x}_{2},\bm{x}_{3}),$$

(5)

$$\begin{array}[]{c}\hat{H}=\displaystyle\sum_{i=1}^{3}-\frac{\nabla_{i}^{2}}{2M_{i}}+\sum_{j>i}V_{ij}\left(\bm{r}_{ij}\right),\end{array}$$

(6)

where $\displaystyle\psi(\bm{x}_{1},\bm{x}_{2},\bm{x}_{3})=\sum_{i=1}^{3}\psi_{i}\left(\bm{x}_{i}\right)$ is total wave function of three-body system, $\psi_{i}(\bm{x}_{i})$ and $\bm{x}_{i}$ are the wave function and coordinate of $i$-th particle, respectively, $\bm{r}_{ij}$ is relative coordinate between $i$-th and $j$-th particle defined as $\bm{r}_{ij}=\bm{x}_{i}-\bm{x}_{j}$. The potentials between $\Omega$-$N$ and $\Omega$-$\Omega$ are the fit results from the HAL-QCD simulation Iritani et al. (2019); Gongyo et al. (2018), and the $N$-$N$ potential is taken as the Malfliet-Tjon potential Malfliet and Tjon (1969):

		$\displaystyle V_{NN}(r)=\sum_{i=1}^{2}C_{i}\frac{e^{-\mu_{i}r}}{r},$		(7)
		$\displaystyle V_{N\Omega}(r)=b_{1}e^{-b_{2}r^{2}}+b_{3}\left(1-e^{-b_{4}r^{2}}\right)\left(\frac{e^{-M_{\pi}r}}{r}\right)^{2},$
		$\displaystyle V_{\Omega\Omega}(r)=\sum_{i=1}^{3}C_{i}e^{-\left(r/d_{i}\right)^{2}},$

where $M_{\pi}$ is taken as $146$ MeV (near the physical mass $140$ MeV). The parameters are listed in Table 2.

There are many methods to solve this kind of three-body equations, such as the Faddeev equation Faddeev (1961); Thompson et al. (2004); Kovalchuk (2014) and the variation method. One kind of the variation methods is mainly based on the hyperspherical-harmonics (HH) method Kievsky et al. (1993, 1994); Kievsky (1997); Kievsky et al. (1997), in which the coordinates are transformed into center-of-mass frame by using the Jacobi transform,

$$\begin{array}[]{c}\left(\begin{array}[]{c}\bm{R}\\ \bm{r}_{1}\\ \bm{r}_{2}\end{array}\right)=J\cdot\left(\begin{array}[]{c}\bm{x}_{1}\\ \bm{x}_{2}\\ \bm{x}_{3},\end{array}\right)\\ \left(\begin{array}[]{l}\bm{P}\\ \bm{q}_{1}\\ \bm{q}_{2}\end{array}\right)=\left(J^{-1}\right)^{T}\cdot\left(\begin{array}[]{l}\bm{p}_{1}\\ \bm{p}_{2}\\ \bm{p}_{3},\end{array}\right)\end{array}$$

(8)

where $J$ is the Jacobi matrix, it reads

$$J=\begin{pmatrix}\frac{M_{1}}{M_{tot}}&\frac{M_{2}}{M_{tot}}&\frac{M_{3}}{M_{tot}}\\[5.0pt] 0&-\sqrt{\frac{M_{2}M_{3}}{M_{23}Q}}&\sqrt{\frac{M_{2}M_{3}}{M_{23}Q}}\\[5.0pt] -\sqrt{\frac{M_{1}M_{23}}{M_{tot}Q}}&\sqrt{\frac{M_{2}^{2}M_{1}}{M_{tot}M_{23}Q}}&\sqrt{\frac{M_{3}^{2}M_{1}}{M_{tot}M_{23}Q}}\end{pmatrix},$$

(9)

where $M_{i}$ is the mass of $i$th particle, $M_{tot}=M_{1}+M_{2}+M_{3}$ is the total mass, $M_{23}=M_{2}+M_{3}$ is the total mass of particles 2 and 3, $Q=\sqrt{(M_{1}M_{2}M_{3})/M_{tot}}$ is the reduced mass which normalizes the Jacobi matrix. For simplicity, the indexes of particles are chosen as symmetric as possible. In this article, the particles 2 and 3 prefer to be identical and particle 1 is different for a three-body nucleus. Sequentially the three-body Schrödinger equation separates into the center of mass motion (no effect on binding energy and relative wave function) and the relative motion Chattopadhyay et al. (1996); Khan (2012),

$$\widehat{T}\psi(\vec{r})+\sum_{j>i}V_{ij}\left(\bm{r}_{ij}\right)\psi(\vec{r})=E_{b}\psi(\vec{r}),$$

(10)

where $\vec{r}=\left(\bm{r}_{1},\bm{r}_{2}\right)=\left(\rho,\alpha,\theta_{1},\phi_{1},\theta_{2},\phi_{2}\right)$ is defined in a six-dimensional hypersphere coordinate, $\rho=\sqrt{r_{1}^{2}+r_{2}^{2}}$ is the hyperradius, $\alpha=\arctan\left(r_{2}/r_{1}\right)$ is the hyperpolar angle which ranges from 0 to $\pi/2$ Chattopadhyay et al. (1996); Khan (2012); He et al. (2015), $\theta_{i},\phi_{i}$ are the azimuth angles of $\bm{r}_{i}$, and the volume element is $d^{6}\vec{r}=\rho^{5}\sin^{2}{\alpha}\cos^{2}{\alpha}\sin{\theta_{1}}\sin{\theta_{2}}d\rho d\alpha d\theta_{1}d\phi_{1}d\theta_{2}d\phi_{2}$. The momentum and angular momentum operators are defined as Kievsky et al. (1993, 1994); Chattopadhyay et al. (1996); Khan (2012); He et al. (2015),

$$\hat{T}=\frac{1}{2Q}\left(-\frac{\partial^{2}}{\partial\rho^{2}}-\frac{5}{\rho}\frac{\partial}{\partial\rho}+\frac{\hat{L}^{2}}{\rho^{2}}\right),$$

(11)

where

$$\hat{L}^{2}=-\frac{\partial^{2}}{\partial\alpha^{2}}-4\cot 2\alpha\frac{\partial}{\partial\alpha}+\frac{\hat{l}_{1}^{2}}{\cos^{2}\alpha}+\frac{\hat{l}_{2}^{2}}{\sin^{2}\alpha}.$$

(12)

The eigen function of $\hat{L}^{2}$ is a hyperspherical harmonic function Kievsky et al. (1993, 1994); Kievsky (1997); Kievsky et al. (1997):

$\displaystyle\begin{aligned} \mathcal{Y}_{K,\kappa}(\alpha,\theta_{1},\phi_{1},\theta_{2},\phi_{2})=N_{kl_{1}l_{2}}\cos(\alpha)^{l_{1}}\sin(\alpha)^{l_{2}}P_{k}^{l_{2}+1/2,\ l_{1}+1/2}(\cos(2\alpha))\times\\ \left\{\left\{Y_{l_{1}}\left(\theta_{1},\phi_{1}\right)Y_{l_{2}}\left(\theta_{2},\phi_{2}\right)\right\}_{L}\left\{s_{i}s_{jk}\right\}_{S_{a}}\right\}_{JJ_{z}}\left\{t_{i}t_{jk}\right\}_{TT_{z}},\end{aligned}$

(13)

and

$\displaystyle\hat{L}^{2}\mathcal{Y}_{K,\kappa}=(K+4)K\ \mathcal{Y}_{K,\kappa},$

(14)

where $K=2k+l_{1}+l_{2}$ is the total hyperangular momentum number, $q$ is a nonnegative integer, $l_{i}$ and $m_{i}$ is the orbital angular momentum number of $\bm{r}_{i}$ direction, $\kappa$ represents the $L$-spin-isospin state defined as $\kappa=\{JJ_{z}(L(l_{1}l_{2})S_{a}(s_{i}s_{jk}))TT_{z}(t_{i}t_{jk})\}$, $N_{kl_{1}l_{2}}$ is a normalization factor Kovalchuk (2014),

$$N_{kl_{1}l_{2}}=\sqrt{\frac{2\ k!\ (K+2)\ (k+l_{1}+l_{2}+1)!}{\Gamma(k+l_{1}+3/2)\ \Gamma(k+l_{2}+3/2))}}$$

(15)

and $P_{k}^{a,b}(x)$ is the Jacobi polynomial and $Y_{m}^{l}(\theta,\phi)$ is the Spherical Harmonic function. The orthogonal basis radial function can be chosen as

$$u_{n}^{[\lambda]}(\rho)=\sqrt{\left(\frac{2\lambda}{n}\right)^{3}\frac{(n-2)!}{2n\ (n+1)!}}e^{-\lambda\rho/n}\left(\frac{2\lambda\rho}{n}\right)L_{n-2}^{3}\left(\frac{2\lambda\rho}{n}\right)\rho^{-\frac{5}{2}}\ (n\geq 2),$$

(16)

in which $\lambda$ is a variation parameter, $n$ is the radial basis index, $L_{a}^{b}(\rho)$ is the associated Laguerre polynomial. Then the orthogonal basis function can be constructed as,

$$\left\langle\vec{r}\ |n,K,\kappa\right\rangle=u_{n}^{[\lambda]}(\rho)\mathcal{Y}_{K,\kappa}(\alpha,\theta_{1},\phi_{1},\theta_{2},\phi_{2}).$$

(17)

Then the relative motion Hamiltonian $\hat{H}$ can be expanded into matrix form $\left\langle n,K,\kappa\left|\hat{H}\right|n^{\prime},K^{\prime},\kappa^{\prime}\right\rangle$. The following assumptions are taken to reduce the dimensions of the matrix: 1) assume that the nucleus is spherical by setting the total $L$ = 0, corresponding to the ground state; 2) the $(I)J^{P}$ is fixed as the same as Garcilazo and Valcarce used Garcilazo and Valcarce (2019); 3) if particle 2 and 3 are identical, the parity between them must be odd Kievsky et al. (1993); 4) $l_{1}$, $l_{2}\leq 6$ is enough for required precision Kievsky (1997), and the number $n$ ranges from 2 to 11 and $k$ up to 45. The matrix elements have been calculated numerically by a Laguerre-Gauss quadrature for the integrals in the hyperradius $\rho$ and a Legendre-Gauss quadrature for the hyperangle $\alpha$  Abramowitz and Stegun .

But the elements of Hamiltonian matrix need six dimensional integral and the complex expressions of $\bm{r}_{12}$ and $\bm{r}_{31}$ for the hypersphere coordinate are based on the transforms of (8) and (9), where $\bm{r}_{2}=-\sqrt{\frac{M_{2}M_{3}}{\left(M_{2}+M_{3}\right)Q}}\ \bm{r}_{23}$. In order to simplify the calculation, Raynal and Revai Raynal and Revai (1970) put forward the RR coefficient which is similar to Clebsch-Gordan coefficient. For example as shown in Fig. 1, it is convenient to calculate $V_{23}(\bm{r}_{23})$ when the hypersphere is based on $\bm{r}_{1}^{I}$ and $\bm{r}_{2}^{I}$ in Coordinate I but hard to calculate $V_{12}(\bm{r}_{12})$ and $V_{31}(\bm{r}_{31})$ for the complex expressions of $\bm{r}_{12}$ and $\bm{r}_{31}$. By using RR coefficient, the hyperspherical harmonic function $\left|I;n,K,\kappa\right\rangle$, defined in the coordinate I, can be expanded by $\left|II(III);n,K,\kappa^{\prime}\right\rangle$ in coordinate II (III),

$$\left|I;n,K,\kappa\right\rangle=\sum_{\kappa_{k}}\left\langle l^{I}_{1}l^{I}_{2}|l^{j}_{1}l^{j}_{2}\right\rangle_{K,L}\left\langle s_{1}s_{23};S|s_{1_{j}}s_{23_{j}};S\right\rangle\left\langle t_{1}t_{23};T|t_{1_{j}}t_{23_{j}};T\right\rangle\left|j;n,K,\kappa_{j}\right\rangle,$$

(18)

where $\left\langle l^{I}_{1}l^{I}_{2}|l^{j}_{1}l^{j}_{2}\right\rangle_{K,L}$ is the RR coefficient which requires that $K$ and $L$ are same in transformation, $\left\langle s_{1}s_{23};S|s_{1_{j}}s_{23_{j}};S\right\rangle$ and $\left\langle t_{1}t_{23};T|t_{1_{j}}t_{23_{j}};T\right\rangle$ are Clebsch-Gordan coefficients, $j$ represents the coordinate II or III, $1_{j}$ and $23_{j}$ represent the particle 2 (3) and the pair of particle 3 (1) and 1 (2) when $j=II(III)$. It is clear that the definition of $\rho$ is same in all coordinates, so the index $n$ does not need to change in the transform (18). After the transformation, $\bm{r}_{23_{j}}$ only relates to the $\rho$ and $\alpha_{j}$, which means the six dimensional integral is simplified into a double integral and a sum of $\kappa_{j}$.

After the calculation of Hamiltonian matrix, it is natural to calculate the minimum eigenvalue of the matrix as the binding energy $B[\lambda]$ of the three-body system and the corresponding eigenvector is the list of coefficients for the basis functions. And the binding energy $B[\lambda]$ requires $\delta B[\lambda]/\delta\lambda=0$, which means that the binding energy is also the minimum point of variation parameter $\lambda$.

Garcilazo and Valcarce Garcilazo and Valcarce (2016) solved three-body amplitudes by the Faddeev equations Faddeev (1961) with considering the spin and isospin freedom. They assumed that three particles were in $S$-wave by which the spin-isospin state was constructed and two-body amplitudes with the Legendre polynomials were expanded to solve the Faddeev equations.

Table 3 shows the calculated binding energy of ${}^{3}H$, ${}_{\Lambda}^{3}H$ and $\Omega pn$ and the comparison with other theoretical results as well as experimental results. The potential between $N$ and $\Lambda$ used in ${}_{\Lambda}^{3}H$ binding energy calculation is YNG-ND interactions Yamamoto et al. (1994); Hiyama et al. (1997) with $k_{F}=0.84fm^{-1}$Hiyama et al. (2002). It can be seen that this calculation of $\Omega pn$ is consistent with the results from Garcilazo and Valcarce’s results  Garcilazo and Valcarce (2019). The error of $pn\Omega$ binding energy is estimated from the fitting errors of the $N-\Omega$ potential. The results of ${}^{3}H$ and ${}_{\Lambda}^{3}H$ are close to experimental results Davis (1986); Adam et al. (2020) and theoretical calculations as well Egorov and Postnikov (2021). Like ${}_{\Lambda}^{3}H$ consisting of spin $\frac{3}{2}$ and $\frac{1}{2}$, one is the ground state (spin $\frac{1}{2}$) and one is thought as a virtual state (spin $\frac{3}{2}$) near the $\Lambda d$ threshold Schäfer et al. (2022), $\Omega pn$ can also be mix of spin $\frac{5}{2}$, $\frac{3}{2}$ and $\frac{1}{2}$. According to the HAL QCD’s calculation Gongyo et al. (2018), the ${}^{3}S_{1}$ $\Omega N$ interaction is too weak to form a bound $\Omega N$ with spin 1. So the ratio of lower spin in $\Omega pn$ is small. In this paper $\Omega pn$ is considered as spin $\frac{5}{2}$.

There is another method for few-body system interaction by which the system is not deeply bound, called as the folding model Watanabe (1958); Etminan and Firoozabadi (2019). The folding model assumes that nucleus is bound as a molecular state like dibaryon-baryon state. For the ${}^{3}S_{1}\ \Omega N$ interaction, it is not as strong as ${}^{5}S_{2}\ \Omega N$, the folding model can be applied in $\Omega nn$, $\Omega pp$, $\Omega\Omega n$ and $\Omega\Omega p$, which softens the ${}^{5}S_{2}\ \Omega N$ interaction. The model uses the free dibaryon wave function to average the potential between dibaryon and baryon,

		$\displaystyle U_{F}(\mathbf{R}_{F})=\int d^{3}\mathbf{r}_{d}\ \psi_{d}^{*}\left(\mathbf{r}_{d}\right)\psi_{d}\left(\mathbf{r}_{d}\right)\times$		(19)
		$\displaystyle\left[V_{12}\left(\mathbf{R}_{F}-\frac{M_{3}\ \mathbf{r}_{d}}{M_{2}+M_{3}}\right)+V_{13}\left(\mathbf{R}_{F}+\frac{M_{2}\ \mathbf{r}_{d}}{M_{2}+M_{3}}\right)\right],$

where $\psi_{d}$ is the dibaryon wave function which is consisted of particle 2 and 3, $U_{F}(\mathbf{R}_{F})$ is average potential and the $\mathbf{R}_{F}$ is relative coordinate between the dibaryon and baryon. This method also simplifies the three-body bound state into two two-body bound states (dibaryon and dibaryon-baryon). The total wave function is $\Psi(\mathbf{R}_{F},\mathbf{r}_{d})=\psi_{d}(\mathbf{r})\ \psi_{mole}(\mathbf{R}_{F})$ and total binding energy is $E=E_{d}+E_{mole}$, where $\psi_{mole}(\bm{R}_{F})$ is the molecular state wave function calculated with the average potential $U_{F}(\bm{R}_{F})$ and $E_{mole}$ is the binding energy of molecular state. The binding energies of $\Omega nn$, $\Omega pp$, $\Omega\Omega n$ and $\Omega\Omega p$ are calculated by the folding model and their errors are estimated from the fitting error of $N-\Omega$ and $\Omega-\Omega$ potential. The results are listed in Table 4. It can be found that different combinations of dibaryon in three-body systems result in different binding energies which are corresponding to different decay channels and will be discussed later.

The Wigner function introduced in Eq. (1) is written as Chen et al. (2003); Zhang et al. (2010); Sun and Chen (2017),

	$\displaystyle\rho^{W}(\vec{r},\vec{q})=\int\psi\left(\vec{r}+\frac{\vec{R}}{2}\right)$	$\displaystyle\psi^{*}\left(\vec{r}-\frac{\vec{R}}{2}\right)\times$		(20)
		$\displaystyle\exp(-i\vec{q}\cdot\vec{R})d^{6}\vec{R},$

where $\vec{r}=(\bm{r}_{1},\bm{r}_{2})$, $\vec{q}=(\bm{q}_{1},\bm{q}_{2})$ are the relative coordinate and momentum, and $\psi(\vec{x})$ is the relative wave function. For the three-body system it is expressed in six dimensions, the Wigner function will be 12 dimensions, which is impossible to draw a picture and hardly calculated. After performing the calculation of eigenvector of Hamiltonian matrix, the major contribution of total wave function comes from a few bases which contribute more than $94\%$ to total amplitude for the parameters of them are large (larger than 0.08). With considering the fitting errors of potential, the total relative errors of such simplified wave functions are about $10\%$, So this kind of simplification retains most information of origin wave function. If the selected bases are only radial related, the total wave function can be simplified as the sum of these bases with weights of their parameters. And then the simplified wave function is only radial related. The Wigner function can be simplified as,

	$\displaystyle\rho_{3}^{W}(r,q,\theta)=$	$\displaystyle\int\psi\left(\sqrt{r^{2}+\frac{R^{2}}{4}+rR\cos\theta_{1}}\right)\times$		(21)
		$\displaystyle\ \ \ \ \ \psi^{*}\left(\sqrt{r^{2}+\frac{R^{2}}{4}-rR\cos\theta_{1}}\right)\times$
		$\displaystyle\ \ \ \ \ \ \ \exp\left(-iqR\cos\theta\cos\theta_{1}\right)\times$
		$\displaystyle\ \ \ \ \ \ \ \ \exp\left(-iqR\sin\theta\sin\theta_{1}\cos\theta_{2}\right)\times$
		$\displaystyle\ \ \ \ \ \ \ \ \ \ 2\pi^{2}R^{5}\sin^{4}\theta_{1}\sin^{3}\theta_{2}dRd\theta_{1}d\theta_{2}.$

A Laguerre-Gauss quadrature is applied for the integrals of hyperradius $R$ and $\theta_{1},\ \theta_{2}$ is integrated by a Legendre-Gauss quadrature Abramowitz and Stegun . The coordinate is defined in a six-dimensional spherical coordinate as $\vec{R}=\left(R,\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5}\right)$, which can be transformed into the six-dimensional Cartesian coordinate:

	$\displaystyle\vec{R}=$	$\displaystyle(\bm{R}_{1},\bm{R}_{2})=(R_{1x},R_{1y},R_{1z},R_{2x},R_{2y},R_{2z})=$		(22)
		$\displaystyle\ (R\cos{\theta_{1}},R\sin{\theta_{1}}\cos{\theta_{2}},R\sin{\theta_{1}}\sin{\theta_{2}}\cos{\theta_{3}},$
		$\displaystyle\ \ R\sin{\theta_{1}}\sin{\theta_{2}}\sin{\theta_{3}}\cos{\theta_{4}},$
		$\displaystyle\ \ R\sin{\theta_{1}}\sin{\theta_{2}}\sin{\theta_{3}}\sin{\theta_{4}}\cos{\theta_{5}},$
		$\displaystyle\ \ R\sin{\theta_{1}}\sin{\theta_{2}}\sin{\theta_{3}}\sin{\theta_{4}}\sin{\theta_{5}}),$

and $0\leq\theta_{i}\leq\pi$ ($i=1,2,3,4$), $0\leq\theta_{5}\leq 2\pi$, the volume element is $d^{6}\vec{R}=R^{5}dR\prod_{i=1}^{5}\sin^{5-i}{\theta_{i}}d\theta_{i}$. The $\vec{r}$ in Eq. (21) is set at $\left(r,0,0,0,0,0\right)$, the $\vec{q}$ is set at $(q,\theta,0,0,0,0)$. By integrating out the angle, the probability to find the $pn\Omega$ ground bound state can be obtained at six-dimensional hyperspherical radius $r$ and at six-dimensional hyperspherical momentum $q$ He et al. (2015),

$$P(r,q)=\frac{1}{24\pi}r^{5}q^{5}\int_{0}^{\pi}\rho_{3}^{W}(r,q,\theta)\sin^{4}\theta d\theta,$$

(23)

which is shown in Fig. 2. The Wigner probability is similar to a Gaussian distribution with tails in both coordinate and momentum space. The most probable position in the coordinate-momentum phase space is located at $(r,q)\sim(2\text{ fm},200\text{ MeV})$. And the normalization of the probability,

$$\int^{\infty}_{0}P(r,q)drdq=1.$$

(24)

If the wave function relates to not only $\rho$ but also $\alpha$, in other word, the wave function relates to both $r_{1}$ and $r_{2}$ which are defined in Fig. 1. Wigner transformation is more complex. $\psi(\rho,\alpha)=\left.\left\langle\rho\alpha\right|\psi\right\rangle$ can be simplified into $\sum_{n_{1},n_{2}}\left.\left\langle r_{1}r_{2}\right|n_{1}n_{2}\right\rangle\left\langle n_{1}n_{2}\left|\psi\right\rangle\right.$. $\left\langle r_{i}\left|n_{i}\right\rangle\right.$ is a 3-dimension radial orthogonal basis which is the same as (16) but the last term is $r_{i}^{-1/2}$ with the same variation parameter $\lambda$ for different $n_{i}$. Here $n_{i}$ ranges from 2 to 26 with $\lambda=10000$. By this way, Wigner transformation can be rewritten as:

	$\displaystyle\rho^{W}(\vec{r},\vec{q})$	$\displaystyle=\int d^{6}\vec{R}\left\langle\psi\left|\vec{r}-\frac{\vec{R}}{2}\right\rangle\right.\left.\left\langle\vec{r}+\frac{\vec{R}}{2}\right|\psi\right\rangle e^{-i\vec{q}\cdot\vec{R}}$		(25)
		$\displaystyle=\sum_{n_{1},n_{2},n^{\prime}_{1},n^{\prime}_{2}}\left.\left\langle\psi\right|n_{1}n_{2}\right\rangle\left\langle n^{\prime}_{1}n^{\prime}_{2}\left|\psi\right\rangle\right.\times$
		$\displaystyle\prod_{i=1,2}\int d^{3}\mathbf{R_{i}}\left\langle n_{i}\left|\mathbf{r}_{i}-\frac{\mathbf{R}_{i}}{2}\right\rangle\right.\left.\left\langle\mathbf{r}_{i}+\frac{\mathbf{R}_{i}}{2}\right|n^{\prime}_{i}\right\rangle e^{-i\mathbf{q}_{i}\cdot\mathbf{R}_{i}}.$

A complex Wigner transformation is simplified by a series of three-dimension Wigner transform.

For the folding model, $\rho_{3}^{W}=\rho_{di}^{W}\times\rho_{di-b}^{W}$, where $\rho_{di}^{W}$ is the Wigner density function for dibaryon and $\rho_{di-b}^{W}$ is the Wigner density function for the pair of dibaryon and third baryon. Both of these two Wigner density functions can be calculated as did in our previous work Zhang and Ma (2020) for two-body systems.

The main errors of Wigner function are from the errors of wave functions. From the relationship between Wigner function and the wave function, the errors of Wigner function are estimated to be about $20\%$.