Symmetry Energy and the Pauli Exclusion Principle

: In this article we present a classical potential that respects the Pauli exclusion principle and can be used to describe nucleon-nucleon interactions at intermediate energies. The potential depends on the relative momentum of the colliding nucleons and reduces interactions at low momentum transfer mimicking the Pauli exclusion principle. We use the potential with Metropolis Monte Carlo methods and study the formation of ﬁnite nuclei and inﬁnite systems. We ﬁnd good agreement in terms of the binding energies, radii, and internal nucleon distribution of ﬁnite nuclei, and the binding energy in nuclear matter and neutron star matter, as well as the formation of nuclear pastas, and the symmetry energy of neutron star matter. Author Contributions: Conceptualization, C.O.D. and G.F.; methodology, C.O.D. and G.F.; software, C.O.D. and G.F.; validation, C.O.D., G.F. and J.A.L.; formal analysis, C.O.D., G.F. and J.A.L.; inves-tigation, C.O.D., G.F. and J.A.L.; resources, C.O.D. and G.F.; data curation, G.F.; writing—original draft preparation, J.A.L.; writing—review and editing, C.O.D., G.F. and J.A.L.; visualization, G.F.; supervision, C.O.D.; project administration, C.O.D.; funding acquisition, All authors


Antecedents
In 1953, Bethe stated [1] that the Nucleon-Nucleon (NN) interaction was the problem most studied in the history of the world, and almost 70 years later the statement is still true. Since the discovery of the neutron by Chadwick in 1932, many experiments, labor and mental works, have been devoted to study the most fundamental problem in nuclear physics, the NN interaction. The goal of describing the atomic-nuclei properties in terms of the interactions between pairs of nucleons is the main objective of nuclear physics, unfortunately, because of the complexity of the quantum mechanical problem, no complete theory has been developed, and progress has been based on approximations.
So far it has been a trade-off; while most researchers have decided to study nuclear systems using mean field methodologies that only respect quantum mechanics approximately, a few others have overlooked some quantum aspects to use classical dynamics that preserves nucleon-nucleon correlations, statistical fluctuations, clusterization phenomena, phase changes, and critical phenomena; all features of the upmost importance in the later stages of fragmentation reactions, and that are not present in the quantum models.

Quantum Caveats
Quantum effects affect the behavior of many body systems on, at least, two fronts: Energy distribution in bound systems and wave mechanics.
From the point of view of the energy, in bound clusters the energy of individual nucleons becomes discrete and the distribution of energy levels is ruled by Fermi-Dirac statistics with the Pauli exclusion principle further regulating the occupation of such levels. Although for high excitation energies the phase space available for nuclei is ample enough to render Pauli blocking practically obsolete [44], at lower energies it is bound to have an impact on the nucleon dynamics prohibiting energy and momentum transfer in some nucleon-nucleon collisions; such a problem is expected to be relevant whenever the number of quantum states available to a nucleon at a given temperature is comparable to the number of nucleons.
A second problem of the validity of the classical approach is connected to the wave nature of the nucleons. Wave mechanics leads to classical mechanics as soon as the distance between particles is larger than their de Broglie wavelengths. This poses a problem when using classical dynamics at high densities, but not at the low sub-saturation densities [44].
Other quantum effects, such as collective excitations, superfluidity, and superconductivity, etc., occur at much lower temperatures and are unlikely to play a major role in the dynamics of medium energy nuclear collisions.
In intermediate-energy heavy-ion collisions, the nuclear systems will first be compressed and heated up to temperatures of several MeV, and then expanded and cooled down as they undergo a phase transition into a liquid-gas mixed phase. Since the first stage occurs at high energies, Pauli blocking is not expected to play a major role in the nucleon-nucleon interactions. In the second stage, however, as the system cools and expands, Pauli exclusion will regulate the final-stage interactions. Likewise, in cold stable nuclei, the occupacy exclusion in energy levels causes nuclei in their ground state to have a non-zero kinetic energy known as Fermi motion.

Pauli Blocking
To mimic Pauli blocking in classical collisions one must introduce a way to regulate the nucleon-nucleon collisions at low energies and low momentum transfer. A plausible approach to this problem is to introduce a momentum-dependent blocking or repulsion to simulate the Pauli exclusion principle. Such an approach was first taken by Wilets et al. [29], and followed up by Dorso and Randrup, who improved the approximation to the nucleon phase-space distribution. Some other models were based on similar approaches with density-dependent potentials and quasiparticles [47,48].
In the first step, Dorso and Randrup used a momentum-dependent repulsive potential to simulate the Pauli exclusion in a nuclear-like Fermi gas over the temperatures and densities attained in nuclear collisions at intermediate energies [49]. To prevent nucleons identical in spin and isospin getting too close in phase space, a repulsive, momentumdependent two-body potential was used.
As a next step, a modified Lennard-Jones-like two-body interaction was added to mimic nuclear interactions [30], along with a standard Coulomb repulsion between protons. The resulting model yielded a system with the proper saturation energy and compressibility of nuclear matter over a broad range of temperatures and densities, as well as the proper sizes and binding energies of finite nuclei [50].
Although the potentials of [50] were very promising in simulating the Pauli exclusion principle and endowing classical models with a realistic ground-state Fermi motion, the models lack a proper nucleon-nucleon repulsion. Specifically, the model allowed the formation of dineutron structures that are not to be expected, despite the short time period during a collision. Thus the motivation of this study: To develop a model that maintains the advantages of classical models but, at the same time, respects the Pauli exclusion principle while preventing the formation of nonphysical states. The structure of the article is as follows: Section 2 presents the new potential introducing its various components. This is followed by applications to model the finite nuclei in Section 3, and infinite systems in Section 4. The model is then used to obtain the symmetry energy of neutron star matter in Section 5. A discussion of the results in Section 6 closes the article.

The Model
To mimic the Pauli exclusion principle, we designed potentials following the procedure of Dorso and Randrup [49]. In a nutshell, a dimensionless distance in phase space was defined as s 2 ij = p 2 ij /p 2 0 + r 2 ij /r 2 0 , where, for any pair of nucleons i and j, r 2 = | r i − r j | 2 , and p 2 = | p i − p j | 2 , and with the parameters p 0 and q 0 determining the size of the effectively excluded volume around each particle in phase space, and adjusted such that the uncertainty relation, p 0 q 0 ≈ 2h, is satisfied.
Modifying the Dorso and Randrup potential by trial and error we arrived at: where the scale factors are V P = 10.33 MeV, p 0 = 61.969 MeV/c, and q 0 = 6.0 fm, and r c = 9.0 fm is a cut-off distance. In our case we find that p 0 q 0 = 1.88h, very close to the uncertainty relation, and in agreement with the expectation that the product be slightly smaller than 2h due to geometric considerations [49]. Furthermore, the present potential improves over the Dorso and Randrup potential by forbidding the formation of di-neutrons. Notice that the added modulating term imposes a repulsive force between equal nucleons at short values of p and r, basically forbidding their interaction. This repulsion ceases for values of r > r c , i.e., the potential is truncated at r = r c . Also notice that for r = 0 and p = 0, the potential becomes V nn = V pp = V p . The potential (1) generates a kinetic energy (i.e., a "Fermi energy") that is in agreement with the results obtained by Dorso and Randrup in [49]. Figure 1 shows the potential for the case of zero momentum transfer, p = 0. The next step is to incorporate a potential to describe the nuclear force. For the interaction potential between two nucleons, we use: And for the Coulomb potential, we use the point-charge repulsion: The parameters of the nuclear potential are V np = 34.8 MeV, V NN = 16.0 MeV, and r c = 5.4 fm, and the rest were estimated based on the experimental value of the nuclear binding energy of the 4 He nucleus: where a is the radius of an α particle, a = 2.3 fm, and σ is twice the proton radius, σ = 1.625 fm. The compressibility (K = 9ρ 2 0 d 2 V/dρ 2 ) yields the condition n × m = − K E(ρ 0 )/A , which, for E(ρ 0 ) = −16 MeV, and compressibility K = 288 MeV, is satisfied by n = 6 and m = 3. The terms V np (r c ), V np (r c ), and V Coulomb (r c ) are the values of their respective accompanying terms at r c to make the potentials smoothly go to zero and avoid the extraneous force that an  Figure 1. The upper blue line corresponds to the spatial factor of the Pauli potential, V q = V Pauli (r, p = 0), plotted as a function of r in units of q 0 = 6 fm; notice the truncation at r = 9 fm. The lower orange curve is the exponential reduction imposed by the momentum dependent factor, V p = exp(−p 2 /2p 2 0 ), as it goes from p = 0 to p = 2p 0 .

Finite Nuclei
The model created with the potentials (1)-(3), can be used to construct finite nuclei. Since the potentials are not separable, symplectic integrators of the equations of motion cannot guarantee the conservation of certain quantities, such as the energy, and the behavior of such systems is best calculated using Metropolis-Monte Carlo methods (MMC); see Appendix A for a synoptic review of the method.
Nuclei were prepared by caging nucleons in a potential well. The systems were then cooled down employing an MMC procedure, until reaching the "ground state" (say, T = 0.01 MeV), where they were kept in thermal equilibrium through a heat bath. Figure 2 shows the radii obtained for several nuclear-like clusters bound with the potentials (1)-(3), calculated with a Metropolis-Monte Carlo method [51], and compared to the experimental values of the radii. The radii were computed as the r.m.s value of the nucleons positions in the ground state. The curve R = c A 1/3 corresponds to the best fit of the radii as a function of A 1/3 (see caption for details).
Likewise, Figure 3 shows the binding energy per nucleon obtained for the nuclei of Figure 2. Although the trend follows the same trajectory of the experimental values, the simulated nuclei are about 50% over-bound, close to the results of the Simple Semiclassical Potential of Horowitz and coworkers [33], see Figure 2 in [52] for a comparison with other models.  To get an idea of the internal structure of the clusters, Figure 4 shows the average distance between nucleons for several nuclei in the ground state (at T = 0.01 MeV). Also shown in the insets, is the un-normalized pair correlation function for np, nn, and pp of 4 He and 137 Cs.

Infinite Nuclear Matter
The model presented by Equations (1) and (2), can also be used to construct infinite systems; notice that the Coulomb interaction, Equation (3), is excluded. This is achieved by placing a large number of nucleons in a cubic cell with 26 replicas surrounding the central cell. Once again, utilizing the MMC method, systems of varying density and isotopic content can be constructed at specific temperatures, and used to determine, e.g., energydensity isotherms, pasta-like structures (at low temperatures), and the corresponding symmetry energy. Infinite systems can be constructed without and with an embedding electron gas, which are known as nuclear matter and neutron star matter, respectively.

Nuclear Matter
Systems composed solely of uncharged protons (otherwise the system is thermodynamically stable) and neutrons are known as nuclear matter. Figure 5 shows the energy per nucleon obtained for nuclear matter systems constructed with 1000 nucleons with equal number of protons and neutrons (with their corresponding spin), for temperatures ranging from T = 0.01 MeV to 1.0 MeV. The minimum of the peculiar "∪" shape corresponds to the saturation density, which is observed to be close to ρ = 0.18 fm −3 for σ = 1.625 fm in Equation (2). The flattening of the energy-density isotherms at sub-saturation densities indicates a departure from the liquid-like phase that exists around saturation density to a mixed liquid-gas phase; at T = 0.10 MeV the state corresponds to a frozen medium, and to a pasta structure at low densities.
The plot of the energy per nucleon as a function of the temperature can indicate changes of phase. Figure 6 shows how the energy per nucleon of a system at a fixed density varies with the temperature. The systems shown were constructed with 5832 nucleons with equal numbers of protons and neutrons and with densities in the range ρ = 0.04 to 0.14 fm −3 . Notice that the system with ρ = 0.12 fm experiences a sharp change around T ≈ 1.5 MeV, denoting a change of phase at that density; previous studies with other potentials [31] connect these discontinuities with a liquid-to-solid change of phase.
Using the present model, structures of nuclear matter with 1,000 nucleons were found at sub-saturation densities and low temperature (T = 0.01 MeV), see Figure 7. Structures at ρ = 0.14 and 0.12 fm −3 depict continuous liquid-like arrangements, while that at ρ = 0.08 fm −3 show pasta-like structures. Figure 8 shows the structure obtained with 5,832 nucleons, equal number of protons and neutrons, at ρ =0.04 fm −3 and T = 0.8 MeV; the structure on the left shows the simulation cell and its replicas.

Neutron Star Matter
Similarly, the model can be used to study a medium generically known as "Neutron Star Matter". It should be noted that crusts of neutron stars, i.e., the outermost layer of about 1 km of depth, are composed of neutrons, protons and electrons with a varying ratio of protons and neutrons; here we study this type of systems with different ratios without claiming that they exist in actual neutron stars. The medium is governed by means of equations (1) and (2), and an all embedding electron gas.
The electron cloud introduces an screening effect on the Coulomb potential of the The snapshot on the left shows the simulation cell and its images. Detailed snapshots of the simulation cell can be seen on the right.

Neutron Star Matter
Similarly, the model can be used to study a medium generically known as "Neutron Star Matter". It should be noted that crusts of neutron stars, i.e., the outermost layer of about 1 km of depth, are composed of neutrons, protons, and electrons with a varying ratio of protons and neutrons; here we study this type of systems with different ratios without claiming that they exist in actual neutron stars. The medium is governed by means of Equations (1) and (2), and an all embedding electron gas.
The electron cloud introduces a screening effect on the Coulomb potential of the protons, which in turn modifies the pasta structures. Such screening between the protons and the electron gas can be treated by means of a Thomas-Fermi screening potential of Symmetry 2021, 13, 2116 9 of 20 the form V TM (r) = e 2 r e −r/λ − V rc , with a screening length λ and V rc = e 2 r c e −r c /λ . The effect of the screening on the formation of the pasta has been studied extensively [31,41,60], the value we use, λ = 10 fm, was found to be the minimal length at which the effects of the finite cell size are eliminated and the pasta morphology stabilizes and ceases to depend on λ. Likewise, the cutoff distance, r c , is set to the value of the screening length, r c = 10 fm. Figure 9 shows the energy per nucleon as a function of the density for an isospin symmetric system with 5832 nucleons (see caption for details). The dashed lines correspond to the a quadratic fit of data between ρ = 0.13 fm −3 and 0.19 fm −3 . The resulting compressibility for neutron star matter turns out to be κ = 523.3 MeV and 572.3 MeV (not to be confused with that for nuclear matter which is about 288 MeV).  Figure 10 shows the energy per nucleon as a function of the temperature for values of the density between ρ = 0.04-0.14 fm −3 , for a system with 5832 nucleons with equal numbers of protons and neutrons (with the corresponding spins). Notice that the energy at ρ = 0.12 fm −3 experiences a sharp change around T ≈ 1.5 MeV, denoting a change of phase. It should be remarked that, at a difference from molecular dynamics calculations, MMC calculations do not explore the energy landscape profusely, and arrive at configurations of minimum energy that correspond to meta-stable liquid states which, at slightly lighter densities drop in internal energy corresponding to non-uniform states; these changes appear as sudden drops in the E − ρ plots at around ρ ≈ 0.13 fm −3 , and tend to diminish for T 0.1 MeV, and are less noticeable in nuclear matter (cf. Figure 5).
The behavior of E can be examined in terms of the kinetic and potential energies. Figures 11 and 12 show, respectively, the kinetic and potential components of the energy of Figure 10 as a function of the temperature. While the potential energy does not show a significant change at around T ≈ 1.5 MeV, the kinetic energy does show a change in the slope for low densities, especially for ρ = 0.12 fm −3 . Figure 11 also shows that the potential energy changes its trend at T ≈ 3 MeV, presumably due to the Pauli potential at low and high momentum transfer.  To visualize changes in the internal structure we recur to the Kolmogorov statistic, which computes the discrepancy between the sampled spatial distribution of the nucleons and a uniform distribution. This statistic, plotted in Figure 13, shows changes at 1.5 T 3 MeV for densities 0.11 fm −3 ρ 0.12 fm −3 , denoting the changes from uniformity to pasta-like structures at low densities, but not at ρ ≥ 0.13 fm −3 .  Figure 14 depicts better the phenomenon by plotting the "net" kinetic energy, after subtracting the linear contribution in Figure 11 (shown in gray). The sharp slopes in Figure 14 correlate with sharp topological changes in the system, as can be noticed from successive snapshots around T ≈ 1.4 MeV (see Figure 15). The divergence from uniformity can also be seen in Figure 15, which shows the structures formed in systems with ρ = 0.12 fm −3 , with a "bubble" appearing in the simulation cell appearing at T = 1.4 MeV.  Figure 14 depicts better the phenomenon by plotting the "net" kinetic energy, after subtracting the linear contribution in Figure 11 (shown in gray). The sharp slopes in Figure 14 correlate with sharp topological changes in the system, as can be noticed from successive snapshots around T ≈ 1.4 MeV (see figure 15).
The divergence from uniformity can also be seen in Figure 15, which shows the structures formed in systems with ρ = 0.12 fm −3 ), with a "bubble" appearing in the simulation cell appearing at T = 1.4 MeV.

Pastas of Neutron Star Matter
The structures of neutron star matter at low temperatures are shown in Figure 16. Again, the structures appear homogeneous at ρ = 0.16 fm −3 , while the first bubble can be seen at 0.12 fm −3 . At ρ = 0.08 fm −3 a "pasta-like" structure is formed. Examining the low density region more closely, Figure 17 shows the pseudo-pasta structures obtained at ρ = 0.04, 0.08 and 0.10 fm −3 for a system with 2,916 protons and 2,916 neutrons.

Pastas of Neutron Star Matter
The structures of neutron star matter at low temperatures are shown in Figure 16. Again, the structures appear homogeneous at ρ = 0.16 fm −3 , while the first bubble can be seen at 0.12 fm −3 . At ρ = 0.08 fm −3 a "pasta-like" structure is formed. Examining the low density region more closely, Figure 17 shows the pseudo-pasta structures obtained at ρ = 0.04, 0.08, and 0.10 fm −3 for a system with 2916 protons and 2916 neutrons.

Symmetry Energy
The nuclear symmetry energy became a topic of intense investigations when radioactive beam facilities began to study nuclei away from the stability valley; understanding isospin asymmetric nuclear matter is needed in areas of nuclear physics and astrophysics [22].
In 1935, Weizsäcker introduced an asymmetry term in a parametrization of the nuclear binding energy to provide favorable binding to those nuclei with a similar number of protons than neutrons [62]. Such a term was later modified to include the role of isospin in a generalized density-dependent asymmetry term [63]. For the case of microscopic models, the symmetry energy can be evaluated numerically following the procedure outlined in Appendix B, and introduced in [44,64], see also [31].
To calculate the symmetry energy it is necessary to know the energy per nucleon at various densities and temperatures, and for various values of the isospin asymmetry, x. Figure 19 shows the behavior of the energy per nucleon as a function of x for temperatures T = 0.5, 1.0, and 2.0 MeV, and densities ρ = 0.04, 0.06, 0.08, and 0.10 fm −3 . Following the procedure of [31,44,64] (cf. Appendix B), one can compute the symmetry energy both of nuclear matter and of neutron star matter. Figure 20 shows the variation of E sym with the temperature for four sub-saturation densities.
The behavior of E sym for neutron star matter is reminiscent to that obtained with molecular dynamics with the "new medium" potential for nuclear matter. Relatively constant values of E sym are observed at temperatures below 2.0 MeV, corresponding to the pasta regime, and with smaller values at higher temperatures corresponding to more uniform systems. See [31,38] for a related analysis.
The observed values of E sym , namely, ranging between 20 MeV and 50 MeV, are reminiscent of previous results obtained with microscopic field theories, see e.g., [23]. Similarly, the trend of E sym , decreasing for larger temperatures, is similar to the behavior found using an equation of state obtained through a virial expansion at T = 2, 4, and 8 MeV [65], and through a self-consistent model using various effective interactions [66]. Comparing these results with previous computations of the symmetry energy obtained with classical potentials, e.g., those of Figure 68b of [31], we can see a major change in the behavior of E sym at low temperatures; surely due to the effect of the repulsion introduced at low momentum transfer.

Discussion
In this article we presented a nucleon-nucleon potential whose strength depends on the relative momenta of the colliding nucleons and, thus, has the ability to reduce interactions at low momentum transfer mimicking the Pauli exclusion principle.
Comparing to previous attempts at introducing Pauli blocking to a classical potential, we improved on the best attempt (i.e., that of Dorso and Randrup from 1987 [49]) by constructing a potential that respects an excluded volume in phase space, respects the uncertainty relation, and inhibits the production of di-neutrons.
As a first test, the potential was used with Metropolis Monte Carlo simulations to study the formation of finite nuclei, their binding energies, radii, and internal nucleon distribution. It was also used to study the binding energy in infinite systems such as nuclear matter and neutron star matter, as well as the structures formed at sub-saturation densities and cold temperatures, known as nuclear pastas. Finally, the symmetry energy was also obtained for the case of neutron star matter.
In general terms, the potential functioned as expected in most areas reproducing binding energy and radii trends of finite nuclei. The saturation density and compressibility of infinite systems were observed near the expected values. Pasta structures were observed to appear at low temperatures and sub-saturation densities. Finally, the magnitude and trends of the symmetry energy were in the ball park predicted by mean field theories.
An observed weakness of the potential is the over-binding found in the synthetic nuclei (cf. Figure 3). By inspecting the radii of the nuclei in Figure 2, and the distance between nucleons in Figure 4, it is easy to see that the model produces nuclei that are too tight, with nucleons closer to one another more than expected, which leads to over-binding. Since the over-binding and smaller radii are a problem of finite nuclei, we do not expect that they will affect the results of infinite systems, vis-à-vis the pasta structures and symmetry energy as they are obtained from systems with fixed densities.
Ongoing investigations are addressing these problems, especially as future studies will concentrate on using finite nuclei in collisions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Metropolis Monte Carlo
The simulations were performed by means of the Metropolis Monte Carlo procedure. It consists essentially on building an appropriate Markov chain by sampling the configurational space. The system was assumed to be in thermal equilibrium with an external heat bath, and thus, the sampled energy states were those of the canonical ensemble.
The nucleons were first placed at the vertices of a cubic lattice, in order to avoid overlapping. The momenta were further set to random with a vanishing mean value. Then, the system was driven to the thermal equilibrium following the same procedure as in the simulation stage (see below), but fixing the thermal bath at T = 5 MeV. The thermalization required 1000 steps, that is, 1000 times the degrees of freedom of the system to arrive to equilibrium. The equilibrium temperature, T = 5 MeV, was warm enough to mimic a liquid state. This can be verified through the radial pair distribution and the distribution of the components of the momentum; Figures A1 and A2 show the radial distribution function and the momentum distribution, respectively, for the case of ρ = 0.04 fm −3 with 1000 nucleons. All the simulations ran over 1000 nucleons for nuclear matter and 6000 for neutron star matter.
The simulation stage was as follows: At each step, any nucleon was selected randomly, and either moved to a new position, or changed its momentum. This yields an energy jump ∆E. Such an energy jump corresponds to the following transition probability in the canonical ensemble. where β = 1/KT (K meaning Boltzmann constant and T is the associated bath temperature). Thus, accepting the transition with probability p attains the right samples for the canonical space. The above "trials" alternate between the positions and momenta of randomly selected particles. The fraction of accepted trials commonly lies in the interval 0.35-0.75. A control mechanism was implemented to keep this fraction between these limits. The controlling magnitudes were the degree of the perturbation of the positions and momenta. The final Markov chain was afterward sub-sampled, in order to get non-correlated observables.