Nucleon clustering at kinetic freezeout of heavy-ion collisions via path-integral Monte Carlo
aa r X i v : . [ nu c l - t h ] S e p Nucleon clustering at kinetic freezeout of heavy-ioncollisions via path-integral Monte Carlo
Dallas DeMartini and Edward Shuryak
Department of Physics and Astronomy, Stony Brook University, Stony Brook NY 11794-3800, USA
Clustering of the four-nucleon system at kinetic freezeout conditions is studied using path-integralMonte Carlo. This method seeks to improve upon recent calculations which relied on approximatesemiclassical methods or few-body quantum mechanics. Finite temperature and finite density effectsare included using data from recent heavy-ion experiments. The effects of a possible modification ofthe inter-nucleon interaction due to the presence of the QCD critical point on the clustering is studiedas well. From this, estimates are given for the decay probabilities of the 4N system into various lightnuclei decay channels and the strength of spatial correlations is characterized. Additionally, a simplemodel is presented to describe the impact of this clustering on proton multiplicity distributions.
I. INTRODUCTIONA. Nucleon clustering
One of the major goals of current heavy-ion collisionexperiments has been the search for signals of the theo-rized QCD critical point. The most prominent of these isthe Beam Energy Scan (BES) at the Relativistic Heavy-Ion Collider (RHIC), which looks at Au + Au collisions atenergies √ s NN = 5 - 39 GeV. Several observables signal-ing the existence of the critical point have been put for-ward. Two experimental observations of particular rele-vance to this work are the increase in ratios of light nucleiat √ s NN = 20 - 30 GeV by the STAR collaboration [1]and the modification of proton multiplicity distributions(kurtosis) [2].Recent work has suggested that modifications to theproduction of light nuclei may be a signal of correlationsand fluctuations induced by critical phenomena [3]. Sim-ple coalescence models [4] suggest that the productionratios for light nuclei such as N t N p /N d should be inde-pendent of beam energy in the absence of strong spatialcorrelations between nucleons. It is precisely such spatialcorrelations that are the focus of this work.In Ref. [5] classical molecular dynamics and theLangevin equation had been used. Strong dependence ofclustering on the inter-nucleon potential was observed. Itelucidates a special role of the 4-nucleon system, as thefirst cluster which can withstand the relevant tempera-tures of kinetic freezeout T kin ∼
110 MeV.In Ref. [6] two approximations had been developed,to calculate thermal/quantum density matrix of severalnucleons: (i) semiclassical method based on “flucton” so-lution; and (ii) the solution of hyperradial Schrodingerequation. The latter leads to explicitly finding the secondbound state, even in the lowest angular momentum case,and, more generally, recognized the existence of multi-ple excited states of 4N system ( He nuclei or alpha-particle). Certain experimental signatures of preclusterformation, leading to those excited states, has been pro-posed, such as possible two-body decay channels d + d , t + p , or He + n . However, both approximate methodstreat the 4-nucleon thermal states focusing on the low- est hyper-harmonics (Schrodinger equation depending onthe hyper-radial coordinate ρ only), while among the ex-perimentally known resonances one finds clear dominanceof states with orbital momentum L = 1 and L = 2.Therefore, these calculations perhaps include only a partof occurring multibody preclustering effects.To clarify the issue completely, it is desirable not to useany approximations at all, treating the few-body quan-tum/thermal dynamics from first principles. This is thepurpose of this paper, in which we make further studiesof the 4-nucleon systems at relevant temperatures, nowusing the most straightforward numerical method, thepath-integral Monte Carlo (PIMC).We focus on N = 4 system – alpha-particle – for threereasons. First of all, it is the largest N correspondingto nucleon spin-isospin states: therefore all four nucle-ons can be in different states p ↑ p ↓ , n ↑ n ↓ and thereforedistinguishable, which greatly simplifies PIMC applica-tion. That is to say that the particles may be treated asquantum boltzmannons, particles which obey quantummechanics and classical statistics. It is for this reasonthat PIMC was first used for this system by one of uslong ago [7].The second reason: it is the first precluster bounddeeply enough, to withstand such temperatures. Let usremind the reader that the binding energy of He is E B =28.3 MeV, while the binding energy of d is only E B =2.2 MeV.The third reason: unlike d, t systems, it has sufficientlylarge number of excited states. Their theoretical andeventually experimental studies can be used as a veryuseful tool, highly sensitive to the inter-nucleon poten-tial under thermal conditions. Additionally, rather thanassuming the cluster decays statistically into the groundstate and the excited state ’precluster,’ the PIMC ap-proach allows us to separate the ground state from thedynamically-generated ’precluster’.While the binding energy of the 4N cluster is toolow to remain bound at T kin , it is still strong enoughto induce significant clustering due to the maximumdepth of the potential V min being comparable to thetemperature. Much stronger clustering can be induced( | V min | > T kin ) near the critical mode where the effec-tive mass of the (attractive) σ -exchange is reduced. Thisleads to large Boltzmann factor that plays an even largerrole when more nucleons are involved; the 4N system hassix nucleon-nucleon interactions. B. Freezeout in heavy-ion collisions
Ultra-relativistic heavy-ion collisions produce hot, de-confined QCD matter, the Quark-Gluon Plasma (QGP).The QGP formed in these collisions quickly expands andcools. As it cools, the quarks and gluons become confinedand hadronize before they reach detectors. In additionto the confinement phase transition (which is a crossoverat small µ B and first-order at larger µ B &
400 MeV),the system undergoes two other transitions: the chemi-cal and kinetic freezeouts.The chemical freezeout occurs when the system hascooled enough for inelastic collisions between hadrons tocease. At this point, the total yields all produced par-ticles are fixed. The temperature T ch and baryon chem-ical potential µ Bch of the are found by use of a statisti-cal hadronization model [8], which fits the parameters toparticle yield ratios.For some time after this freezeout, elastic collisionscontinue. Once the system has expanded so much thatelastic collisions cease, it undergoes kinetic freezeout. Be-cause this is the end of interactions for these particles, thetransverse momentum p T is fixed beyond this time. Thetemperature of the kinetic freezeout T kin can be calcu-lated from the hydrodynamically-motivated ”blast wave”model [9, 10].For collision energies above 7 . T ch is onlyweakly dependent on collision energy and remains nearthe deconfinement temperature T c . Similarly, T kin re-mains on the order of 100 MeV for the same range ofcollision energies. For lower collision energies (higherchemical potential), both of these temperatures decreaserapidly and T ch ≃ T kin . See Table V for the specificvalues in the relevant collisions. At the largest collisionenergies, T kin begins to decrease due to the very largemultiplicity of the fireball at high energies. Indeed, ininfinite nuclear matter, there is no kinetic freezeout aselastic collisions never cease.This paper is laid out as follows: In Section II, we dis-cuss the basic details of the Monte Carlo simulation andsome preliminary calculations. Section III lays out theidea of the nucleon ’precluster’ and results of simulationsof the kinetic freezeout conditions. Section IV looks atthe effects of modified inter-nucleon interactions and Sec-tion V describes the observables that can be estimatedfrom this work. II. TECHNICAL DETAILSA. The setting: nucleon interactions and the pathintegral
The fundamental quantity of interest in quantum sta-tistical mechanics is the path integral: h ~x ( t f ) | e − ˆ H ( t f − t i ) | ~x ( t i ) i = Z Dx ( t ) e − S [ ~x ] , (1)where R Dx ( t ) is the integral over all possible particlepaths x ( t ). For finite temperature, this is the integralover all paths with a period β such that ~x ( t + β ) = ~x ( t ).The classical Euclidean action S [ ~x ( t )] of a path ~x ( t ) is S [ ~x ( t )] = Z β dt (cid:18) m ˙ ~x ( t ) + V ( ~x ( t )) (cid:19) . (2)(Note that the sign of the potential V is inverted in Eu-clidean notation.) In order to perform a numerical eval-uation of this path integral, Euclidean time must be dis-cretized into N t time slices with spacing a = β/N t . Thediscretized action then becomes S [ ~x ] = a N t X j =1 (cid:16) m a ( ~x j +1 − ~x j ) + V ( ~x j ) (cid:17) . (3)Here ~x j is the position of the particle on the j-th timeslice and the periodicity condition is imposed by setting ~x N t +1 = ~x . This discretized action is exact in the limit a → N t → ∞ .For the 4-particle, 3-dimensional system we are inter-ested in, the path integral is a 3 · · N t -dimensional in-tegral. This integral is evaluated by way of the standardMetropolis algorithm, where the change in the statisticalweight e − S , calculated as in Eq. 3, serves as the prob-ability for accepting or rejecting updates to particle po-sitions. From the configurations generated this way, wemay calculate any desired quantities: average energies,correlation functions, etc. While the calculation of mostquantities is straightforward, it can be seen in Eq. 3 thatthe kinetic energy term diverges in the limit a →
0. Inorder to avoid this when calculating the average kineticenergy, it is more useful to use the virial estimator [11] h KE i = 32 N T − h N X i =1 ( ~x ij − ~ ¯ x i ) · ~F ij i . (4)Here x ij is the position of particle i on the j-th time sliceand F ij is the classical force on particle i at the j-th timeslice. ~ ¯ x i is the time-averaged center of mass of particlei. This form avoids the a − divergence as a → β , which disappears in the high temperaturelimit as | ~x ij − ~ ¯ x i | → He nucleus or four nucleons suchthat n N V = 4. Simply placing the system in a hard-walled box of volume V would introduce unphysical ef-fects. Rather, we include a set of images on each face ofthe box. These are six sets of nucleons shifted in x, y, orz by a distance L = V / . This gives us a total of sevenboxes with the main simulation box in the center. Whenthe center of mass of a nucleon travels outside the sidesof the box, it and all of its images are shifted by L to theother side of the box. See Appendix C for a discussionof the periodic box and images.The main inter-nucleon interaction that is consideredis the Serot-Walecka potential [12] V SW ( r ) = − α σ e − m σ r r + α ω e − m ω r r , (5)with α σ = 6 . m σ = 500 MeV, α ω = 15 .
17, and m ω = 782 MeV. While this potential describes proper-ties of infinite nuclear matter well, it does not possess abound state in the 4 N system. We therefore modify therepulsive term, using α ω = 11 .
02 which approximatelyreproduces the desired binding energy. In Section IV, weconsider the effects of modifying the potential with a re-duced value of m σ . Note that this simplified potential isiso-scalar, so that it is the same for nn, np, pp pairs. Ad-ditionally we consider the repulsive Coulomb interactionbetween the two protons V C ( r ) = 7 . · − /r . (Notethat these energies are in units of f m − and that 1 fm − ≃
197 MeV.)If there exists the QCD critical point, a near-masslesscritical mode should be strongly fluctuating in its vicin-ity. These fluctuations are nonlinear, inducing certainthree- and four-nucleon forces. Calculation of those willbe a subject of our subsequent publication [13].
B. The ground states of two and four nucleons
The first task is to tune PIMC parameters and to testwhether the ground state binding and the wave func-tions are correctly reproduced. Naturally, one wouldlike to make sure that the PIMC code reproduces resultswhich are well-known. We tested the code by producingthe ground state probability distribution of the Waleckadeuteron - two nucleons interacting via the Serot-Waleckapotential. Being a two-body system with only a radial in-teraction, the Schrdinger equation can be solved withoutapproximation.Calculations for all ground state distributions of the4N system were done at T = 0 . He (see Table I). This leads to aMatsubara time of 394 fm, discretized into N t = 3000time steps. FIG. 1: Normalized probability distribution for the Waleckadeuteron. The curve is the numerical solution of theSchrdinger equation.
The system of four interacting nucleons exists in a 9-dimensional space. In order to study the dynamics of thesystem, it is useful to use the hyperdistance, defined asthe quadrature sum over the distances between all (six)pairs of nucleons i and j ρ = 14 X i Calculated paths are used to show distributions in hy-percoordinate ρ , at different temperatures. We used the FIG. 2: Normalized probability distribution of the groundstate | χ ( ρ ) | of He over hyperdistance ρ . unmodified Serot-Walecka model potential, as well as themodified ones.The PIMC simulations performed for the kineticfreezeout conditions use the temperatures listed in Ta-ble V. For the 2.4 GeV collisions, we see from Ref. [14]that T ch = 65 ± T kin = 71 ± T kin ≤ T ch , we use the temperature and nucleon den-sity calculated from the condition of chemical freezeoutfor this energy. All of these runs have Euclidean timediscretized into N t = 40 time slices.Simulations were performed for the six collision ener-gies in Table V from √ s = 2 . ρ clustering region.For the temperatures and nucleon densities corre-sponding to the kinetic freezeout conditions, we find thedistribution in ρ . This distribution can be split into threedistinct, dynamical parts: the ground state, the ’preclus-ter’ consisting of all excited states, and the long-distancetail corresponding to propagating positive-energy states.The ’precluster’ can be identified by subtracting out theother two contributions from the total probability distri-bution. This ’precluster’ and its relative size are the mainobject of study in this work; it is the distribution whichcan be decomposed into the excited states of He thatwill then decay into the various two-body states. Herewe will discuss step-by-step how the distribution is de-composed into the three essential parts for √ s NN = 19 . FIG. 3: (Color online) Probability distributions | ψ ( ρ ) | of theground state (blue • ) and kinetic freezeout distribution (red (cid:4) ) in hyperspace ρ for √ s NN = 19 . tive normalization is arbitrary and they must be nor-malized appropriately. This is done by normalizing thekinetic freezeout distribution so that the ground statedistribution fits underneath it. Fig. 3 shows the resultsof this.Now the ground state distribution can be easily sub-tracted out. Next comes the less-trivial subtraction of thethermal tail. One should notice that the kinetic freezeoutdistribution does not become flat at long distances, as oneshould expect for such a distribution. This is an effectof the finite size of the box. In infinite space, one wouldexpect | χ ( ρ ) | to go as ρ at large hyperdistance (giventhat this is a 9-dimensional system with 8 hyperangularcoordinates), meaning that | ψ ( ρ ) | should be constant.One finds that this distribution does rise almost as ρ forintermediate values of ρ , but as ρ becomes comparableto the length of the box L , the phase space of possibleconfigurations decreases and eventually goes to 0 at themaximum value of ρ that fits in the box.This all means that the thermal tail is not just a con-stant at large ρ , but rather a ρ -dependent distributionto be subtracted out. This distribution can be found bysimply generating random configurations of four nucle-ons in a box of the appropriate size. These distributionstypically have poor statistics in the small- ρ regime due tothe probability going as ρ for small ρ . Because of this,the distributions are fitted to a function of the form f ( ρ ) = Ae ( ρ − B ) /C + 1 (9)to correct the smallest- ρ data points as seen in Fig. 4.The form of this fit is not physically motivated, rather itis a simple form that describes the distributions reason-ably well. Fortunately, the determination of the preclus-ter is insensitive to the small- ρ behavior of the randomdistribution. FIG. 4: Probability distributions | ψ rand ( ρ ) | of fourrandomly-placed nucleons in hyperspace ρ at a density equiv-alent to kinetic freezeout conditions at √ s NN = 19 . By dividing the kinetic freezeout distribution by therandom distribution, we effectively eliminate the mod-ification of the phase space due to the geometry of thesystem. This reveals a distribution much more similar anthe infinite-space distribution (see e.g. Fig. 5 in Ref. [6]).An important point to note here is that while dividing bythe random distribution removes the effects of the geom-etry, it does not remove the effects of the images. Thisexplains why the distribution does not completely flattenand, in fact, drops off quickly at large ρ : At the largestvalues of ρ the nucleons are near the sides of the box andmay be interacting very strongly with the images.Similarly to the ground state, the random distributionmust be appropriately normalized relative to the kineticfreezeout distribution. Unlike the ground state distri-bution, the random distribution should not fit entirelyunder the kinetic freezeout distribution. The questionis, to which region of the distribution should the twobe equal. The answer is that it should be fit to the re-gion of the distribution where the inter-nucleon potential h V NN i ≃ ρ means that they make up a significantly largercontribution to the density matrix when the ρ factor isincluded in the integration. The ground states typicallymake up just a few percent of the cluster, which is tobe expected given the ∼ 50 excited states that make upthe precluster. Similarly, the cluster makes up a smallfraction of the overall thermal distribution. Despite sig-nificant clustering being possible at these temperatures,the 4N system is still fragile, being dominated by con-figurations where at least one nucleon is well separatedfrom the rest.For the five beam energies with nearly the same tem-perature, one can see that the relative sizes of the groundstate and precluster distributions remain roughly thesame. The main trend that is seen in increasing the beamenergy is the lengthening of the thermal tail due to thedecreased density and larger box size. In the limit thatthe length of the box is much larger than the values of ρ we are interested in, the thermal tail should flattencompletely. FIG. 5: Ratio of kinetic freezeout distribution to random dis-tribution. Dashed line is the average value of the flat region3 fm < ρ < The set of known excited four-nucleon states is repro-duced in the Table I. Note that the energy here is theexcitation energy, counted from the ground state. Sincebinding is -28.3 MeV, the highest state (at the bottomof the table) has absolute energy of about 1 M eV abovezero. Note however, that its width is an order of mag-nitude larger than that, and also larger than ∼ FIG. 6: (Color online) Normalized probability distribution in hyperspace ρ of the ground state (blue • ), kinetic freezeoutdistribution (red (cid:4) ) and the precluster (green N ) for the six beam energies considered B. Angular distributions One may also gain information about the density ma-trix by looking at the distribution in angles formed bythree nucleons cos ( α ) ≡ ( ~r i − ~r j ) · ( ~r i − ~r k ) | ~r i − ~r j || ~r i − ~r k | (10)where i = j = k are three different nucleons. Note thatclassical minimum of the potential energy corresponds totetrahedral shape of 4-nucleon cluster, in which all sidesare equilateral triangles and therefore all angles are suchthat cos ( α ) = 1 / 2. This tetrahedral minimum-energyconfiguration is modified only slightly by the relatively-weak Coulomb interaction. There was considerable dis-cussion of α -particle-like clustering in the ground statesof light nuclei, e.g. C , O .In Fig. 7 we see, as before, that the angular distribu-tion of the ground state is both broad and not centeredat cos ( α ) = 1 / 2. The broad distribution is the resultof quantum/thermal fluctuations around the minimum-energy configuration stemming from the kinetic term inthe action. The shift in the peak to an angle near 36 ◦ isdue partially to the Coulomb interaction of the protonsand also due to the fact there is simply more phase spaceto produce a triangle with small angles than large ones,as can be seen by the random distribution.The kinetic freezeout distribution is nearly identical FIG. 7: (Color online) Normalized probability distribution ofinternal angles P (cos α ) for the ground state (blue • ), kineticfreezeout distribution (red (cid:4) ), and the random distribution(green N ) for conditions at √ s NN = 19 . to that of randomly-placed nucleons. This makes senseas the system spends most of its time at large valuesof ρ outside of the clustering region where the average TABLE I: Low-lying resonances of the He system, from BNLproperties of nuclides J P are the total angular momentumand parity, Γ is the decay width. The last column is thedecay channel branching ratios, in percent. p, n, d correspondto the emission of proton, neutron, or deuterons respectively. E (MeV) J P Γ (MeV) decay modes, in %20.21 0 + p = 10021.01 0 − n = 24, p = 7621.84 2 − n = 37, p = 6323.33 2 − n = 47, p = 5323.64 1 − n = 45, p = 5524.25 1 − n = 47, p = 50, d = 325.28 0 − n = 48, p = 5225.95 1 − n = 48, p = 5227.42 2 + n = 3, p = 3, d = 9428.31 1 + n = 47, p = 48, d = 528.37 1 − n = 2, p = 2, d = 9628.39 2 − n = 0.2, p = 0.2, d = 99.628.64 0 − d = 10028.67 2 + d = 10029.89 2 + n = 0.4, p = 0.4, d = 99.2 inter-nucleon potential h V NN i ≃ 0. The cluster at kineticfreezeout, which is a sum of states of various angular mo-menta and the thermal tail, shows no angular correlationbetween the nucleons. The distributions in Fig. 7 arequalitatively identical for all beam energies tested.The wide distribution in internal angles, even in theground state, shows the importance of an approximation-free method for few-body quantum systems such asPIMC. The previously-used methods, reducing the sys-tem to a one-dimensional quantum system or the semi-classical ”flucton” method, assume a symmetry in inter-nal angles which is not found in the majority of configu-rations. IV. MODIFICATION OF THEINTER-NUCLEON POTENTIAL In heavy-ion collisions the medium is expected to mod-ify the parameters of the inter-nucleon potential. This isespecially true near the critical point, where long-rangecorrelations should arise. The most dramatic effect sug-gested of this form is the reduction of the σ mass m σ near T c [25] m σ ∼ (cid:18) | T − T c | T c (cid:19) ν . (11)Such a reduction in the σ mass, which drives the attrac-tive portion of the interaction, will greatly increase thestrength and range of the attraction. A similar decreasein the ω mass is not expected, meaning the repulsion FIG. 8: Binding energy E B of the 4N system as a function ofthe σ mass m σ should not increase to compensate. Clearly such a mod-ification will greatly modify the clustering dynamics ofthe 4N system of interest in this work. The most obvi-ous question to ask is then: how much modification of m σ should be expected near kinetic freezeout? The an-swer is that it should be only a small modification for twomain reasons. The first is that the reduction of m σ to 0 isexpected at T c which is some ∼ 40 MeV above T kin for awide range of beam energies. The second is that even atthe critical point, the finite size and finite lifetime of theQGP system prevent the correlation length from gettingtoo large, and thus, m σ from getting too small.The strengthening of the attractive force is seen clearlyin Fig. 8. Here the binding energy of the 4N system isshown for the σ mass reduced down to 300 MeV. Re-ducing the mass down to this value increases the bindingenergy from its physical value of 28 . m σ studied, 300 MeV, thebinding energy is greater than a GeV. Clearly such adeeply-bound state is unrealistic. This is a limit of thesimplicity of our model. In addition to the binary inter-action included here, one should also include three- andfour-body interactions as well. These interactions shouldalso be sensitive to m σ and the three-body interactionin particular should be strongly repulsive, which shouldprevent the system from becoming so deeply bound atsmall m σ . These interactions, while not included here,will be the subject of our future studies.The effect of the modified interaction on clustering isseen directly in Fig. 9. Reducing the standard σ massby 50 MeV causes the peak correlation to jump from ∼ ∼ 20. Similar values are seen for all but the lowestbeam energy, which all have approximately equal T kin .For √ . T kin is reduced, the peak cor- FIG. 9: (Color online) Probability distribution relative to ran-dom distribution P ( ρ ) at the kinetic freezeout conditions for √ s = 19 . m σ = 500 MeV (blue • ) and m σ = 450MeV (red (cid:4) ) relations observed are ∼ ∼ m σ ≤ 400 MeV for kinetic freezeout at √ s = 2 . m σ = 400 MeV. V. ESTIMATES OF THE OBSERVABLESA. Decay channels of the 4N system The main experimental observable of the clustering ofthe 4N system is the yield of He and its decay products.One of the goals of this work is to determine dynamicallyrather than statistically the fraction of these correlatedclusters which decay into the ground state. Withoutknowledge of the wave functions of the various excitedstates, the precluster is assumed to populate the vari-ous states in Table I statistically. These excited statesthen decay into light nuclei with probabilities that areknown experimentally. Note that only one decay prod-uct is known for each decay mode. We assume all decayslisted in Table I are two-body decays. In principle, thesystem could decay into three nuclei ( p + n + d ) or fournucleons ( p + p + n + n ).With this, the probability of the correlated cluster todecay into either one of the three possible two-body decaychannels or to populate the He ground state is deter-mined for all of the beam energies studied.At all but the lowest energy studied, the ground statemakes up only a few percent of the total cluster. Only TABLE II: Decay probability (in percent) of the 4N clusterinto He or two-body states for the different collision energies.Three- and four-body decays are not considered here. √ s (GeV) He p + t n + He d + d at √ s = 2 . He comparable to the individual two-body decay probabili-ties. This is due to the significant reduction in T kin atthis energy increasing the relative statistical weight ofthe ground state. Finite-density effects are much smallerthan the effect of changing temperature. For the fiveenergies with roughly equal T kin increasing the densitymodestly increases the ground state probability. Thismay be attributable to the size of the box being com-parable to or smaller than the size of some of the ex-cited states, removing them from the spectrum of allowedstates, while the smaller ground state does not feel theeffects of the boundary.If one considers the high temperature and low densitylimit, the ground state should have a probability of just2% given that there are 49 (known) excited states of He .For the highest energies studied, in which the density islow and T kin > ∆ E , we find the ground state probabilityto be just above 2%. This suggests that our system,which does not include fermion statistics and uses thesimple Serot-Walecka potential, has a reasonable numberof excited states compared to the known He states.Decays of the 4N system should contribute only a smallamount to the total yield of these particles. A completepicture of light nuclei production should include statisti-cal hadronization as well as feed-down from excited nu-cleon states (such as N ∗ or ∆), which should greatlyoutnumber 4N states. B. Virial expansion and clustering The potential part of the partition function (of singlespecies system) of N particles can be re-written in theform Z pot = 1+ 1 V N Z d x ... Z d x N (cid:2) e (cid:0) − P i>j V ( ~x i − ~x j ) /T (cid:1) − (cid:3) (12)by adding and subtracting 1. Since we focus on clustersof 4 particles, coordinates of all others can be integratedout, as well as the coordinate of its center of mass. Onemay study the 2 N and 3 N systems to determine the sizeof the effect proportional to n and n , respectively. How-ever, given the small binding of d and t , one should expectsuch contributions to be small. What is left is Z pot = 1+ N ( N − N − N − V (9) cor V ) ≈ n V (9) cor N V (9) cor is V (9) cor = 32105 π Z dρρ ( P ( ρ ) − , (14)Here P ( ρ ) is the probability distribution relative to anideal gas in 9-dimensional radius ρ shown in Fig. 5, andthe factor in front is the solid angle in 9 dimensions. Weneglect repulsion and integrate over the region in whichthe integrand is positive. The addition to free energyis then ∆( − T log ( Z )) = − T n V (9) cor N , same as to thegrand partition sum. Differentiating it over µ , presentin each N , one finds the addition to the particle number∆ N/N = n V (9) cor / 3! (a factor of 4 cancels out).As a check on the numerical factor in the thermody-namic expression, we can compare it to a more ’direct’method of computing the ratio of clusters. We can de-fine the total volume of the entire distribution V (9) tot anal-ogously and then compute to ratio R of clusters to un-clustered configuration R = V (9) cor / ( V (9) tot − V (9) cor ).The expressions are a bit modified in the case of sev-eral particle species. In the problem at hand, nucleonshave spin 1/2 and isospin 1/2, so the number of distinctspecies is 4. If n s is density per species , the total densityof symmetric matter is simply n B = 4 n s . In the caseof particular clusters we actually simulate, made of fourdistinct species, there is no need for symmetrization andthere is no 4! = 24 in denominator. However, the densityin front is in this case n s , not total n B , and the numericalsuppression factor is 1 / = 1 / 64. Therefore, the clus-tering contribution to h N i is small, at the percent levelor less. TABLE III: Correlation volume V (9) cor of the 4N system at allbeam energies for m σ = 500 , , 400 MeV. √ s (GeV) V (9) cor (fm ) m σ (MeV) 500 450 4002.4 8 . · . · . · . · . · . · . · . · . · . · . · . · 27. 9 . · . · . · 39. 9 . · . · . · Furthermore, the n -th derivative of log ( Z ) over chem-ical potential, called K n , can also be calculated. Onefinds that K , with extra three derivatives compared to K = h N i , does not have this 1 / numerical factor, andso, in the same approximation as above, K h N i − n B V (9) cor (15)The values of the r.h.s. are shown in Fig. 10. As one cansee, the predicted cluster contribution to the 4th momentare no longer small, for the two left points, correspondingto HADES and the lowest BES-I energy √ s = 7 . GeV .This is precisely what is observed. FIG. 10: The 4th cumulant deviation (Eq. (15)) versus √ s ,using the 9-dimensional correlated volume V (9) cor determinedfrom the PIMC simulations. Another perspective at the observed cumulants K n canbe obtained using f actorial cumulants C n . In particular K − h N i = 7 C + 6 C + C . (16)According to [15], BES-I data show small values for C , C and K is in fact dominated by the factorial cu-mulant C . It correlates well with our general finding,that clustering starts from the 4N systems, and with thedependence of K , C on collision energy.Characterizing the strength of spatial correlations inthe 4N system is necessary for making predictions of theoverall magnitude of the feed-down contributions. Mod-els of light nuclei production, such as the previously-mentioned coalescence model [4] show explicit depen-dence on spatial correlations of nuclei. Our correlatedvalue is a measure of such correlations.From Tables II and IV one can then estimate thatabout 25% − 30% of the clusters decay into t . The re-sulting ratio of tritium to proton yields t/p is then about0 . 4% at √ s = 7 . V (9) cor for all six beam ener-gies and for three values of m σ . It must be noted that0these values should only be taken with a moderate levelof uncertainty. The correlated volume is sensitive to themethod of subtracting out the long-distance thermal tailfrom the distribution. Even a change of a few percentin the normalization factor (which comes from averagingover the flattest portion of the distribution P ( ρ )) can leadto a factor of 2 change in V (9) cor . Such a change results in aslight modification of the length of the largest- ρ portionof the cluster, which when multiplied by ρ drasticallychanged the integral. It is clear that the largest source ofuncertainty (and room for greatest improvement) comesfrom this subtraction of the thermal tail. Further stud-ies should not solely rely on hyperdistance ρ , but ratherdistributions in the full set of 9-dimensional hypercoor-dinates. TABLE IV: Ratio of clusters to unclustered configurationscalculated from the thermodynamic contribution expected inEq. (13) and the ratio R computed from direct integrationover distributions. √ s (GeV) n B V (9) cor R Knowing this, let us briefly discuss some qualitativefeatures of these values. For fixed values of m σ values atdifferent energies vary by up to a factor of ∼ 5, with thefour highest energies being roughly equal for all values of m σ . These four energies have nearly equal temperatureand the lowest densities, meaning they are affected leastby the finite system size. Decreasing m σ by 50 MeVincreases V (9) cor by ∼ − VI. NUCLEON MULTIPLICITYDISTRIBUTION While the total baryon number is conserved in colli-sions, the observed multiplicity distribution shows fluc-tuations which may be caused by many effects: (i) wan-dering inside and outside the detector acceptance; (ii)turning of the observable protons into unobservable neu-trons; (iii) turning into light nuclei species, and, last butnot least, (iv) the preclustering phenomenon. The firsttwo are well studied with various event generators. Anti-correlations between proton and light nuclei multiplicitiesin individual bins are still to be experimentally studied.The reason why preclustering phenomenon affects thenucleon multiplicity distribution can be explained as fol-lows. Suppose a certain number of preclusters of four or more nucleons are formed at freezeout. Instead of focus-ing on their two-body decays to light nuclei as we didbefore, let us think about feed down into four individualnucleons.On one hand, the number of preclusters is relativelysmall, and their branching ratio to such modes is alsonot large, so one may expect that such feed down is com-pletely negligible. This is indeed the case near the max-imum of the multiplicity distribution N p ≃ h N p i . Yetthe effect should be large and observable at the tails ofthe multiplicity distribution. The main multiplicity dis-tribution is, to zeroth approximation, just a Poisson dis-tribution, with the mean defined by detector acceptance.Its tails, at small and large multiplicity, are the proba-bility that a large number of protons happen to cross theacceptance boundary and be either outside or inside it.For example, according to STAR data [17, 18] for central Au − Au collisions at √ s = 7 . h N p i ≃ 40 and the tails we discuss are at N p ∼ 10 and N p ∼ the mean number of clusters issmaller .Indeed such deviations from Poisson are observed. Inorder to quantify those, one now uses a sensitive statisti-cal tool, calculating the factorial cumulants of the distri-butions. Recall that for the Poisson distribution C i = 0for i > 1. For the STAR data in question, they are in-stead (the maximum values taken from all energies andcentralities) C ≈ − , C ≈ − , C ≈ . While the exact reason for these large values of C and C remain unknown, we present here an ad-hoc model,which quantifies the effect and qualitatively reproducesthe cumulants. (Its discussion was inspired by anothermodel proposed in Ref. [15].) It has the distribution P ( N p ) = 0 . P P ( N p , 40) + 0 . P P ( N c , . 1) (17)where P P ( N, h N p i ) is the Poisson distribution and N c isthe number of clusters of four protons. This particularexample gives C = − . C = 190.Of course, this is a schematic model. An assumptionthat a cluster decays into exactly 4 protons is unrealisticas some of them are neutrons. That is why we took thecluster probability normalization to be so low, 0.005.Another effect so far ignored is that clusters are treatedlike point-like objects, so they are either in or out of theobserved phase space for all nucleons. The clusters dohave certain coordinate and momentum spreads in thedensity matrix, so they may be outside of the volume partially . This effect can be quantified using ensemblesfrom our simulations. Dividing the volume in half along1 FIG. 11: (Color online) The two Poisson components of Eq.(17), proton component (blue • ) and 4 N cluster component(red (cid:4) ) any axis, we record the probability for observing N nu-cleons (0 − 4) in said half. The result is shown in Fig 12,showing that the probability distribution for interactingnucleons is wider than randomly-placed nucleons. FIG. 12: (Color online) Probability distribution of number ofnucleons N in half of simulation box for a random distribu-tion (blue • ) and for interacting nucleons at kinetic freezeoutconditions (red (cid:4) ) at √ s = 2 . σ mass. This, of course, modifies the nucleon multiplicity dis-tribution by increasing the probability that either all ornone of the nucleons in the 4N system are in the phasespace accessible to the detector. This effect should be greatly enhanced for the modified interactions resultingfrom the reduction of the σ mass. In fact, it is preciselysuch modification of the NN interaction which has gen-erated such interest in proton multiplicity distributions.The definitive feature of the QCD critical point is thelarge increase in correlation length ξ of the QGP sys-tem. Given that the correlation length may only increasemodestly at freezeout in heavy-ion collisions, much workhas been done on identifying observables which dependmore strongly on ξ . It has been proposed [19, 20] thatnon-Gaussian moments of the proton multiplicity distri-bution display such behavior. In particular, the kurtosis κ is expected to be sensitive to ξ . Much experimentaleffort has been exerted to study these distributions in-cluding at the programs discussed in this work, HADESand STAR [21, 22]. VII. SUMMARY A program of quantum/statistical mechanical studiesof few-nucleon clustering has been started by Refs. [5, 6].The theoretical methods used in those papers includeclassical molecular dynamics, semiclassical ”flucton” ap-proach at finite temperatures, and hypercoordinates in3( N − f our − nucleon systems, as the lightest one possessingmultiple states/resonances.The goal of this paper is to check calculations usingthis set of approximate methods by a direct first-principlecalculation based on path-integral Monte Carlo. We hadshown that the method works reliably, and, after fine tun-ing the potential, we also reproduce the binding of theground state of He . After that we performed calcula-tions at finite temperature and density. Our main resultsare shown in Fig. 6, in which we plotted the density ma-trix in terms of the hyper-radial coordinate ρ (Eq. 6). Tothe extent we can compare those results with those of ap-proximate methods mentioned, we see a rather consistentpicture.One can see that the precluster shape is different fromthat of the ground wave function squared, so that af-ter freezeout quantum mechanical decomposition of thepreclusters should produce not only He , but also a su-perposition of (near-zero-energy) bound and resonancestates. Those have close energies but different quantumnumbers (in particular, angular momenta), and thereforehave different decay widths and branching ratios. Fortu-nately, these states and their decays were experimentallystudied long ago, so in principle one can evaluate the feed-down from them into yield of light nuclei, such as d , t , He , and He . Recent summaries of experimental situa-tion can be found in [3, 24]. The most intriguing experi-mental observation is that the yield ratio N t N p /N d seemsto non-monotonously depend on collision energy, withapparent maximum at BES-1 mid-range √ s ∼ GeV .2The value of the ratio, especially from low energy HADESdata, are very different from ratio of statistical weights,indicating presence of strong feed-down into tritium.Another manifestation of clustering can potentially bestudies of moments of proton multiplicity distribution.Using the 9-dimensional correlation volume evaluatedfrom PIMC data, we calculated the fourth-order virialcoefficient of nucleon matter at kinetic freezeout. Whileits contribution to particle number is small, at the per-cent level, we found that its contribution to C / h N i be-comes of order one at collision energy √ s = 7 . N tail of themultiplicity distribution and propose a model for quali-tatively reproducing STAR data.If there exists a QCD critical point, one expects ad-ditional contributions to those from fluctuations of thecritical mode [25] when the freezeout occurs close to itslocation on the phase diagram. Nonlinear features ofsuch fluctuations are complicated, but fortunately wellknown e.g. from lattice studies of the Ising model. Asoutlook, we mention that our next publication [13] willdiscuss three- and four-nucleon forces induced by suchfluctuations. Acknowledgments This work is supported by the Office of Science, U.S.Department of Energy under Contract No. DE-FG-88ER40388. We also thank the Stony Brook Institutefor Advanced Computational Science for providing com-puter time on its cluster. Appendix A: Bulk properties of matter at freezeouts The RHIC beam energy scan, suggested to look forQCD critical point and whose first results were reportedin Ref. [26], took data at the energies listed in Table V.As one can see from the table, in the scan energy re-gion the temperature of chemical freezeout T ch is growing(errors ± M eV ) to a constant, while the baryon chem-ical potential µ Bch strongly decreases. The fireball vol-ume within this RHIC energy scan range approximatelydoubles. The temperature of the kinetic freezeout T kin defined from “blast wave” fit, is, on the other hand, con-stant within errors ( ± M eV ). The mean flow velocityis also constant h β i ≈ . ± . µ ( T ) = µ ch TT ch + m (1 − TT ch ) (A1)The corresponding values for nucleons at freezeout aregiven in the table, together with the thermal baryon den- TABLE V: The collision energies of RHIC in low energy scan;Fitted chemical freezeout parameters, T ch , µ Bch , and R ch , fromGrand Canonical Ensemble fit to particle yields [23], for themost central bins [0 , T kin from blast wave fit to STAR spectra [23]; Parame-ters for 2 . n Bkin . √ s (GeV) 2.4 7.7 11.5 19.6 27. 39. T ch (MeV) 65. 143.8 150.6 157.5 159.8 159.9 µ Bch (MeV) 784 398.2 292.5 195.6 151.9 104.7 R ch (fm) 9.6 5.89 6.16 6.04 6.05 6.27 T kin (MeV) 71. 116. 118. 113. 117. 117. µ Bkin (MeV) - 503. 432. 406. 363. 329. n Bkin (fm − ) 0.051 0.035 0.021 0.013 0.011 0.0082 sities at kinetic freezeout. Those are the ones used insimulations described in the main text.For completeness, let us mention that for pion thechemical potential at kinetic freezeout at T kin ≈ M eV is µ π ≈ M eV , by a curve given in Ref. [28].The modification of thermal pion spectrum at small p t in-duced by pion chemical potential (of similar magnitude)has been demonstrated already in the original paper [27],using the pion spectra from SPS NA44 experiment. Thiseffect was recently reconfirmed in LHC ALICE pion spec-tra, see Ref. [29]. Appendix B: Feed down from excited nucleon states The thermal conditions of the fireball produce, in prin-ciple, all species of hadrons with some non-zero density,most of which decay long before reaching detectors. Thisincludes excited states and resonances. Of particular in-terest to this work are the excited nucleons states whichdecay into p and n . At chemical freezeout, inclusion ofsuch states, as well as even weak decays, is necessary toget accurate fits of particle yield ratios [8].After chemical freezeout the numbers of individualspecies of hadrons are fixed as inelastic collision cease.Due to the system expanding however, the density de-creases. Decays from excited states may increase thenumber of p and n , however. The time at which theseresonances decay after a collision is still an open ques-tion and different models assuming different ordering ofevents (resonance decay then nuclei coalescence or vice-versa) give different prediction for the light-nuclei ratio N t N p /N d [3]. In this work, it has been assumed thatfeed down from these states is not significant by the timekinetic freezeout occurs and the nucleon densities usedare computed directly from T kin and µ Bkin . One shouldexpect that if such feed down substantially increases the3nucleon density at kinetic freezeout, the role of clusteringwould be more significant.At the lower temperatures of kinetic freezeout, onewould expect a reduced contribution from these > TABLE VI: Ratio of densities of nucleon excited states to thenucleon density at kinetic freezeout conditions for all N ∗ and∆ states with mass ≤ . J P is total angular momentum and parity. Listof nucleon states taken from Ref. [30]. N ∗ state J P N (1440) 12 + × − . 023 0 . N (1520) − × − . 025 0 . N (1535) − × − . 011 0 . N (1650) − × − . 005 0 . N (1675) − × − . 011 0 . N (1680) 52 + × − . 011 0 . N (1700) − × − . 006 0 . 32 + . 033 0 . 23 0 . 32 + × − . 014 0 . − × − . 006 0 . − × − . 006 0 . There is a clear √ s -dependence, with the excited stateshaving a much-reduced relative density at 2.4 GeV dueto the reduced temperature. In all cases, the lightest res-onance considered, ∆(1232) has a higher density than allother excited states considered here combined. Through-out the range √ s ∼ . − 39 GeV, the total density ofthe excited states should be about 30% of the nucleondensity. For comparison (see Table IV), the statistically-correlated clusters studied in out simulation make upabout 1% of the total configurations sampled over thesame energy range. This indicates that the total amountof feed down from the 4 N system should be compara-tively small. At √ s = 2 . Appendix C: Convergence and isotropy of theperiodic box setup Perfect periodic boundary conditions can be imposedin numerical simulation by an infinite number of boxwith image particles identical to the main simulation box. This, of course, would require infinite computa-tional power and thus one must include only a smallnumber of boxes. The question becomes then: what isthe smallest number of boxes one should include to accu-rately include the effects of the periodic boundaries? Themost important measure of this is the convergence of thedesired observables - how the output valuables vary withthe number of boxes. Here, we consider three configura-tions: a single periodic box with no images, a box withsix images - one attached to each face of the box (7 to-tal boxes), and a box enclosed by boxes touching everyface, edge, and corner (27 total boxes). The most obvious FIG. 13: (Color online) Normalized probability distributions P d ( r ) of the deuteron system with 1 box (blue • ), 7 boxes(red (cid:4) ), and 27 boxes (green N ) consideration of the effects of the periodic boxes is thecorrection to the distributions of the system by the in-teractions of nucleons with images in other boxes, as thepotential energy due to inter-box interactions affects theMetropolis updates. Additionally we look for anisotropyintroduced by the finite number of images, which breaksspherical symmetry of the system.To test these properties of the setups, we considerthe simplest system, the Walecka deuteron, in a boxsuch that the nucleon density n N = 0 . 054 fm − , thelargest density considered in the main work. To testfor anisotropy, the distribution of the nucleons’ position P ( θ, φ ), where θ and φ are the standard angles in spher-ical coordinates relative to the z-axis at the center of thebox is measured and decomposed into real spherical har-monics P ( θ, φ ) = ∞ X ℓ =0 ℓ X m = − ℓ C ℓm Y ℓm ( θ, φ ) . (C1)In the case of a perfectly isotropic distribution C ℓ,m = 0for all ℓ = 0.4 FIG. 14: (Color online) Coefficients of spherical harmonicexpansion C ℓ,m for ℓ = 1, 2 for the single-box setup (blue • )and the 7-box setup (red (cid:4) ). We find, for our 1-box and 7-box setup, all coefficientsof the expansion C ℓ,m = 0, except for C , . The fact thatthe coefficient is nearly the same for both setups suggeststhat this is an artifact of the boundaries of the box it-self rather than the images. The fact that nucleons areonly moved to the opposite side of the box when theircenter of mass crosses the boundary may cause such ef-fects. This may be due to the fact that the Y , has noprobability along the ’diagonal’ directions and the cubicgeometry differentiates the directions along the diago-nals and axes. Using a system with more images in thediagonal directions, such as the 27-box setup, may alle-viate this anisotropy, but it would increase computationtime and result in reduced statistics, so we do not do so.However, the main point of the comparison is to checkthat the 7-box setup used throughout this work does notintroduce any additional anisotropy to the system. Thisseems to be the case as the coefficients are all equal withinuncertainty in both setups. [1] P. Liu [STAR Collaboration], Nucl. Phys. A , 811(2019). doi:10.1016/j.nuclphysa.2018.10.023[2] J. Thder [STAR], Nucl. Phys. A , 320-323 (2016)doi:10.1016/j.nuclphysa.2016.02.047 [arXiv:1601.00951[nucl-ex]].[3] D. Oliinychenko, [arXiv:2003.05476 [hep-ph]].[4] K. Sun, L. Chen, C. M. Ko, J. Pu andZ. Xu, Phys. Lett. B , 499-504 (2018)doi:10.1016/j.physletb.2018.04.035 [arXiv:1801.09382[nucl-th]].[5] E. Shuryak and J. M. Torres-Rincon, Phys. Rev. C ,no. 2, 024903 (2019) doi:10.1103/PhysRevC.100.024903[arXiv:1805.04444 [hep-ph]].[6] E. Shuryak and J. M. Torres-Rincon, arXiv:1910.08119[nucl-th].[7] E. V. Shuryak and O. V. Zhirov, Nucl. Phys. B , 393(1984). doi:10.1016/0550-3213(84)90401-2[8] A. Andronic, P. Braun-Munzinger andJ. Stachel, Nucl. Phys. A , 167 (2006)doi:10.1016/j.nuclphysa.2006.03.012 [nucl-th/0511071].[9] E. Schnedermann, J. Sollfrank and U. W. Heinz, Phys.Rev. C , 2462 (1993) doi:10.1103/PhysRevC.48.2462[nucl-th/9307020].[10] D. Teaney, J. Lauret and E. V. Shuryak, Phys. Rev.Lett. , 4783 (2001) doi:10.1103/PhysRevLett.86.4783[nucl-th/0011058].[11] M. F. Herman, E. J. Bruskin, and B. J. Berne, J. Chem.Phys. , 10 (1982).[12] B. D. Serot and J. D. Walecka, Adv. Nucl. Phys. , 1(1986).[13] D. DeMartini and E. Shuryak , Many-body forces and nu-cleon clustering near the QCD critical point, in progress[14] M. Lorenz (For the HADES collaboration),Strangeness production at HADES, 3rd EMMI Work-shop (2019), [15] V. Koch, A. Bzdak, D. Oliinychenko and J. Steinheimer,[arXiv:2001.11079 [nucl-th]].[16] D. Zhang [STAR], [arXiv:2002.10677 [nucl-ex]].[17] X. Luo and N. Xu, Nucl. Sci. Tech. , no.8, 112 (2017)doi:10.1007/s41365-017-0257-0 [arXiv:1701.02105 [nucl-ex]].[18] A. Bzdak, V. Koch, D. Oliinychenko and J. Stein-heimer, Phys. Rev. C , no.5, 054901 (2018)doi:10.1103/PhysRevC.98.054901 [arXiv:1804.04463[nucl-th]].[19] M. Stephanov, Phys. Rev. Lett. , 032301 (2009)doi:10.1103/PhysRevLett.102.032301 [arXiv:0809.3450[hep-ph]].[20] M. Asakawa, S. Ejiri and M. Kitazawa, Phys. Rev. Lett. , 262301 (2009) doi:10.1103/PhysRevLett.103.262301[arXiv:0904.2089 [nucl-th]].[21] L. Adamczyk et al. [STAR], Phys. Rev. Lett. , 032302 (2014) doi:10.1103/PhysRevLett.112.032302[arXiv:1309.5681 [nucl-ex]].[22] J. Adamczewski-Musch et al. [HADES],[arXiv:2002.08701 [nucl-ex]].[23] L. Adamczyk et al. [STAR Collaboration],Phys. Rev. C , no. 4, 044904 (2017)doi:10.1103/PhysRevC.96.044904 [arXiv:1701.07065[nucl-ex]].[24] E. Shuryak and J. M. Torres-Rincon, [arXiv:2005.14216[nucl-th]].[25] M. A. Stephanov, K. Rajagopal andE. V. Shuryak, Phys. Rev. Lett. , 4816-4819 (1998) doi:10.1103/PhysRevLett.81.4816[arXiv:hep-ph/9806219 [hep-ph]].[26] B. Mohanty [STAR], J. Phys. G , 124023 (2011)doi:10.1088/0954-3899/38/12/124023 [arXiv:1106.5902[nucl-ex]].[27] C. M. Hung and E. V. Shuryak, Phys. Rev. C , 1891 (1998) doi:10.1103/PhysRevC.57.1891 [hep-ph/9709264].[28] D. Teaney, nucl-th/0204023.[29] V. Begun, W. Florkowski and M. Rybczyn-ski, Phys. Rev. C , no. 1, 014906 (2014)doi:10.1103/PhysRevC.90.014906 [arXiv:1312.1487 [nucl-th]].[30] M. Tanabashi et al. [Particle Data Group],Phys. Rev. D98