Ab initio study of nuclear clustering in hot dilute nuclear matter

We present a systematic ab initio study of clustering in hot dilute nuclear matter using nuclear lattice effective field theory with an SU(4)-symmetric interaction. We introduce a method called light-cluster distillation to determine the abundances of dimers, trimers, and alpha clusters as a function of density and temperature. Our lattice results are compared with an ideal gas model composed of free nucleons and clusters. Excellent agreement is found at very low density, while deviations from ideal gas abundances appear at increasing density due to cluster-nucleon and cluster-cluster interactions. In addition to determining the composition of hot dilute nuclear matter as a function of density and temperature, the lattice calculations also serve as benchmarks for virial expansion calculations, statistical models, and transport models of fragmentation and clustering in nucleus-nucleus collisions.

The equation of state of nuclear matter and its composition are topics that lie at the heart of nuclear physics and determine important properties of neutron star evolution, binary neutron star mergers, core-collapse supernovae, and other astrophysical processes.Many terrestrial studies have used particle yields and cluster distributions from heavy-ion collisions to infer the equation of state and composition of nuclear matter as a function of density, temperature, and isospin [1][2][3][4][5][6][7][8][9][10][11][12].In this work, we focus on the properties of hot dilute nuclear matter at densities below saturation and temperatures below the pion production threshold.This regime is relevant to the liquid-gas phase transition in nuclear matter [13] and corecollapse supernovae [14,15].While previous ab initio lattice calculations have determined the equation of state and phase diagram [16], in this work we present the first ab initio calculations of the cluster composition of hot dilute nuclear matter.
Nuclear lattice effective field theory (NLEFT) [38,39] is a powerful numerical method which uses finite volume Monte Carlo (MC) methods and is formulated in the framework of nuclear chiral EFT [40].As the NLEFT simulations include many-body correlations to all orders, effects such as deformation and clustering appear naturally.NLEFT defines a systematically improvable ab initio nuclear theory starting from bare nuclear forces and has favorable computational scaling with the number of nucleons.NLEFT has been successfully applied to probe alpha clustering in finite nuclei [41][42][43][44][45][46][47][48].In Ref. [16], it was extended to calculations of the thermodynamics of nuclear systems using the pinhole trace algorithm (PTA).
In the present work, we use NLEFT simulations to determine the fractions of nucleons forming two-body, three-body, and four-body (alpha) clusters in symmetric nuclear matter as a function of temperature T and density ρ.We use the Wigner SU(4)-symmetric nuclear force introduced in Ref. [49], where the interaction is independent of spin and isospin.While the extension to high-fidelity chiral nuclear forces is straightforward, the calculations require new algorithms to perform perturbation theory for wave functions currently being developed [50] and will therefore be presented in future work.It should also be noted that SU(4)-symmetric forces already provide a highly successful description for the ground-state properties of many light and medium-mass nuclei [49], for pure neutron matter and symmetric nuclear matter [16] as well as for the low-lying excited states of the 12 C nucleus [47,48], where alpha clustering effects are extremely important.Details about the Wigner SU(4)-symmetric nuclear interactions and the PTA can be found in the Supplemental Material [51] and Refs.[16,49].We note that the SU(4)-invariant deuteron is degenerate with the di-neutron and di-proton ground state and has less than half of the physical deuteron binding energy.The SU(4) symmetry makes the enumeration of light clusters simple and therefore useful for introducing the light-cluster distillation method.
We regard the clustering of nucleons in nuclear matter as a spatial localization of nucleons.As a detailed probe of the nucleon correlations, we compute the following set of correlation functions: Here, n is a site on a cubic spatial lattice with lattice spacing a and volume L 3 , n is the radial distance from the origin, the columns :: denote normal ordering, and is the nucleon density operator for the spin σ i and isospin τ i .By {σiτi} , we denote summation over all spin-isospin indices without repetition of indices.This avoids the strong short-range exclusion effects due to the Pauli principle.The correlation functions in Eq. ( 1) are computed using the rank-one operator method introduced in Ref. [50].To illustrate the characteristic features of the correlation functions given in Eq. (1) we consider a system with A = 16 nucleons in a cubic box of length L = 10 corresponding to a density ρ = 0.007 fm −3 and perform lattice simulations for temperature T = 10 MeV.In Fig. 1 the lattice results are shown by the solid lines and red-square points.We find that all the correlation functions are nearly flat for large spatial separations and have a strong peak at shorter distances.As we will now see, this signal is consistent with a gas of small clusters.
In order to quantitatively determine the cluster abundances, we use the method of light-cluster distillation.We write the correlation functions for hot dilute nuclear matter as a weighted sum of correlation functions for individual light clusters, Here we narrow our focus to light clusters such as dimers ( 2 H, pp, nn), trimers ( 3 H, 3 He) and alpha particles ( 4 He).As we consider an SU(4)-invariant Hamiltonian and symmetric nuclear matter, we need not distinguish different types of dimers and trimers.In Eq. ( 2), ⟨• • • ⟩ k denotes the computed correlation functions for dimers (k = 2), trimers (k = 3), and alphas (k = 4) in their ground states.The semi-positive valued parameters w k are interpreted as the number of clusters and determined by a least-square fitting procedure.Finally, G l ij denotes the long-range parts of the correlation functions defined by, e.g., where n m = 8 fm/a separates the short-range and long-range physics.For sufficiently large values of n m , the results are independent of n m .Although the character of the density correlations at short distances depend on the regulator used for the nuclear interactions, the distillation coefficients into clusters should be independent of the regulator when the low-energy nuclear interactions remain the same.In Fig. 1, the results from fitting lattice data to Eq. ( 2), shown by the dashed lines and blue-circle points, indicate that the correlation functions are accurately represented by lightcluster distillation with the three parameters, w 2 , w 3 and w 4 .This finding suggests that the system is a gas of nucleons and light clusters.After determining the w k parameters from light-cluster distillation, we compute the mass fractions given by X k = kw k /A for clusters composed of k nucleons as well as the total fraction of nucleons bound in light clusters defined as We perform lattice simulations at temperatures ranging from 10 to 20 MeV for densities ρ = 0.007 and 0.014 fm −3 , and we employ light-cluster distillation.The heavier-cluster mass fraction with A > 4 is estimated to reach about 10% at density 0.01 fm −3 and temperature 10 MeV [17].The abundance of heavier clusters with A > 4 are much smaller for lower densities and/or higher temperatures.Correlation functions for ρ = 0.007 and 0.014 fm −3 at several representative temperatures are shown in the Supplemental Material [51].
In Fig. 2, we show the mass fraction and total fraction results for ρ = 0.007 fm −3 and 0.014 fm −3 .We find that the fraction of nucleons bound in light clusters decreases monotonously with increasing T .This observation strongly suggests the thermal evaporation of nucleons from the clusters.On the other hand, the mass fractions of dimers, trimers and alphas show somewhat complementary behaviors.For fixed ρ, X 3 and X 4 decrease with temperature from 10 MeV to 20 MeV, while the behavior for X 2 is slightly increasing or relatively flat.For fixed T , X 3 and X 4 increase with density from 0.007 fm −3 to 0.014 fm −3 , while X 2 is decreasing.
We have also perform lattice simulations at the same temperatures and densities using a larger volume of length L = 12.We find that the mass fraction results are nearly the same, and the results are shown in the Supplemental Material [51].This indicates that the finite size effects are small and the lattice results are close to the thermodynamic limit.This is not unexpected since the box length L is significantly larger than the thermal wavelength λ = 2π/(mT ), which is a measure of the correlation length at non-zero temperature.
As seen in Fig. 2 at lower temperatures and higher densities, the uncertainties on the determination of the light-cluster mass fractions become larger.This is caused by the increasing proportion of clusters composed of more than four nucleons.The determination of cluster abundances at lower temperatures and higher densities requires enumerating clusters with A > 4.This can be performed by calculating and analyzing correlation functions beyond those in Eq. (1) and will be addressed in future studies.See the Supplemental Material for further discussion of heavier clusters [51].
In Fig. 3, we compare our results for X k (ρ, T ) with an ideal gas model, where the system is assumed to consist of free nucleons, dimers, trimers, and tetramers (alpha clusters), governed by the grand canonical ensemble [51,52].
To ensure clarity and avoid any potential ambiguities arising from the presence of heavier clusters, we focus on densities ρ ≤ 0.007 fm −3 in our analysis.Notable deviations become apparent around ρ ≃ 0.007 fm −3 for all types of clusters.Nevertheless, at lower densities, specifically around ρ ≃ 0.001 fm −3 , the agreement is excellent.These observations suggest that the deviations at higher ρ arise from nonnegligible cluster-nucleon and cluster-cluster interactions.
In order to quantitatively clarify such effects, we perform lattice simulations at temperature T = 5 MeV and for various numbers of nucleons A in a cubic box of length L = 10.
For each A, we also vary the proton number Z from 0 to A/2 to step-by-step activate cluster-nucleon and cluster-cluster interactions.The results are presented in Fig. 4. The ideal gas results are predictions from the canonical ensemble with the same particle content and physical volume (for more details, see [51]).
We first consider the dimer case.For A ≤ 4 and Z ≤ 1, the deviations between lattice and ideal gas results are negligible, which implies that nucleon-nucleon, nucleon-dimer, and dimer-dimer (nn-nn and nn-np) interactions have no significant impact on the dimer abundance.For A = 4 and Z = 2, lattice results give slightly lower abundances than that of the ideal gas.This is suggestive of dimer-dimer (nn-pp or npnp) interactions reducing the dimer abundance.The observed discrepancies become significant as the values of A and Z increase.
For the case of the trimer, noticeable differences first appear for A = 4, which suggests that the trimer abundance is increased due to nucleon-trimer interaction.We also observe a similar behavior in the case of the tetramer or alpha, and we find that the nucleon-alpha interactions similarly enhance the alpha abundance.From the consistency of the results in Fig. 3 and Fig. 4, we infer that the effects of clusternucleon and cluster-cluster interactions are responsible for the observed deviations from the ideal gas abundances.
In summary, we have used NLEFT simulations to investigate nuclear clustering in hot dilute nuclear matter.We have introduced a theoretical framework called light-cluster distillation to determine the abundances of light nuclear clusters.The corresponding abundances of clusters are well determined for the cases where heavier cluster abundances are negligible.The comparisons with ideal gas predictions show excellent agreement at low densities ρ ≃ 0.001 fm −3 .At higher densities, we observe deviations from ideal gas abundances due to cluster-nucleon and cluster-cluster interactions.Since we perform fully non-perturbative ab-initio calculations starting from bare nucleonic interactions, our results contain manybody correlations to all orders and can thus serve as benchmarks for approximate calculations and models such virial expansions, statistical models, and transport models.In future work, we will perform analogous studies using high-fidelity chiral nuclear forces combined with state-of-art many-body methods such as perturbative quantum Monte Carlo [53] together with the rank-one operator method [50] and wave function matching [54].

SUPPLEMENTAL MATERIALS Lattice Hamiltonian with Wigner SU(4)-symmetric nuclear force
We define the Hamiltonian in a cubic box with lattice coordinates n = (n x , n y , n z ).We employ an SU(4)-symmetric interaction, and our Hamiltonian reads where K is the kinetic term with nucleon mass m = 938.9MeV, and the colons indicate normal ordering.The smeared density operator ρ(n) is defined as where i is a collective spin-isospin index, and the smeared creation and annihilation operators are defined as The summation over the spin and isospin implies that the interaction is SU(4) invariant.Here, the nearest-neighboring smearing parameter s L controls the strength of the local part of the interaction, while the nearest-neighboring smearing parameter s N L controls the strength of the nonlocal part of the interaction.The couplings C 2 and C 3 give the overall strength of the two-body and three-body interactions, respectively.In this work, we perform our calculations at a lattice spacing a = 1.32 fm, which corresponds to a momentum cut-off Λ = π/a ≈ 471 MeV, and we utilize averaged twisted boundary conditions [1].The parameters values of our interaction are set as s L = 0.061, s N L = 0.5, C 2 = −3.41× 10 −7 MeV −2 , and C 3 = −1.4× 10 −14 MeV −5 .These parameters are proposed in Ref. [2] by fitting the properties of light and medium-mass nuclei.Moreover, it has been observed that this specific set of parameters effectively captures the thermal properties exhibited by nuclear matter [1].

Pinhole trace algorithm
In canonical ensemble calculations, we fix the parameters of nucleon number A and temperature T , yielding the expectation value of any observable O as follows, where Z(β) is the canonical partition function, β = T −1 is the inverse temperature, H is the Hamiltonian, and Tr A is the trace over the A-nucleon states.Throughout this work, we use units where ℏ = c = k B = 1.In Ref. [1], the pinhole trace algorithm (PTA) is proposed to evaluate the Eq.(S4).The explicit form of the canonical partition function Z(β) is expressed in the single particle basis as follows, where the basis states are Slater determinants composed of A nucleons.c i = (n i , σ i , τ i ) denotes the quantum numbers of the ith nucleon, where σ i is the spin and τ i is the isospin.On the lattice, the components of n i take integer numbers from 0 to L − 1, where L is the box length in lattice units.The neutron number N and proton number Z are separately conserved, and the summation in Eq. (S5) is limited to the subspace with specified values of N and Z.The inverse temperature β is divided into L t slices with temporal lattice spacing a t such that β = L t a t , where a t = (2000 MeV) −1 is taken in this work.We use the Hubbard-Stratonovich transformation [3,4] and a discrete auxiliary field which couples to the nucleon density, where the ω k 's and ϕ k 's are real numbers, and we require ω k > 0 for all k.Taking D = 3 and expanding Eq. (S6) up to O(ρ 4 ), the constants ϕ k 's and ω k 's can be determined by comparing both sides order by order, see Ref. [2] for details.Now we can rewrite partition function Z(β) in terms of the transfer matrix operator with discrete auxiliary field, where is the normal-ordered transfer matrix, and s nt is our shorthand for all auxiliary fields at n, n t [5,6].Note that for the decomposition in Eq. (S6), the s nt can only take several discrete values {ϕ k ln ω k }, so the path integral over real variables means the summations over indices.For a given configuration s nt , the transfer matrix M (s nt ) consists of a string of one-body operators which are directly applied to each single-particle wave function in the Slater determinant.Similar to the previous work [1], the abbreviations To evaluate the Eq.(S7), the projection Monte Carlo methods with importance sampling are used to generate an ensemble Ω of ⃗ s, ⃗ c of configurations according to the relative probability distribution The expectation value of any operator Ô can be expressed as where To generate the ensemble Ω, ⃗ s and ⃗ c are updated alternately.First for a fixed nucleon configuration ⃗ c the auxiliary fields ⃗ s are updated through the shuttle algorithm [2], then the nucleon configuration ⃗ c is updated using the Metropolis algorithm.The details for updating ⃗ s and ⃗ c can be found in Ref. [1].

Correlation functions
In this section, we present the correlation functions G 11 , G 21 , G 31 , and G 22 alongside the results from light-cluster distillation for various temperatures and densities.Figs.S1 and S2 depict the results for a density of ρ = 0.007 fm −3 at temperatures of T = 14 and 18 MeV, respectively.Also, Figs.S3, S4, and S5 show the results for a density of ρ = 0.014 fm −3 at temperatures of T = 10, 14, and 18 MeV, respectively.At T = 10 MeV and ρ = 0.014fm −3 , we observe a noticeable decline in the agreement of light-cluster distillation for G ij due to the non-negligible fractions from heavier clusters with A > 4. As the temperature increases, the contribution of heavier clusters becomes less significant, resulting in an improved agreement between light-cluster distillation and the correlations.
Considering that heavier clusters induce spatial localizations of two nucleons with the same spin and isospin, they can be probed by examining the correlation function στ : ρ στ (0)ρ στ (r) :.Fig. S6 illustrates this correlation as a function of distance at three different temperatures.The correlation function is zero at the origin due to the Pauli principle and increases up to its maximum value at around r = 3 fm.For ρ = 0.007 fm −3 , at all temperatures the correlation function remains nearly constant for larger r.However, for ρ = 0.014 fm −3 , a significant peak is observed at temperature of T = 10 MeV, which provides a direct evidence for the emergence of heavier clusters.The diminishing importance of heavier clusters with increasing temperatures is further substantiated by the attenuation of this peak.

Finite volume effects
To check the impact of finite volume effects on the mass factions X k 's, a larger box size of L = 12 is used in the calculations.Total nucleon numbers are adjusted to A = 28 and 56 in order to achieve the densities of 0.007 and 0.014 fm −3 , respectively.In Fig. S7, a  density ρ = 0.007 fm −3 , almost perfect consistency is observed for all shown temperatures.Some discrepancies of the center values are observed for ρ = 0.014 fm −3 at temperatures of 10 and 12 MeV due to the heavier clusters.Nevertheless, good consistency is still found with increasing temperature, indicating that the obtained X k values from light-cluster distillation are quite reliable once the fractions of heavier clusters are negligible.

Ideal gas estimated from the grand canonical ensemble
The system is approximated as an ideal gas composed by noninteracting nucleons, dimers, trimers, and tetramers (alphas).We simply denote the corresponding chemical potentials as µ 1 , µ 2 , µ 3 , and µ 4 with considering that the interactions in our lattice calculations are SU(4) invariant and the calculations are restricted to symmetric nuclear matter.In the thermodynamic equilibrium with a temperature T , the average number for each species of particle are described via, where g k is the degeneracy from spin and isospin.For nucleons, g 1 = 4.For dimers, g 2 = 6, where 2 from nn and pp dimers, and 4 from np dimers.For trimers, g 3 = 4 from nnp and npp trimers.This can be seen from Fig. S8, where E coll.kin.(p)'s on the lattice are obtained from boosted ground-state wave functions.The summation over energy in Eq. (S12) can be transferred into integration over momentum and spatial spaces,

FIG. 1 .
FIG.1.Correlation functions in Eq. (1) versus relative distance r in units of fm.The solid lines and red squares show our results for the density ρ = 0.007 fm −3 (total particle number A = 16 for the box size L = 10) and temperature T = 10 MeV.The dashed lines and blue circles show the fitted results of light-cluster distillation using Eq.(2).The lines connecting the data points are intended to guide the eye.

FIG. 3 .
FIG. 3. The mass fraction X k for temperatures T = 14 (top) and 10 MeV (bottom) as a function of density.Results are shown for A = 16, and L = 10 . . .20 for steps of 2. The lattice results are shown by squares, circles, and triangles, while the dashed lines are the predictions from an ideal gas model.The error bars reflect the fitting uncertainty for the w k in light-cluster distillation.

FIG. 4 .
FIG. 4. The number of dimers (top panel), trimers (middle panel), and alphas (bottom panel) as a function of A and proton number Z.The lattice simulations are performed at T = 5 MeV and for L = 10.Circles with dashed lines show the predictions of an ideal gas model in the canonical ensemble.For each A, the proton number Z is varied from 0 (left) to A/2 (right).
FIG.S1.Correlation functions G11, G21, G31, and G22 for relative distance r in units of fm.The solid lines and red squares show lattice MC results at temperature of T = 14 MeV using A = 16 and L = 10 which corresponds to the density ρ = 0.007 fm −3 .The dashed lines and blue circles show the fitted results of light-cluster distillation.The lines connecting the data points are intended to guide the eye.
FIG.S3.Correlation functions G11, G21, G31, and G22 for relative distance r in units of fm.The solid lines and red squares show lattice MC results at temperature of T = 10 MeV using A = 32 and L = 10 which corresponds to the density ρ = 0.014 fm −3 .The dashed lines and blue circles show the fitted results of light-cluster distillation.The lines connecting the data points are intended to guide the eye.

For alpha, g 4 = 1 . 2 2km + a 2 p 2 2km 2 .
FIG.S5.Correlation functions G11, G21, G31, and G22 for relative distance r in units of fm.The solid lines and red squares show lattice MC results at temperature of T = 18 MeV using A = 32 and L = 10 which corresponds the density ρ 0.014 fm −3 .The dashed lines and blue circles show the fitted results of light-cluster distillation.The lines connecting the data points are intended to guide the eye.