Production of (cid:2) N N and (cid:2)(cid:2) N in ultra-relativistic heavy-ion collisions

Even though lots of (cid:2) -hypernuclei have been found and measured, multi-strangeness hypernuclei consisting of (cid:3) are not yet discovered. The studies of multi-strangeness hypernuclei help us further understand the interaction between hyperons and nucleons. Recently the (cid:3) N and (cid:3)(cid:3) interactions as well as binding energies were calculated by the HAL-QCD’s lattice Quantum Chromo-Dynamics (LQCD) simulations and production rates of (cid:3) -dibaryon in Au + Au collisions at RHIC and Pb + Pb collisions at LHC energies were estimated by a coalescence model. The present work discusses the production of more exotic triple-baryons including (cid:3) , namely (cid:3) N N and (cid:3)(cid:3) N as well as their decay channels. A variation method is used in calculations of bound states and binding energy of (cid:3) N N and (cid:3)(cid:3) N with the potentials from the HAL-QCD’s results. The productions of (cid:3) N N and (cid:3)(cid:3) N are predicted by using a blast-wave model plus coalescence model in ultra-relativistic heavy-ion collisions at √ s NN = 200 GeV and 2 . 76 TeV. Furthermore, plots for baryon number dependent yields of different baryons ( N and (cid:3) ), their dibaryons and hypernuclei are made and the production rate of a more exotic tetra-baryon ( (cid:3)(cid:3) N N ) is extrapo-lated.


Introduction
Hypernucleus consisting of hyperons and nucleons is described by not only mass and charge but also hypercharge. Danysz and Pniewski first discovered the 3 H from cosmic rays in 1952 [1]. Since then more attention has been paid to hypernuleus research and many -hypernuclei were discovered in cosmic rays as well as by accelerator beams [2,3]. Recently a e-mail: song_zhang@fudan.edu.cn b e-mail: mayugang@fudan.edu.cn (corresponding author) the observation of − -14 N was also reported by the J-PARC laboratory [4]. Nowadays, relativistic heavy-ion collisions can produce a large number of strange hyperons [5][6][7][8], which provides a venue to discover the hypernucleus even antihypernucleus. The research on hypernuclei is becoming an important direction in heavy-ion collision experiments [9]. On the other hand, multi-quark exotic hadrons or hadronic molecules are also in current focus in particle and heavy-ion physics [10][11][12][13][14][15][16]. The HAL-QCD Collaboration reported the most strangeness dibaryon candidates, N and [17,18] by the Lattice Quantum Chromo-Dynamics (LQCD) simulations. Based on their results, our previous work calculated the production of and N dibaryons and gave the yields of -dibaryon by the blast-wave model or A Multiphase Transport (AMPT) model coupling with a coalescence model in relativistic heavy-ion collisions at √ s N N = 200 GeV and 2.76 TeV [19].
The attractive nature of the N interaction leads to the possible existence of an N dibaryon with strangeness = −3, spin = 2, and isospin = 1/2, which was first proposed in Ref. [20]. Later on the HAL-QCD Collaboration calculated the -N and -interaction by the LQCD simulations near the physical point and the LQCD potentials are fitted by Gaussians and (Yukawa) 2 . The results lead to the binding energy B n = 1.54 MeV, B p = 2.46 MeV and B = 1.6 MeV [17,18]. The STAR Collaboration made a first measurement of momentum correlation functions of p − for Au + Au collisions at √ s N N = 200 GeV [21] which indicates that the scattering length is positive for the proton-interaction and favors the proton-bound state hypothesis by comparing with the predictions based on the proton-interaction extracted from (2 + 1)-flavor LQCD simulations [22]. Later on the ALICE collaboration measured the momentum correlation function of p − in pp collision at √ s = 13 TeV [23] which supports the HAL-QCD result [17]. The poten-tials given by the HAL-QCD [17] are also used to calculate the binding energy of -hypernuclei with A = 3. Garcilazo and Valcarce [24] calculated the bound states of threebody -hypernuclei, namely N N and N , by solving the Faddeev equations [25] with the HAL-QCD potentials and obtained their binding energies ranging from 2 MeV to 20 MeV.
In this work the productions of N N and N are calculated by a coalescence model in which the nucleon and hyperon phase space distributions are given by a blast-wave model [19,[26][27][28]. The potentials from the LQCD are taken into account to obtain the relative wave functions and binding energies of N N and N by solving Schrödinger equation via a variation method. The estimation of the yields of N N and N will shed light on searching for -hypernuclei in experiment, such as at LHC-ALICE.
The calculation of production is introduced in Sect. 2, which includes the brief introductions of the blast-wave model and the coalescence model [19,[26][27][28], simplification of Wigner function as well as the variation calculation method of three-body bound state. It is compared with the results from the Faddeev equations used by Garcilazo and Valcarce [24]. In Sect. 3, productions of N N and N are reported for Au + Au collisions at √ s N N = 200 GeV and Pb + Pb collisions at √ s N N = 2.76 TeV. The decay channels of N N and N are also discussed in this section.

Blast-wave model and coalescence model
Cluster formation in heavy-ion collision can be realized by the coalescence model [27,[29][30][31][32] or other methods like kinetic approaches [33][34][35]. 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, is the Wigner density function which describes the coalescence probability, and g 3 = (2S + 1)/ ((2s 1 + 1)(2s 2 + 1)(2s 3 + 1)) is the coalescence statistical factor [36][37][38][39][40], 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 -hypernucleus (A = 3) and triton. In this work the particle emission distribution, S i (x i , p i ), is given by the blast-wave model [19,[26][27][28] 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 [41]. In heavyion collisions, the freeze-out time is considered following a Gaussian distribution [19,[26][27][28]42]. The blast-wave model is formalized as, where M T and y p are the transverse mass and the rapidity of a single particle, r and ϕ s are the radius and azimuthal angle of coordinate space, τ and η s are proper time and space pseudorapidity. 2 is the Gaussian distribution of freeze-out proper time, where τ 0 and τ are the mean value and dispersion of this distribution.
Bose distribution of a single particle boosted into the centerof-mass frame, where s is the spin of the particle, u μ 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, where ϕ p is the azimuthal angle in momentum space and ρ ⊥ is the transverse rapidity of fireball with a transverse radius R 0 , defined as ρ ⊥ = vρ 0 ⊥ r R 0 . If the parameters (τ 0 , τ , ρ 0 ⊥ , R 0 and T kin ) are fixed, the transverse momentum distribution is given as [19]:

Solving three-body bound state
In order to obtain the Wigner function, bound state wave functions of N N and N need to be calculated. The non-relativistic Schrödinger equations of N N and N 's bound state can be written as, where ψ(x 1 , is total wave function of three-body system, ψ i (x i ) and x i are the wave function and coordinate of i-th particle, respectively, r i j is relative coordinate between ith and jth particle defined as r i j = x i − x j . The potentials between -N and -are the fit results from the HAL-QCD simulation [17,18], and the N -N potential is taken as the Malfliet-Tjon potential [43]: where M π 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 threebody equations, such as the Faddeev equation [25,44,45] and the variation method. One kind of the variation methods is mainly based on the hyperspherical-harmonics (HH) method [46][47][48][49], in which the coordinates are transformed into center-of-mass frame by using the Jacobi transform, ⎛ where M i is the mass of ith particle, 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 [50,51], where r = (r 1 , is the hyperradius, α = arctan (r 2 /r 1 ) is the hyperpolar angle which ranges from 0 to π/2 [50][51][52], θ i , φ i are the azimuth angles of r i , and the volume element is d 6 r = ρ 5 sin 2 α cos 2 α sin θ 1 sin θ 2 dρdαdθ 1 dφ 1 dθ 2 dφ 2 .
The momentum and angular momentum operators are defined as [46,47,[50][51][52], wherê The eigen function ofL 2 is a hyperspherical harmonic function [46][47][48][49]: and 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 r i direction, κ represents the L-spin-isospin state defined as κ = {J J z (L(l 1 l 2 )S a (s i s jk ))T T z (t i t jk )}, N kl 1 l 2 is a normalization factor [45], and P a,b k (x) is the Jacobi polynomial and Y l m (θ, φ) is the Spherical Harmonic function. The orthogonal basis radial function can be chosen as in which λ is a variation parameter, n is the radial basis index, L b a (ρ) is the associated Laguerre polynomial. Then Table 2 Parameters of potentials V N N (r ) [43], V N (r ), V (r ) [17,18] (5) 0.305 (29) 0.949 (58) the orthogonal basis function can be constructed as, Then the relative motion HamiltonianĤ can be expanded into matrix form n, K , κ Ĥ n , K , κ . 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 [24]; (3) if particle 2 and 3 are identical, the parity between them must be odd [46]; (4) l 1 , l 2 ≤ 6 is enough for required precision [48], 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 ρ and a Legendre-Gauss quadrature for the hyperangle α [53].
But the elements of Hamiltonian matrix need six dimensional integral and the complex expressions of r 12 and r 31 for the hypersphere coordinate are based on the transforms of (8) and (9), where r 2 23 . In order to simplify the calculation, Raynal and Revai [54] 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 (r 23 ) when the hypersphere is based on r I 1 and r I 2 in Coordinate I but hard to calculate V 12 (r 12 ) and V 31 (r 31 ) for the complex expressions of r 12 and r 31 . By using RR coefficient, the hyperspherical harmonic function |I ; n, K , κ , defined in the coordinate I, can be expanded by I I (I I I ); n, K , κ in coordinate II (III), where l I is the RR coefficient which requires that K and L are same in transformation, s 1 s 23 ; S|s 1 j s 23 j ; S and t 1 t 23 ; T |t 1 j t 23 j ; T 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 It is clear that the definition of ρ is same in all coordinates, so the index n does not need to change in the transform (18). After the transformation, r 23 j only relates to the ρ and α j , which means the six dimensional integral is simplified into a double integral and a sum of κ j . After the calculation of Hamiltonian matrix, it is natural to calculate the minimum eigenvalue of the matrix as the binding energy B[λ] of the three-body system and the corresponding eigenvector is the list of coefficients for the basis functions. And the binding energy B[λ] requires δ B[λ]/δλ = 0, which means that the binding energy is also the minimum point of variation parameter λ.
Garcilazo and Valcarce [55] solved three-body amplitudes by the Faddeev equations [25] 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 , 3 H and pn and the comparison with other theoretical results as well as experimental results. The potential between N and used in 3 H binding energy calculation is YNG-ND interactions [56,57] with k F = 0.84 fm −1 [58]. It can be seen that this calculation of pn is consistent with the results from Garcilazo and Valcarce's results [24]. The error of pn binding energy is estimated from the fitting errors of the N − potential. The results of 3 H and 3 H are close to experimental results [59,60] and theoretical calculations as well [61]. Like 3 H consisting of spin 3 2 and 1 2 , one is the ground state (spin 1 2 ) and one is thought as a virtual state (spin 3 2 ) near the d threshold [62], pn can also be mix of spin 5 2 , 3 2 and 1 2 . According to the HAL QCD's calculation [18], the 3 S 1 N interaction is too weak to form a bound N with spin 1. So the ratio of lower spin in pn is small. In this paper pn is considered as spin 5 2 . There is another method for few-body system interaction by which the system is not deeply bound, called as the folding model [64,65]. The folding model assumes that nucleus is bound as a molecular state like dibaryon-baryon state. For the 3 S 1 N interaction, it is not as strong as 5 S 2 N , the   folding model can be applied in nn, pp, n and p, which softens the 5 S 2 N interaction. The model uses the free dibaryon wave function to average the potential between dibaryon and baryon, where ψ d is the dibaryon wave function which is consisted of particle 2 and 3, U F (R F ) is average potential and the 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 (R F , r d ) = ψ d (r) ψ mole (R F ) and total binding energy is E = E d + E mole , where ψ mole (R F ) is the molecular state wave function calculated with the average potential U F (R F ) and E mole is the binding energy of molec-ular state. The binding energies of nn, pp, n and p are calculated by the folding model and their errors are estimated from the fitting error of N − and − 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.

Wigner function
The Wigner function introduced in Eq. (1) is written as [27,31,32], where r = (r 1 , r 2 ), q = (q 1 , q 2 ) are the relative coordinate and momentum, and ψ( 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, A Laguerre-Gauss quadrature is applied for the integrals of hyperradius R and θ 1 , θ 2 is integrated by a Legendre-Gauss quadrature [53]. The coordinate is defined in a sixdimensional spherical coordinate as R = (R, θ 1 , θ 2 , θ 3 , θ 4 , θ 5 ), which can be transformed into the six-dimensional Cartesian coordinate: R sin θ 1 sin θ 2 sin θ 3 cos θ 4 , R sin θ 1 sin θ 2 sin θ 3 sin θ 4 cos θ 5 , R sin θ 1 sin θ 2 sin θ 3 sin θ 4 sin θ 5 ), (22) and 0 ≤ θ i ≤ π (i = 1, 2, 3, 4), 0 ≤ θ 5 ≤ 2π , the volume element is The r in Eq. (21) is set at (r, 0, 0, 0, 0, 0), the q is set at (q, θ, 0, 0, 0, 0). By integrating out the angle, the probability to find the pn ground bound state can be obtained at six-dimensional hyperspherical radius r and at six-dimensional hyperspherical momentum q [52], 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) ∼ (2 fm, 200 MeV). And the normalization of the probability, If the wave function relates to not only ρ but also α, 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. ψ(ρ, α) = ρα| ψ can be simplified into n 1 ,n 2 r 1 r 2 | n 1 n 2 n 1 n 2 |ψ . r i |n i is a 3-dimension radial orthogonal basis which is the same as (16) but the last term is r −1/2 i with the same variation parameter λ for different n i . Here n i ranges from 2 to 26 with λ = 10000. By this way, Wigner transformation can be rewritten as:  P(r, q) of pn , which represents the probability to find the pn ground bound state at binding energy 21.7 MeV at six-dimensional hyperspherical radius r and at six-dimensional hyperspherical momentum q A complex Wigner transformation is simplified by a series of three-dimension Wigner transform. For the folding model, where ρ W di is the Wigner density function for dibaryon and ρ W di−b 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 [19] 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%.

Result and discussion
In blast-wave model, the parameters (τ 0 , τ , ρ 0 ⊥ , R 0 and T kin ) are fitted with experimental transverse momentum spectra of proton and by Eq. (4) and adjusted with the results of triton for different collisions, as shown in Fig. 3. Table 5 listed the parameters used in this work.
The transverse momentum spectra of pn is calculated by using the blast-wave model coupled with coalescence model (BLWC) as Eq. (1) and shown in Fig. 3a for Au+Au collisions at √ s N N = 200 GeV and Fig. 3b for Pb + Pb collisions at √ s N N = 2.76 TeV. The results of nn, pp, n and p with the relative wave function Table 5 The blast-wave model parameters for proton ( ) in Au + Au collisions at √ s N N = 200 GeV [66], which is fitted from the RHIC data [67,68] as well as in Pb+Pb collisions at √ s N N = 2.76 TeV [19] fitted from the ALICE data [69][70][71] T (MeV) ρ 0 R 0 (fm) τ 0 (fm/c) τ (fm/c)       [67,72,73] and ALICE [69][70][71]. The lines are calculated by our previous work [19] from the folding model are shown in Fig. 4. The p T spectra of p and from our previous work [19] as well as this work are also presented in Fig. 4. The p T spectra of n is not shown here because it is almost as same as p .
To further investigate the productions of -dibaryons and hypernuclei, the p T integrated yields d N/dy at midrapidity are given in Tables 6 and 7. The predicted results show N ∼ ×10 −4 [19], ∼ ×10 −7 , N N ∼ ×10 −7 and N ∼ ×10 −9 . The uncertainties of the integrated yields are directly from the Wigner functions, whose relative errors are about 20%. So the relative errors of yields are considered as 20%. Though the uncertainties from the blast-wave parameters are also important, which have been discussed by other model work [28], it will not be discussed in this paper. And  With the growing of constituents number A such as → N → N N and N → N → N , the production rates appear to follow the exponential function exp(−b A), here b is the so-called reduction factor [75][76][77], as shown in Fig. 5 for Pb + Pb collisions at 2.76 TeV. This A-dependent trend is similar to that for light nuclei of p → d → t ( 3 H ) in Fig. 5. However, it can be seen that n -n ( pp) slightly deviate from the trend in → N → N N . Keep in mind that the treatment of interaction is slightly different between pn and dibaryon-baryon via the folding method, which results in the slight deviation. In general, we have two classes for these production chains. One is for N → d → t ( 3 H ), → N → N N and → N (solid lines), they are almost parallel with the increase of N constituent number. Another is for N → N → N , → and d → N N chains (dash lines), they are almost parallel with the increase of number. Obviously much larger reduction factor b for the second class than the first class, indicating that much less yield for adding one more than one more nucleon. The different reduction factor b results from the different interactions between N − and − as well as the difference of productions of N and . Inspired by this, the production of hypernuclei N n H m (N for nucleons and H for one kind of hyperons) can be estimated by the intersection of N i H m and N n H j chains (i( j) is smaller than n(m)). Even if there is one point on the chain of N n H j , the reduction factor b of this chain is similar to the chain of H j or other chains whose b is known in the same class with N n H j . From Fig. 5, the prediction of the N N production is about 10 −10 . It implies that the production of hypernuclei is sensitive to the interaction among the constituents in the coalescence framework and then the systematic measurement of hypernuclei can shed light on the production mechanism and the baryon interaction.  N N can weak decay through an decay, which decays into -hypernuclei (A = 3) or -hypernuclei (A = 3), note that -hypernuclei (A = 3) can not be formed according to the HAL-QCD's results but might be formed under the ESC08c potential [78]. N N can also strong decay into N or N which is based on the interaction N − and N − reported by the HAL-QCD [79]. As for N , it can decay into N or N and mesons from the weak decay of . It can also decay into or by strong interaction. All here mentioned three-baryon group, such as N and N , may not be bounded. From Fig. 4, it is hard to figure out the difference between -N and N -. Although the p T spectra of -N and N -are almost the same, the strong decay channels are different in the folding model. For N -, it would decay into a or through the N − or N − channel and an , while the -N can hardly decay into or , since the N and are not bound directly in this folding model.

Summary
The three-body bound state problem can be solved through a variation method coupled with an eigenvalue problem. For weakly-bounded triple-particle system, the folding model is applied. The N -and -potentials used in this work are fitted from the lattice QCD's simulation near the physical point, which was reported by HAL-QCD collaboration. In coalescence model, the phase-space information of nucleons and are generated by the blast-wave model and the particles are coalesced into N N and N by using the Wigner density function from the simplified three-body wave function. The production of N N is about 10 −7 and N is about 10 −9 . There are also A-dependent trends similar with that for p → d → t ( 3 H ). The production rates follow the exponential function exp(−b A). With adding different baryon, the reduction factor b is different. Due to different factor b, two classes of hypernuclei chains will intersect at certain points where the production rate of new hypernucleus could be estimated. And the decay modes of N N and N are briefly discussed in order to search for such exotic triplebaryons (hypernuclei) in future experiments, which could provide a method to understand the Y N and Y Y interactions for multi-strangeness hadrons. The systematic measurements of hypernuclei can definitely shed light on the production mechanism and baryon interactions.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study and no experimental data has been listed.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permit-ted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .