Properties of nuclear pastas

i

This article has been updated. A Correction to this article was published on .

i

The Original Article was published on

i

The Original Article was published on

i

This article has been updated

i

This article was retracted on

i

This article has been updated

Abstract

In this review we study the nuclear pastas as they are expected to be formed in neutron star crusts. We start with a study of the pastas formed in nuclear matter (composed of protons and neutrons), we follow with the role of the electron gas on the formation of pastas, and we then investigate the pastas in neutron star matter (nuclear matter embedded in an electron gas).

Nuclear matter (NM) at intermediate temperatures (1 MeV \lesssim T\lesssim 15\text{MeV}), at saturation and sub-saturation densities, and with proton content ranging from 30% to 50% was found to have liquid, gaseous and liquid–gas mixed phases. The isospin-dependent phase diagram was obtained along with the critical points, and the symmetry energy was calculated and compared to experimental data and other theories. At low temperatures (T\lesssim 1\text{MeV}) NM produces crystal-like structures around saturation densities, and pasta-like structures at sub-saturation densities. Properties of the pasta structures were studied with cluster-recognition algorithms, caloric curve, the radial distribution function, the Lindemann coefficient, Kolmogorov statistics, Minkowski functionals; the symmetry energy of the pasta showed a connection with its morphology.

Neutron star matter (NSM) is nuclear matter embedded in an electron gas. The electron gas is included in the calculation by the inclusion of an screened Coulomb potential. To connect the NM pastas with those in neutron star matter (NSM), the role the strength and screening length of the Coulomb interaction have on the formation of the pastas in NM was investigated. Pasta was found to exist even without the presence of the electron gas, but the effect of the Coulomb interaction is to form more defined pasta structures, among other effects. Likewise, it was determined that there is a minimal screening length for the developed structures to be independent of the cell size.

Neutron star matter was found to have similar phases as NM, phase transitions, symmetry energy, structure function and thermal conductivity. Like in NM, pasta forms at around T ≈ 1.5 MeV, and liquid-to-solid phase changes were detected at T ≈ 0.5 MeV. The structure function and the symmetry energy were also found to depend on the pasta structures.

Graphical abstract

Keywords

nuclear pasta / neutron star matter / nuclear symmetry energy / molecular dynamics / nuclear phase transitions

1 Introduction

Neutron stars are the final product of the death of a massive star. Stars die when the internal thermonuclear fusion can no longer balance gravitational compression, and a supernova shock ejects most of the mass of the star leaving behind a dense core. If the star is big enough as to generate a remnant core with a mass larger than the Chandrasekhar limit, the gravity will be so strong that the electron degeneracy pressure will not be able to support the system, and the star will compress down to nuclear density. Since during the collapse a big part of electrons and protons turn into neutrons and escaping neutrinos through electron capture, the produced core tends to have an excess of neutrons over protons, thus justifying the name of a neutron star. A neutron star is a neutral object composed of neutron, protons and electrons, with the charge of the protons balanced by the one of electrons.

Neutron stars range in mass between 1 and 3 solar masses, with a radius is of the order of about 10 km, and the average density of a degenerate Fermi gas composed of nucleons, 10¹⁵ g/cm³ or about 20 times that of normal nuclei. Energy considerations indicate [1, 2] that the cores of neutron stars have a 1-km crust where the β decayed-produced neutrons form neutron-rich nuclear matter immersed in a sea of electrons. The density of neutron star crusts goes from normal nuclear density (≈3×10¹⁴ g/cm³, or ≈ρ₀=0.16 nucleons/fm³) at a depth of ≈1 km to the neutron drip density (≈4×10¹¹ g/cm³) at ≈1/2 km, to, finally, an even lighter mix of neutron-rich nuclei also embedded in a sea of electrons with densities decreasing down to practically zero in the neutron star envelope.

It is known that nuclear systems at low temperatures and sub-saturation densities can exist in non-uniform structures, known as “pastas”, and it is suspected that the crusts of neutron stars may conform to such configurations [1–9]. Knowing the structures attained by neutron crusts is an important factor in the understanding of star-quakes, pulsar frequencies and neutron star evolution. Indeed since the cooling of non-pulsar neutron stars is due mostly to neutrino emission from the core, the interaction between neutrinos and the crust structure is relevant to their thermal evolution [10]. The study of such structures is the main goal of this review.

1.1 The pasta

The pasta structures of nuclear systems have been studied using different models [1–9]. Such structures appear to form due to the interplay between attractive-repulsive nuclear and Coulomb forces. As the density, temperature and proton fraction vary, the structure changes from a uniform phase at the core to configurations with voids filled with small clusters “gnocchi-like”, to “lasagna-like” layers of nuclear matter and gas, to “spaghetti-like” rods of matter embedded in a nuclear gas, to “gnocchi” clumps, to a practically uniformly dissolved gaseous phase [2].

To attain these structures, nucleons should reach states of a free energy minimum. In the original studies of the 1980s the energy minima were determined using static methods such as the liquid drop model [1, 9], mean field theories [11] and Thomas–Fermi models [8]. These methods, however, usually work only at zero temperature, locate only the “traditional” global minima and tend to miss the “non-traditional” local minima of energy barriers; to gain access to the complete set of structures (traditional and non-traditional) temperature-dependent dynamical models are needed.

From a different point of view, an approximate scheme known as Energy Density Functional (EDF) has also been used to study the pasta. EDF models the energy of the system either by the use of effective nuclear forces (adjusted to fit the known binding energies of nuclei), or by the use of a function of one-body density and its gradients. Recent calculations [12, 13] apply EDF to predict nuclear masses, radii or deformations along the nuclear chart, as well as the elastic and inelastic parts of the structure factor for uniform matter. A problem for the application of EDF to nuclear matter and neutron star matter lies on extending the functionals into untested regions of non-uniform matter; in the mentioned calculations EDF was coupled to an analytical model for the minimal surfaces (which divides matter from vacuum) for the basic pasta structures, managing to produce an elastic structure factor with rich and distinctive patterns, which can have an effect on neutrino scattering.

Dynamical studies, such as quantum molecular dynamics [3, 6, 14–17], predict the formation of the pasta phases by dynamical means, but tend to be limited to states of global energy minima. To properly obtain the pasta structures, it is best to use models that can achieve phase changes, clusterization and identify local minima by cooling, such as classical molecular dynamics models [18–21].

In this review we will use the classical molecular dynamics model (CMD) to study the pasta in nuclear matter and in neutron star matter; systems which will be defined in the next section. To avoid technical distractions, the description of the CMD model is relegated to Appendix A.

1.2 The pasta in nuclear matter and in neutron star matter

The neutron star crust is composed of protons, neutrons and an embedding gas of electrons. To study such system we systematically divide the review in three parts. First, a study of nuclear matter (NM), e.g., systems composed solely of protons and neutrons; second, an investigation of the role electrons have on the properties of NM; and third, a study of neutron star matter (NSM), i.e., systems of protons and neutrons embedded in an electron gas.

Nuclear matter exhibits fascinating complex phenomena at subsaturation densities and warm and cold temperatures. At densities below the saturation density, ρ₀=0.16 fm⁻³, and temperatures, say, between 1 MeV and 5 MeV, nuclear systems exist in liquid and gaseous phases, as well as in a mixture between the two. At lower temperatures, crystals and structures resembling the so-called nuclear pastas appear with different morphologies depending on the temperature, density and isospin content and, furthermore, nucleons inside such structures can undergo phase transitions.

This review studies nuclear matter first. At intermediate temperatures we review the bulk properties, phases and phase transitions of NM, all with varying percentages of isospin content to allow the study of the symmetry energy. At lower energies, the formation of the pasta is investigated along with its properties, phases, transitions, and symmetry energy.

This is followed by a review of the role the electron gas has on the pasta of nuclear matter. As the electron gas is introduced in the CMD calculation by means of a screened Coulomb potential, the study focuses on the effects the strength and range of the interaction have on the morphology of the pastas.

Finally, the properties of the pastas in neutron star matter, i.e., systems with protons, neutrons and electrons, are reviewed with special attention to their shapes, phases, phase transitions, and symmetry energy.

In understanding the pastas both in NM and in NSM we also study how the energy varies with respect to the isospin content, i.e., the symmetry energy. It is in this type of studies that one can appreciate the value of CMD, as to perform such a study at intermediate and low temperatures and sub-saturation densities, a model capable of exhibiting clustering phenomena and phase changes is needed.

Similarly, as mentioned before, the cooling of neutron stars is related to pasta structure of the crust as the neutrino opacity is enhanced by coherent scattering. Here, we also study the very-long range order of the pasta phases, focusing on its influence on the opacity of the crust to low-momentum neutrinos.

In this review the properties of nuclear matter will be studied in Section 2, the role of the electron gas on the NM will be investigated in Section 3 and, finally, the neutron star matter will be studied in Section 4. Analytical tools used, such as the classical molecular dynamics method used, algorithms for the recognition of clusters, the radial correlation functions, Lindemann coefficient, Kolmogorov statistics, Minkowski functionals, and the procedure used to calculate the symmetry energy are presented in the various appendices.

2 Nuclear matter

In this section we study properties of nuclear matter. NM is composed of protons and neutrons in varying proportions. The properties of NM vary depending on the temperature and density of the system, as well as on its isospin content. At intermediate temperatures, that is between 1 MeV \lesssim T\lesssim T_{critical}\approx 15\text{MeV}, NM can exist in liquid form, gaseous form, and in a mixture between these phases. At lower temperatures crystal structures develop around normal saturation nuclear density, ρ₀=0.16 fm⁻³, and non-homogeneous structures (pasta-like) at sub-saturation densities. We will study NM at intermediate temperatures in Section 2.1, and at lower temperatures in Section 2.2.

For the purposes of the present investigation, nuclear matter is taken as protons and neutrons interacting only through the Pandharipande potentials, i.e., the proton–proton interaction is ignored (since the Pandharipande potentials yield the correct scattering cross sections, it can be though of as containing the proton–proton Coulomb interaction effectively). Later, in Section 3, NM will be embedded in an electron gas and the proton-electron interaction will be taken explicitly through an screened Coulomb Potential; this will effectively constitute neutron star matter as studied in Section 4.

2.1 Nuclear matter at intermediate temperatures

Here we study systems with different values of isospin content, density and temperature to obtain bulk properties such as the energy per nucleon, pressure, equilibrium density, compressibility, symmetry energy, as well as phases and phase transitions. An overview of the bulk properties of nuclear matter at intermediate temperatures is presented in Section 2.1.1, followed by a discussion of the existence of phases in Section 2.1.2, the phase diagram in Section 2.1.3, and an estimation of the nuclear symmetry energy in Section 2.1.4; these findings are summarized in Section 2.3.

2.1.1 Bulk properties of nuclear matter at intermediate temperatures

The studied systems were created with the molecular dynamics code described in Appendix A, implemented according to the details described in Appendix A.2.1. In summary, a total of 2000 nucleons were placed in cubic cells of the appropriate size to have densities ρ = 0.02, 0.03, …, 0.2 fm⁻³, and with periodic boundary conditions to simulate an infinite system. The ratio of protons to neutrons was fixed to values x = Z/A = 0.3, 0.4, and 0.5, where Z is the number of protons, and A is the total number of nucleons. The temperatures of the systems studied are T = 1, 2, 3, 4, and 5 MeV.

The energy per nucleon (kinetic plus potential) is shown in Fig. 1 as a function of the density at various temperatures. The curves show characteristic “∪” shapes around their saturation density (minimum of the ∪); these shapes resemble the predictions of Skyrme–Hartree–Fock and relativistic mean-field calculations for zero temperature [22]. At lower temperatures all three cases presented in Fig. 1 have regions of negative binding energy which correspond to bound (liquid) matter; at high temperatures the systems become unbound. At lower densities the curves depart from the ∪ shapes showing a phase transitions from a uniform medium around saturation densities to non-homogeneous media at sub-saturation densities.

Figure 2(a) shows the pressure–density curves for T = 1 and 5 MeV, and x = 0.3, 0.4, and 0.5 and, corroborating the previous findings, the equilibrium densities (zero pressure points) correspond to the minima of the energy–density curves. It must be remarked that at low densities the CMD pressure is not correct as the system should approximate a free nucleon gas, which is a quantum system (see Appendix B). The corresponding saturation densities can be estimated through a fit of the minima of E(T\mathrm{,}\rho ); Fig. 2(b) shows the saturation density as a function of the temperature. At near zero temperatures ρ₀ varies from 0.16 fm⁻³ at x = 0.5, to a low 0.04 fm⁻³ at x = 0.3. The compressibility at saturation density, obtained from the energy through K(T\mathrm{,}\rho )=9{\rho }^2{\left({\partial }^{2}E/\partial {\rho }^2\right)}_{{\rho }_0}, is presented in Fig. 2(c). Nuclear matter becomes a very soft fluid for x = 0.3 at all temperatures, and it varies drastically with x, and reduces by ≈30%, as T increases from 1 to 5 MeV for x = 0.4 and 0.5. This softening with excitation energy is consistent with IQMD simulations of ¹⁹⁷Au+¹⁹⁷Au at 600 MeV/A [23] and with BUU calculations of similar reactions at 1 GeV/A [24].

2.1.2 Phases of nuclear matter at intermediate temperatures

The existence of phases in nuclear matter can be seen from Fig. 1. In the T = 1 MeV curve the transition from a smooth ∪ shape to an extraneous curve happens at ρ ≈ 0.10 fm⁻³ for x = 0.5, and at ρ ≈ 0.08 fm⁻³ for x = 0.4, but it is not as pronounced for the x = 0.3 case. The smooth ∪ shape corresponds to a uniform crystal-like phase at low temperatures and liquid-like phase at higher temperatures, and the lower-density part corresponds to non-homogeneous crystal–gas mixture at low temperatures or liquid–gas mixture at higher temperatures. As we will see in Section 2.1.3, systems with x<0.3 would appear to be unable to reach such a liquid–gas mixture region and would always stay in a liquid-like continuous medium down to very low densities.

To further investigate these phases, we examine the T = 1 MeV systems through the radial distribution function and their mass distribution. As explained in Appendix G.1, and more specifically in Appendix G.1.1, the radial distribution function g(r) signals the presence of neighbors and their average distance, which is different in liquid and gaseous states. Figure 3 shows g(r) for systems at T = 1 MeV at a liquid density (ρ = 0.16 fm⁻³, bottom panel) and at a liquid–gas mixture density (ρ = 0.06 fm⁻³, top panel). The strengths of the nearest-neighbor peaks show that at low densities nucleons tend to be more correlated than at higher densities indicating that at ρ = 0.06 fm⁻³ the main contribution at short distances is from nucleons in droplets, while at ρ = 0.16 fm⁻³ the larger nucleon mobility reduces such correlation. This is also observed in the second-neighbor peaks which appear at the same distance for all values of x at ρ = 0.06 fm⁻³ but not at ρ = 0.16 fm⁻³ indicating again a reduced mobility of nucleons in the droplets. The growth of the second-neighbors peak over the first peak for the case of x = 0.3, is due to the large number of nn repulsive interactions which exceed the smaller number of np attractive interactions at such large discrepancy between the number of neutrons (70%) and protons (30%).

Further information about the phases found before can be obtained by the mass distributions attained by the systems in the liquid phase and in the liquid–gas mixture; the method used to recognize cluster is outlined in Appendix G.6. Looking at the x = 0.3 cases, Fig. 4 shows the mass distribution in the central cell (i.e., not continuing the fragments into neighboring cells) at ρ = 0.06 fm⁻³ and T = 1 MeV in the liquid–gas mixture at ρ = 0.16 fm⁻³ compared to that in the liquid–gas mixture at ρ = 0.06 fm⁻³. Figure 4 shows that in the mixture zone (ρ = 0.06 fm⁻³) there are more intermediate-mass droplets than in the liquid phase (at ρ = 0.16 fm⁻³). The relatively large number of clusters with A>1000 in the liquid phase, indicates that the system remains in a continuous medium; notice that such large cluster necessarily must have low populations for obvious reasons.

The results obtained so far indicate that liquid and gaseous phases exists in nuclear matter with isospin asymmetries in the range from x = 0.3 to 0.5. Transitions between these phases appear to be are alive and well in these cases of high isospin asymmetry. A complementary observation is that these systems appear become unbound at temperatures as low as T = 2 MeV for x = 0.3. It is known that pure neutron matter is unbound, i.e., it exists only in gaseous phase [25]; in the next section we determine the boundaries of these phases.

2.1.3 Phase diagram

The phase diagram of nuclear matter has been studied for isospin symmetric matter in Refs. [26–28], and evaluated for isosospin asymmetric matter in Ref. [29], but it was obtained dynamically with full inclusion of the isospin degree of freedom the first time in Ref. [30]. The phase diagram of nuclear matter can be obtained from the pressure–density isotherms of Section 2.1, cf. Fig. 2(a), by means of the technique known as the Maxwell construction, described in Appendix B.

Pressure-density isotherms, obtained by CMD simulations of systems with isospin content of x = 0.3, 0.35, 0.4, 0.45, and 0.5, at temperatures varying between T = 1,2,…,15 MeV, and densities ρ = 0.01,0.02,…,0.2 fm⁻³, were used in the Maxwell construction. Figure 5 shows the T = 11 MeV isotherms for x = 0.3, 0.4, and 0.5.

As explained in Appendix B, the application of the Maxwell construction to each pressure isotherm yields the value of the density corresponding to the boundary between the liquid–gas mixture and the gaseous phase, ρ_Gas and that of the limit of the mixture and the liquid phase, ρ_Liq for such isotherm. When plotted together in a density–temperature plane, all such points outline the liquid–gas coexistence region; Fig. 6 shows the points obtained from the Maxwell construction for the cases x = 0.5 and x = 0.3. The area delineated by such points indicates the liquid–gas coexistence region; the area to the left and to the right of the coexistence region correspond, respectively, to gaseous and liquid phases; the continuous and dashed lines are simple approximations to indicate the boundaries of the liquid–gas mixture regions for the two values of x. The maximum of each curve are the respective critical points for each value of x.

Also shown in Fig. 6 are six characteristic configuration of protons and neutrons for two cases on the boundary of the gaseous phase: T = 1 MeV and ρ = 0.02 fm⁻³ for x = 0.3 and 0.5 (left two frames), two cases in the liquid–gas mixed regions: T = 5 MeV and ρ = 0.04 fm⁻³ for x = 0.3 and 0.5 (central top and bottom frames), and two more cases on the boundary of the liquid phase: T = 1 MeV and ρ = 0.06 fm⁻³ for x = 0.3 and 0.5 (right two frames). Each configuration shows the protons with dark spheres and neutrons with white ones. For clarity the boxes are shown approximately of the same size, but in reality they are drawn at different scales. The varying sizes of the spheres are due to the perspective effect of the Visual Molecular Dynamics (VMD) software used, and are not intended to represent clusters.

The right part of each of the six configurations shown in Fig. 6 is an attempt to quantify the connectivity of the phases. The surfaces shown have similar particle densities obtained from a volumetric Gaussian density map computed from particles in a grid. In all of these cases the isosurfaces were calculated with VMD’s QuickSurf [31] application with radius scale of 1.0, density isovalue of 2.5 and grid spacing of 0.5, and clearly show an evolution from low connectivity in the gaseous phase, to mid-connectivity in the liquid–gas mixed phase, to high connectivity in the liquid phase.

Combining the density–temperature points from the Maxwell constructions for the cases x = 0.3, 0.35, 0.4, 0.45, and 0.5, it is possible to construct a three-dimensional plot of the liquid–gas mixture region in the density–x–temperature space. Figure 7 shows the surface of the liquid–gas coexistence region; notice that and the critical points now join into a critical ridge. The surface was obtained with the instruction ListPlot3D of the Mathematica package.

Figure 7 indicates that the liquid–gas mixture region ends at a value of x smaller than 3, but this is due to a lack of data for lower values of x. A more insightful estimate of the lower x boundary of the coexistence region can be obtained by looking at the trend of the critical temperatures.

Figure 8 indicates how the critical temperatures decline for lower values of x. The points in the figure show the critical temperatures obtained from the critical ridge, i.e., from the maxima of the fits used in Fig. 6 for all values of x. The continuous line is a parabolic fit, which has been extended to lower values of x. An artistic touch was added with an Escher-like background indicating the region where bubbles and droplets coexist. As the curve reaches the T = 0 value at x = 0.1324, we conclude that such point is the lower x limit of the liquid–gas mix region, and that outside of this ρ–T–x volume all nuclear systems are expected to be fully unbound. Although this coexistence region has never been calculated before, our findings are in general agreement with the low binding energy of x = 0.3 matter found in Ref. [22], and of unbound pure neutron matter (x = 0) found in Ref. [25] at zero temperature.

2.1.4 Symmetry energy of nuclear matter at intermediate temperatures

The symmetry energy can be obtained from the CMD results of E(T,ρ,x) using a method first introduced in Ref. [32] for nuclear matter around liquid densities, and in Ref. [33] for nuclear matter in the liquid–gas mixed phase. For completeness the method for calculating the symmetry energy of nuclear matter at intermediate temperatures is presented in Appendix C.

Figure 9(a) shows E_{Sym} at low densities (in the liquid–gas coexistence phase). The figure compares the symmetry energy to the experimental points of [34–36]. Notice that the CMD-based symmetry energies approach non-zero values around the binding energy of alpha clusters in the limit of low density, as demanded by Natowitz and coworkers [34]. The reason for this is that nuclear matter, just like CMD, form alpha-like clusters at low densities, and thus the energy tends to this value.

Similarly, Fig. 9(b) presents E_{Sym} expanded up to liquid densities to compare with a relativistic Hartree calculation (dashed line labeled “NL₂”) and other field theories calculations at their corresponding saturation densities; see [37] for complete details.

It should be pointed out that, as discussed in Ref. [32], the numerical procedure used to obtain E_{Sym} (cf. Appendix C) has enough flexibility to make the CMD calculation fully agree with the calculations shown in Fig. 9(b). Likewise, some characteristics, such as the temperature dependence of E_{\text{sym}}(T\mathrm{,}\rho ) are in agreement with the variations of liquid-drop terms and nuclear surface tension in the range of T\lesssim 5\text{MeV} [38] but in disagreement with trends obtained by other theories [39, 40]. Although the CMD symmetry energy is reminiscent of relativistic and non-relativistic Hartree calculations [37], there are notable differences mostly due to the distinct assumptions of the models, see [32] for a more complete discussion.

2.2 Nuclear matter at low temperatures

At low temperatures, T<1 MeV, nuclear matter is known to produce crystal-like structures around saturation densities, and pasta-like structures at sub-saturation densities. In this section we investigate such structures and their properties, namely their morphology, phase changes and their phase transitions. These studies require specialized tools such as the caloric curve, Lindemann coefficient, radial distribution function, Kolmogorov statistic, and Minkowski functionals. The calculation of the symmetry energy also needs a special procedure based on the one used at intermediate energies (cf. Appendix C).

2.2.1 Nuclear matter at low temperatures around saturation density

To study the uniform phase at zero temperature, a simple-cubic crystalline structure at a given density was constructed and the energy per nucleon calculated by direct summation between all nucleons. The dependence of the binding energy of the uniform phase on the density is explored by scaling the lattice parameter. As shown in Fig. 10, this yields the characteristic ∪-shaped curve, with a minimum at the saturation nuclear density.

For comparison, and given that nuclear matter is composed of neutrons and protons, systems with crystal geometries similar to those formed by binary alloys were constructed. Figure 10 shows a simple cubic lattice in which every first neighbor of a proton is a neutron and vice versa (SC), a BCC lattice, and a diamond lattice with nucleons arranged so that every first neighbor of a protons is a neutron and viceversa (FCC). It is clear that the simple cubic is the one with an energy closer to the nuclear binding energy. See [41] for more details.

Since the previous scaling cannot yield non-homogeneous systems, i.e., a phase transition, CMD was used with systems at densities around and below saturation. Starting from a random positioning of nucleons in a cubic cell under periodic boundary conditions, nucleons were endowed with velocities sampled from a Maxwell–Boltzmann velocity distribution corresponding to a given initial temperature. The systems then were evolved until equilibrium is achieved at a high temperature (T\gtrsim 2\text{MeV}) and then brought down to the final desired temperature of T = 0.001 MeV using the Andersen thermostat in small temperature steps. After reaching equilibrium, the analysis tools described in Appendix G were used to characterize and visualize the produced configurations.

Figure 10 shows the near zero temperature results for symmetric (x = 0.5) matter. As it can be seen, CMD reproduces the simple cubic (B1) lattice calculations, corresponding to the uniform phase, up to a density of ρ ≈ 0.13 fm⁻³, while for lower densities the systems goes into a series of non-homogeneous phases. Figure 11 shows visual representations of the structures corresponding to the four densities 0.01 (a), 0.06 (b), 0.13 (c), and 0.16 fm⁻³ (d).

The compressibility can be obtained through K=9{\rho }_0^2{({\text{d}}^{2}E/\text{d}{\rho }^2)}_{{\rho }_0}, where E(ρ,T,x) is a simple degree-two polynomial fit of the B1 curve of Fig. 10 around the saturation density. Figure 12(a) shows the energy, pressure and compressibility for the homogeneous B1 lattice. The bulk modulus for the medium Pandharipande potential is found to be 283 MeV, comparable to the value of 250 MeV quoted by its creators [42], who used different fitting polynomial. The pasta-like structures are found in a mechanically unstable (negative pressure) density region, but well above the divergence in compressibility that would signal the thermodynamic instability.

At slightly higher temperatures, between 0.001 MeV≤T≤1.0 MeV, the NM systems follow the ∪ shape characteristic of the uniform T = 0 crystalline phase at densities \rho \gtrsim 0.13 fm⁻³. Figure 12(b) shows the same type of results as Fig. 10, again, at a density ρ ≈ 0.13 fm⁻³ the systems move away from the uniform phase forming non-homogeneous stable arrangements. Figure 13 shows some of those structures at T = 1.0 MeV for the four densities labeled from (a) to (d) in Fig. 12(b).

These non-homogeneous structures can be characterized using the mean curvature and Euler characteristic. As explained in detail in Appendix G, the Euler characteristic is a topological invariant that describes a structure regardless of the way it is bent. Different structures have distinct values of these variables and, in general, follow the pattern outlined in Table 6. As a reference, the perfect crystals formed at T = 0 and ρ ≈ 0.13 fm⁻³ [v.g. point (A) in Fig. 10 and Fig. 12(b)] are formally uniform and infinite because of the periodic boundary conditions imposed, hence they have no surfaces and null Euler characteristic.

Figure 14 shows the effect of temperature on the curvature–Euler coordinates of the labeled structures of Fig. 12(b). Points on the lines have the same densities but their temperatures vary from T = 0.001, 0.1, 0.5, 0.6, 0.8, and 1.0 MeV. The labels corresponding to those of Fig. 12(b): cases “A-a” at normal density (ρ₀) correspond to uniform crystalline structures so all have zero curvature and Euler number at the temperatures studied; cases “B-b” of density 0.1 fm⁻³ the almost spherical bubbles at T = 0.001 MeV become distorted at higher temperatures; cases is “C-c” at ρ = 0.06 fm⁻³, go from being a perfect “lasagna” to a “jungle-gym”, a complex of lasagna joined by columns (and similar to the Schwarz’s Primitive surface of differential geometry); and cases “D-d”, on the other side, go from spherical “gnocchi” to deformed droplets.

Complementary studies for higher compressibility potentials, the stiff Pandharipande potential and the effect of the periodic boundary conditions can be found in reference [41].

2.2.2 Properties of NM pasta-like structures

A big advantage of using CMD to study the pasta-like structures, is that the nucleon dynamics can drive the system into states with local free energy minima, which abound in complex energy landscapes, and find nontraditional (amorphous, sponge-like) structures. Such states cannot be found by static models due to the energy barriers which usually surround local minima [19].

To characterize the pasta, its phase changes and to calculate the symmetry energy we use the cluster-recognition algorithms, the radial distribution function, the Lindemann coefficient, Kolmogorov statistics, Minkowski functionals and a numerical method to estimate E_{\text{sym}}. These techniques are reviewed in Appendices D and G.

a. Phases

The temperature–excitation energy correlation known as the caloric curve, has been used in nuclear physics to detect first-order liquid–gas phase transitions [43]. The caloric curve simply relates the temperature of the system to the energy input needed to reach such temperature. Whenever energy is used for something other than heating up the system, like breaking bonds while melting a crystal, energy will flow in but the temperature will not increase, such discontinuities in the slope of the caloric curve will signal phase changes. In this section the caloric curve is computed for systems with three different proton ratios to identify phase changes. See Ref. [19] for complete details.

Figure 15(a) shows the caloric curve at ρ = 0.05 fm⁻³ for systems of 6000 nucleons at proton ratios of x = 0.3, 0.4 and 0.5. Two changes of slope are noticeable, a sharp one close to T = 0.5 MeV, and one less conspicuous just before T = 2.0 MeV; the changes are more pronounced in symmetric matter. Similar changes are seen for other cases on the left panels of Fig. 16.

The significance of the slope changes of the caloric curve can be investigated with the Lindemann coefficient. As explained in Appendix G.2, the Lindemann coefficient provides an average estimation of displacement of the nucleons, and it can signal a change of mobility and, hence, a phase transition. Figure 15(b) shows the behavior of the Lindemann coefficient as a function of the temperature superimposed on the caloric curve of Fig. 15(a). It is easy to see that at lower temperatures the nucleon mobility is only a fraction of its higher-temperature value; the change at T = 0.5 MeV appears to correspond to a solid-liquid phase transition within the pasta regime, as concluded in Ref. [19].

The pressure-temperature curves can be used to construct phase diagrams. The right panels of Fig. 16 show the behavior of the system pressure as a function of the temperature. Comparing to the left panels, one can see that the p–T curves experience slope changes at similar temperatures as the E–T curves. Notice that the pressure of the symmetric nuclear matter (x = 0.5) also changes sign, while the pressure of the asymmetric nuclear matter (x<0.5) remains positive until very low temperatures. It is clear that the former enters into the metastable regime, while the latter does not.

The microscopic structure of the different phases can be explored through the radial distribution function g(r), introduced in Appendix G.1.2. To study the caloric curve results, g(r) was calculated in the three temperature regions found, namely, in T<0.5 MeV, 0.5 MeV<T<2.0 MeV, and T>2.0 MeV. Figure 17 shows g(r) for nuclear matter at ρ = 0.085 fm⁻³ and at temperatures in the three different regions. Panel (a) show the results for the symmetric (x = 0.5) case, and (b) the non-symmetric (x = 0.4) case.

The peaks of g(r) in Fig. 17 indicate the average position of the nearest neighbors. At T = 0.2 MeV the peaks of panel (a) show large correlations at 1.85 fm, 2.65 fm, and 3.25 fm; such uniformity of nearest- and second-neighbors is characteristic of a solid state. Furthermore, since the peaks occur at a\mathrm{,}\sqrt{2}a\mathrm{,}\sqrt{3}a (where a is the nearest-neighbor distance), it indicates that the solid is a simple cubic arrangement. The gradual loss of correlation at T = 1 MeV and 2 MeV indicate that the phase is losing correlations, as expected in melting. Similar behaviors are observed in panel (b) for the case of x = 0.4. Notice that both the T = 1 MeV and 2 MeV cases show a rapid drop of correlations at distances larger that, say, 4 fm; this is expected as such distances lie outside the pasta structure. The results from the caloric curves, g(r), Lindemann coefficient, and the p–T curves, together indicate that the nucleons inside the pasta structures exist in different phases.

The argument outlined so far does not hypothesize on the pasta shape. It only compares the average density inside the pasta with respect to N/V. Thus, it is expected to hold on a variety of pastas. This was verified for different densities, as shown in Fig. 18. Figures 19 and 20 show the pasta structures that correspond to Fig. 17(b). Panel (a) shows all nucleons (protons in light color and neutrons in dark), (b) protons only, and (c) neutrons only. Among other things, this indicates that pasta structures exists in non-symmetric nuclear matter.

b. The onset of the pasta

One way of determining the onset of the creation of the pasta is by looking at the spatial distribution of nucleons during the cooling of a nuclear matter system, and by quantifying the geometrical shape of the pasta structure. For these tasks we use the Kolmogorov Statistics and the Minkowski Functionals.

The spatial characteristics can be quantified by means of the Kolmogorov statistics, which measures the difference between a given distribution of nuclei and a homogenous one. As the temperature decreases and the pastas are created, the spatial distribution of nucleons will depart from a uniform distribution. The Kolmogorov statistic is presented in Appendix G.3. In the present case we treat the distribution in separate coordinates (more useful for the case of planar structures, such as “lasagnas”), case for which the Kolmogorov statistic is “distribution-free” or parameter-free. Figure 21 shows the corresponding results for an isospin symmetric system at ρ = 0.05 fm⁻³; these results are from a single simulation, since the Minkowski statistics can only be applied to raw data, not to averaged one; see caption for further details on the simulation procedure.

As can be seen in Fig. 21, the values of the 1D Kolmogorov statistic, i.e., the discrepancies to homogeneous distributions, are negligible for temperatures above 2 MeV, as expected for energetic particles moving around homogeneously. At T ≈ 2 MeV, the three statistics experience a change in the slope, although D_x and D_z return to negligible values as the temperature further decreases. The D_y statistic, instead, attains a definite departure from homogeneity for T<2 MeV. A lasagna or slab-like structure across the y-axis occurs at conditions ρ = 0.05 fm⁻³, T = 0.2 MeV and x = 0.5, as presented in Fig. 82(a).

The 1D Kolmogorov statistic attains the departure from homogeneity at an early stage of the pasta formation. The arrow in Fig. 21(b) points to the most noticeable bubble appearing in the system at T = 2 MeV. The bubble-like heterogeneity also explains the changes in the slope for D_x and D_z at this temperature, as pictured in Fig. 22. For decreasing temperatures, the bubble widens (on the right side of the image due to the periodic boundary conditions), a tunnel appears [Fig. 22(b)], and at T<1.5 MeV it finally splits into two pieces while the (x,y) homogeneity gets restored [Fig. 22(c)] returning D_x and D_z back to their negligible values.

We further applied the 1D Kolmogorov statistic on protons and neutrons separately for non-symmetric nuclear matter systems. Figure 23 exhibits the most significant 1D statistics for the densities ρ = 0.05 fm⁻³ and ρ = 0.85 fm⁻³, respectively. The corresponding spatial configurations can be seen in Figs. 19 and 20 for ρ = 0.085 fm⁻³, and in Figs. 24 and 25 for ρ = 0.05 fm⁻³.

According to Fig. 23, the temperature threshold at which the 1D Kolmogorov statistic becomes significant, decreases with smaller proton ratios (at fixed density), in agreement with the caloric curves of Section 2.2.2a, but more sharply exposed now. Both species, protons and neutrons, seem to depart from homogeneity at the same temperature threshold, but the departure appears more sharply for protons than for neutrons, which attain higher maximum values than the neutrons.

Apparently, neutron-rich systems do not develop bubbles at the same temperature as the symmetric systems. The excess of neutrons inhibits the pasta formation until lower temperatures where protons manage to form the pasta. The released neutrons get distributed along the cell disrupting the pasta structure.

c. Classification of the pasta

The Minkowski functionals supply complementary information on the pasta structure from the early stage to the solid-like stage; see Table 6 for the classification of the pasta in terms of the curvature B and the Euler characteristic χ. The accuracy of this information is, however, conditional to the correct binning of the simulation cell. Tiny “voxels” (that is, a high density binning) may produce fake empty voids (artificial bubbles or tunnels) and, on the contrary, oversized voxels may yield a wrong structure of the system due to the lack of details. Therefore, some effort needs to be spent to determine the correct size for the voxels; Appendix G.5 summarizes this procedure.

The simulation cell was first divided into cubic voxels of edge length d = 2.35 fm. The Euler characteristic χ was computed for symmetric nuclear matter, according to Eq. (G.12) and the results are shown in Fig. 26. As seen in this figure, the Euler characteristic χ for symmetric matter shows three distinct regions as a function of temperature, one at T>2 MeV, one at T<0.5 MeV, and a transition one between these two.

Figure 26 shows that at T>2 MeV χ does not vary much although it attains different values depending on the density. At the low densities of ρ = 0.05 fm⁻³ and ρ = 0.06 fm⁻³ χ exhibits negative values which, according to Eq. (G.12), indicate that the nucleons are sparse enough to form tunnels and empty regions across the cell. At higher densities, however, χ becomes positive indicating that tunnels begin to fill forming cavities and isolated regions. This is confirmed by Fig. 27 which shows an inside view of the discretized nuclear matter at T = 2.5 MeV, x = 0.5 and for the four densities under consideration.

The complementarity of χ over other measures can be seen by comparing it to, for instance, the results of the Kolmogorov statistics (Section 2.2.2b). As seen in Section 2.2.2d, as the nucleons get distributed uniformly at high temperatures neither the 1D nor the 3D Kolmogorov statistic capture the qualitative difference between a tunnel-like and a cavity-like scenarios. Both landscapes may not exhibit noticeable heterogeneities, and thus, they appear to be essentially the same from the point of view of the Kolmogorov statistic.

The energy, however, distinguishes between aforementioned different scenarios. From the comparison between the caloric curves introduced in Section 2.2.2a and the current Euler characteristic, one can see that both magnitudes are density-dependent in the high temperature regime. While χ increases for increasing densities, the (mean) energy per nucleon diminishes (see Figs. 16 and 26). Thus, the less energetic configuration (say, ρ = 0.085 fm⁻³) appears to be a cavity-like (or small bubble-like) scenario from the point of view of χ (with voxels size of d =2.35 fm).

The Euler characteristic χ exhibits a dramatic change at T ≈ 2 MeV. This is associated to the departure from homogeneity at the early stage of the pasta formation, as already mentioned in Section 2.2.2c. It also agrees with the results from quantum models shown in Ref. [44]). Notice that the χ values for the examined densities join into a single pattern for T<2 MeV, in agreement with the behavior of the energy seen in Fig. 16.

It should be emphasized that although all the examined densities share the same χ pattern for T<2 MeV, their current morphology may be quite different. Figure 28 illustrates two such situations, see caption for details. It seems, though, that whatever the morphology, these are constrained to be equally energetic (see Fig. 16).

Extending the χ study for non-symmetric nuclear matter appears to confirm the complexity observed in Sections 2.2.2a and 2.2.2b. The global pressure does not present negative values for x = 0.4 and x = 0.3 at temperatures above the solid-like state. Neither noticeable bubbles nor other heterogeneities could be detected at the early stage of the pasta formation (say, T ≈ 2 MeV). We are now able to confirm these results through the χ functional. Figure 29 shows the Euler characteristic for two different densities and x = 0.3, 0.4, and 0.5.

Figure 29(a) shows three distinct behaviors of χ at ρ = 0.05 fm⁻³. The case of symmetric nuclear matter (x = 0.5) was already analyzed above. The curve for x = 0.4 appears left-shifted with respect to the symmetric case, in agreement with our previous observation that proton ratios of x<0.5 frustrate for a while the pasta formation (see Section 2.2.2b). In spite of that, the pattern for x = 0.4 achieves a higher positive value at lower temperatures than the symmetric case indicating that the tunnel-like scenario (χ<0) switched to a bubble-like or an isolated-structure scenario (χ>0). It can be verified from Fig. 30(a) that this is actually occurring at T ≈ 1 MeV; many isolated structures may be visualized in red, while no tunnels seem to be present in the blue region (see caption for details).

The curve of x = 0.3 in Fig. 29(a) does not change sign nor increases in value at lower temperatures, as opposed to the other two curves. The fact that χ<0 for all of the examined temperatures indicates that tunnel-like landscapes are the relevant ones. Figure 30(b) illustrates this scenario: the heavy tunnel-like region (highlighted in blue color) is mostly occupied by neutrons, as shown in Fig. 25, indicating that repulsive forces between neutrons dominate in a large fraction of the cell, and thus producing a positive global pressure all along the examined temperature range, as noticed in Fig. 16. Figure 29(a) also shows that χ for x = 0.3 decreases in magnitude with at lower temperatures; this corresponds to the departure from the spacial homogeneity detected in Fig. 23(a), and found by the Kolmogorov statistic and Euler functional to occur at T ≈ 1.5 MeV.

At the higher density of ρ = 0.085 fm⁻³ the behavior of χ is substantially different. Figure 29(b) shows the temperature dependence of χ for x = 0.3, 0.4 and 0.5; both isospin asymmetric cases are qualitatively similar to the symmetric one. An inspection of the voxels’ configuration (not shown) confirms that tunnels become relevant at very low temperatures; this can be checked from the configurations presented in Figs. 19 and 20.

Recall from Appendix G.5 that the chosen binning procedure minimizes the existence of spurious voids, that is, voxels that appear to be empty but they are really not. This criterion involves all the nucleons in the simulation cell. Thus, we are not able to extend the Minkowski analysis to protons or neutrons separately, as we did with the Kolmogorov statistic. In this context, this would produce an incorrect counting of voids (or tunnels).

In summary, at T ≈ 2 MeV the Euler characteristic χ of isospin-symmetric low-density systems (ρ<0.06 fm⁻³) shows drastic changes from negative to positive values indicating a transition from a void-dominated regime to one with bubbles and isolated regions. Higher density systems (ρ>0.06 fm⁻³), in spite of always having χ>0, also increase the value of their χ at this temperature, reaching a common maximum for all densities at T ≈ 1.5 MeV. This maximum corresponds to the formation of bubbles or isolated regions, and indicates the formation of the pasta near the solid–liquid transition; recall that the Kolmogorov statistic was able to detect the pasta formation since the bubbles or isolated regions stage.

For isospin asymmetric systems the low-density (ρ ≈ 0.05 fm⁻³) growth of χ is also observed but only for x = 0.4 and 0.5; systems at x = 0.3 have χ<0 at all temperatures. At higher densities (ρ ≈ 0.08 fm⁻³) the Euler characteristic is always positive for all temperatures.

Table 6 classifies the pasta according to the sign of χ and the curvature B, and it was a goal of the present study to extend this classification for isospin asymmetric systems, but our results indicate that this labeling becomes meaningless for the non-symmetric case. For a given temperature, the χ functional attains positive or negative values depending on the isospin content and the density of the system. In general, the excess of neutrons obscures the pasta structures for the protons and, thus, the early stage of the pasta formation (that is, the formation of bubbles or isolated regions) is not detectable. In spite of this, we observe that the system departs from homogeneity at T~1.5 MeV [see, for example, Fig. 23(a)].

d. Symmetry energy in the pasta

We now study the symmetry energy of nuclear matter in the pasta region. As mentioned before, at a given temperature the energy E(\rho \mathrm{,}T) shows three distinct behaviors as a function of the density: the pasta region for densities below 0.085 fm⁻³, the crystal-like region for densities above 0.014 fm⁻³, and an intermediate region in between the first two. In what follows we will focus on the symmetry energy in the pasta region.

Figure 31 shows the computed symmetry energy as a function of the temperature for four densities. The E_{sym} was computed through the fitting procedure outlined in Appendix D; an analysis of the goodness of the fitting can be found in Ref. [45].

Figure 31 shows several distinct regions for E_{sym}. In cooling, a liquid system with T>2 MeV starts with a low value of E_{sym}. Upon entering the T<2 MeV region and until T ≈ 1.5 MeV, the symmetry energy increases in magnitude while the liquid pasta is formed. Its value stabilizes in an intermediate value in the warm-to-low temperature range of 0.5 MeV<T<1.5 MeV where the liquid pasta exists. At T<0.5 MeV, when the liquid-to-solid phase transition happens within the pasta, the symmetry energy reaches its highest value. We now look at these stages in turn.

In the temperature range, T>2 MeV, Fig. 31 shows that E_{Sym} has higher values for higher densities. Since a similar relationship is maintained by the Euler characteristic χ (cf. Fig. 26), it is possible that there may be a connection between the symmetry energy and the morphology of the system. Remembering from Section 2.2.2c that at those temperatures higher densities are associated to cavity-like or isolated regions, and lowest densities with tunnel-like structures, it is possible that E_{sym} increases as the tunnels become obstructed and cavities or isolated regions prevail.

Between 1.5 MeV<T<2 MeV the symmetry energy varies in a way resembling the variation of the 3D Kolmogorov statistic, D, during the pasta formation stage. Figure 32 compares the variation of E_{sym} with that of D as a function of temperature for the four densities of interest, and finds a good match between both quantities; one can then conclude that the variation of E_{sym} can be also be associated to the changes in the morphology of the nuclear matter structure and, furthermore, that as the pasta is formed during cooling, the symmetry energy increases in magnitude.

The region between 0.5 MeV <T<1.5 MeV corresponds to the pasta structures filled with liquid nuclear matter, and E_{sym} changes at around T ≈ 0.5 MeV, at the same temperature at which the caloric curve and the Lindemann coefficient undergoes similar changes [cf. Fig. 15(b)], indicating the phase transition between a liquid pasta (T>0.5 MeV) and a solid pasta (T<0.5 MeV). Based on this, one can conclude that the symmetry energy attains its largest value in the solid crystal-like phase.

Finally, it is worth mentioning that the values of E_{sym} attained in this section are not directly comparable to those of Section 2.1 for the liquid–gas coexistence region and compared to experimental results. This is the case because, first, Fig. 9(a) was obtained at higher temperatures (2 MeV <T<5 MeV) and lower densities (ρ<0.03 fm⁻³), and, second, it corresponds to an homogeneous medium while the present one uses data from a pasta-structured system. In spite of these differences, it is reassuring that the values of E_{sym} obtained in the present study for the highest temperatures used, T>2 MeV, are within the range of values calculated in Section 2.1 for the highest densities considered in such study, namely 0.05 fm⁻³<ρ<0.06 fm⁻³.

2.3 Summary of nuclear matter properties

2.3.1 Summary of nuclear matter properties at intermediate temperatures

In Section 2.1 the bulk properties of nuclear matter were obtained along with their isospin dependence. It was determined that isospin symmetric and asymmetric matter can be self bound for low values of T. As T increased from 1 to 5 MeV their saturation densities varied from ρ₀=0.16 fm⁻³ to 0.12 fm⁻³ for isospin symmetric matter, and from ρ₀≈0.12 fm⁻³ to 0.09 fm⁻³ for matter with x = 0.4. Likewise, isospin asymmetry reduces the compressibility at saturation density by about 30% to 50% as x drops from 0.5 to 0.3 and as T increases from 1 to 5 MeV.

In Section 2.1 it was shown that nuclear matter can exist in liquid and gaseous phases and there are transitions between them. The liquid phase was identified at low temperatures for all isospin asymmetries studied by the ∪ shapes of the energy E(T\mathrm{,}\rho ). At sub-saturation densities the energy–density curves of x = 0.4 and 0.5 signaled a transition to a liquid–gas mixture phase, while x = 0.3 systems maintains phases and transitions among them at T = 1 MeV but become unbound at all densities at T = 3 MeV and higher. The addition of the x dependence to the phase diagram alloed its representation in Section 2.1 in the 3D space of T\mathrm{,}\rho \mathrm{,}x; the liquid–gas coexistence region extends, approximately, to densities of up to {\rho }_0/2, temperatures of up to 16 MeV, and isospin content as low as of 13%.

The symmetry energy of nuclear matter at intermediate temperatures at low densities was calculated in Section 2.1, and compared satisfactorily to experimental data. E_{Sym} is susceptible to the clusterization that occurs in phase transitions.

Besides corroborating previous studies, this work extends some of their results to other values of isospin content and non-zero temperatures. Findings that we believe are new are the temperature variation of the saturation density and compressibility for isospin asymmetric matter, certain details of the existing phases at x = 0.3 and 0.4, the extension of the phase diagram into the isospin axis, as well as a new procedure to estimate E_{sym}(T\mathrm{,}\rho ) from kinetic simulations. Finally, it is worth mentioning that the CMD model indeed helps to understand the role of isospin on several nuclear properties.

2.3.2 Summary of nuclear matter properties at low temperatures

In Section 2.2 the formation of the pasta-like structures in nuclear matter, its phase transitions and its symmetry energy. were investigated. In Section 2.2 the crystalline structure of NM was studied at zero and near-zero temperatures. Using scaling and CMD it was found that the potentials used yield a simple cubic structure. CMD indicated the departure from the crystalline state to a non-homogeneous pasta at sub-saturation densities; such structures are maintained up to T ≈ 1 MeV; this was evidenced by the use of curvature and Euler characteristic.

In Section 2.2 it was determined that the pasta exist in isospin symmetric and non-symmetric systems. Furthermore, solid-to-liquid phase transitions within the pasta were detected. Using the radial distribution function, the Lindemann coefficient, Kolmogorov statistics, Minkowski functionals the morphology of the pastas was succesfully linked to the phases and phase transformations. Nucleons can exist in liquid and crystal phases inside the pasta structures. It was also found that non-isospin-symmetric pasta depart from the classification of symmetric nuclear matter presented in Table 6.

The symmetry energy in pastas was studied in Section 2.2.2d. connection to the morphology of the pasta and to the phase transitions. E_{Sym} showed different values as a function of the temperature and density, i.e., in the different phases of nucleons inside the pasta. E_{sym} is associated to the morphology of the nuclear matter structure, and it attains its largest value in solid pastas.

In conclusion, classical molecular dynamics simulations show the formation of pastas in isospin symmetric and non-symmetric systems. The computational tools developed and applied, although not perfect, demonstrated their usefulness to detect the in-pasta phase transitions first seen in Ref. [19], and to extend the calculation of the symmetry energy of Ref. [33] to lower temperatures, and connect its value to the structure and thermodynamics of the neutron-rich pasta.

The bulk properties of nuclear matter, its phases and phase transitions, along with the behavior of E_{sym} is indicative of similar phenomena present in neutron star matter (NSM). Neutron star matter, however, is a system different than pure nuclear matter; the connection of NM and NSM will be investigated in Section 3, and the behavior of NSM in Section 4.

3 Electron gas: Connecting nuclear matter with neutron star matter

As explained in the introduction (cf. Section 1), nuclear matter is connected to neutron star matter by means of an electron gas that fills all space between nucleons. In addition to neutrons and protons, neutron star crusts also contain electrons which fill the space between nucleons. That is, neutron star matter can be thought of as nuclear matter embedded in an electron gas.

Many authors (see e.g. [15, 20, 39, 46]) believe the pasta is nothing but frustrated structures formed by a competition between nuclear and Coulomb forces; obviously, as seen in the section devoted to nuclear matter (Section 2), this is not necessarily the case; nuclear repulsive nn and pp and attractive np forces form pastas as well. Thus, before embarking on a full study of neutron star matter, we find convenient to dissect the effect an electron gas has on the nuclear matter pastas.

As it will be explained below, to introduce an electron gas in a CMD study of NM pasta, the addition of an screened Coulomb potential is necessary. This section focuses on the effects the strength and range of such interaction have on the morphology of the pastas.

At this point it is worth mentioning that, as it will be seen in Section 3.2, the pasta-like structures calculated without the proper screening length introduce spurious effects due to the size of the cell and form only one “pasta like” structure per cell; such “pastas” are not a true nuclear pasta property, but their use in the analysis does not invalidate the study of the effect of the strength of the screened Coulomb potential.

3.1 The strength of V_C

To mimic neutron star matter it is necessary to include the Coulomb effect of an electron gas in the nuclear matter studied in Section 2. CMD is specially well suited for this study as it provides a microscopic view of the structures formed. The electron cloud introduces an screening effect on the Coulomb potential of the protons, which in turn modifies the pasta structures seen in Section 2.2. The techniques to introduce the Coulomb interaction of an electron gas into CMD are presented in Appendix A.1.2.

The effect of an electron gas on cold nuclear structures has been studied before. In 2003 a static liquid–drop model was used, and it was found that its main effect was to extend the range of densities where bubbles and clusters appear [47]. Later, in a 2005 study a density functional method was used (at zero temperature) to find that the density region in which the pasta exists becomes broader when the electron gas is taken into account [48]. More recently, Monte Carlo simulations, which approximated the Coulomb interaction via an Ewald summation [15, 20, 49, 50], determined the lowest proton fractions that are compatible with β equilibrium in neutron star crust environments (at T = 1.0 MeV).

It is worth mentioning that a system with a coulomb term of just one sign (as would be the case of nuclear matter if protons are considered as carrying charge) is unstable and thermodynamically non additive. But the inclusion of the electrons which render the coulomb term of finite range via screening transforms the system into stable and thermodynamically additive, and as such the energy will be additive as well as the entropy.

In this section we study the effect the strength electron gas has on the formation of the nuclear pasta at non-zero temperatures through the calculation of the pasta structures with and without effect of the electron gas as well as with softened interactions which – of course – do not exist in nature. Full details of this study can be found in Ref. [51].

3.1.1 Simulation procedure

To study nuclear matter embedded in a degenerate electron gas, the electron gas density is first taken as to produce an overall neutral (β-equilibrated) system. The electron gas effectively screens out the positive charges of the protons, resulting in a more local effect. This is taken into account through the use of a Thomas–Fermi screened Coulomb potential, V_C=(q^2/r){\text{e}}^{-r/\lambda }; this potential is introduced in Appendix A.1.2.

We focus on CMD simulations of isospin symmetric systems of 1000 protons and 1000 neutrons, and isospin asymmetric systems comprised of 1000 protons and 2000 neutrons, i.e., of x = 0.3. Nucleons are placed in cubical boxes with densities between 0.01 fm⁻³≤ρ≤ρ₀, temperatures between 0.1 MeV ≤T≤1.0 MeV and periodic boundary conditions; these conditions correspond to semi-frozen states where pastas are known to exist in NM. In particular, we produce pastas as in the previous sections, and compare their structures obtained with and without the electron gas as well as with varying strengths of it. In summary, pasta structures were obtained at subsaturation densities (0.015 fm⁻³≤ρ≤0.072 fm⁻³) and low temperatures (0.1 MeV≤T≤1.0 MeV) for symmetric (x = 0.5) and asymmetric (x = 0.3) cases. Altogether 200 simulations were carried out per each combination of \left\{\rho \mathrm{,}T\mathrm{,}x\right\}.

To see the effect of the electron gas, each case was “cooked” repeatedly with a screened Coulomb potential with a varying amplitude, i.e., with \alpha V_C, where 0\le \alpha \le 1. For a fair comparison, all corresponding cases with different values of α were produced with identical initial conditions of \left\{x\mathrm{,}\rho \mathrm{,}T\right\} and started off from the same initial random configuration. As an illustration, Fig. 33 shows two corresponding structures obtained with and without the electron gas.

We use several tools for inspecting the pastas produced. At the densities and temperatures selected, the nucleons do not have large mobilities and self-bound cumuli can be detected with a simple “Minimum Spanning Tree” (MST) algorithm, which finds particles closer to each other than a clusterization radius set to 3.0 fm, and yields information about the fragment multiplicity. Again, we limit the cluster analysis to the central cell ignoring the neighboring cells of the periodic boundary condition, and thus eliminate possible clusters that may continue into other cells and can be, in principle, infinite in size.

At a microscopic level, the dynamics of the nucleons can be quantified through their average displacement as a function of the “time” steps of the simulation. Likewise, the microscopic stability of the clusters can be gauged through the “persistency” [28] which measures the tendency of members of a given cluster to remain in the same cluster. Another interesting descriptor is the isospin content x of each cluster produced. Other global characterization tool is the pair correlation function, and the Minkowsky functionals.

3.1.2 The effect of the strength of αV_C

The presence of an electron gas affects the fragment size multiplicity. Figure 34 shows the cluster multiplicity observed in 200 configurations of asymmetric matter with and without the Coulomb interaction, at x = 0.3, ρ = 0.015 fm⁻³ and T = 0.1 and 1.0 MeV. The increase of the number of fragments of sizes A\gtrsim 500 as α goes from α = 1 to α = 0 underline the role the electron gas has on determining the mass distribution.

The change of the inner structure can be quantified through the use of the radial distribution function. Figure 35 shows the g(r) of symmetric structures at x = 0.5, ρ = 0.072 fm⁻³ and T = 0.1 MeV and for varying strengths of the Coulomb potential, namely α = 1 (full Coulomb), 0.8, 0.2, and 0 (without Coulomb). The electron gas appears not to have an effect on the nearest neighbor distances.

The isotopic content on the cluster is also affected by the presence of the electron gas. Figure 36 shows the x content of the clusters formed in 200 configurations of asymmetric matter with x = 0.3, ρ = 0.015 fm⁻³ and T = 1.0 MeV. The electron gas clearly enhances the isotopic content from x = 0.3 up to x ≈ 0.4, as well as the maximum fragment size which grows up to A ≈ 800.

As suspected in Ref. [48], the interactions between protons and electrons are responsible for the rearrangement of protons in varying degrees. Figure 37 shows the RMS displacement of nucleons with and without Coulomb during the evolution of the ergodic process. Each curve is made of 200 points, each representing the average displacement with respect to the original configuration of the 2000 nucleons in systems with x = 0.5, ρ = 0.015 and 0.72 fm⁻³ at T = 0.1 and 1.0 MeV). At low temperatures Coulomb enhances nucleon mobility more than at high temperatures.

Figure 38 shows the “persistence” (i.e., the tendency of nucleons to remain in the same cluster), with and without Coulomb, as it evolves through 200 configurations of the stationary ergodic process, measured with respect to the original configuration of 3000 nucleons with x = 0.3, ρ = 0.015 and T = 1.0 MeV. As with the mobility, Coulomb enhances the transfer of nucleons thus decreasing the persistence as shown by the two curves.

The Euler characteristic and mean curvature introduced in Appendix G.4 can demonstrate the effect of the electron gas on the shape of the pastas. The mean curvature and the Euler characteristic were calculated from the digitized polyhedra constructed for each of the 200 nuclear structures obtained with and without Coulomb at each value of \left\{x\mathrm{,}\rho \mathrm{,}T\right\}. The values of the curvature varied approximately from -6000 to 12000 and those of the Euler characteristic from -700 to 1000. Table 1 shows the curvature–Euler number as percentages of the maximum values obtained, and Fig. 39 shows their location on the curvature–Euler plane as well as their changes in position as the Coulomb strength is diminished; the standard deviations are of the size of the points used in the plots.

3.2 The screening length

A crucial parameter of V_C is the screening length, λ. Here we study the effect this parameter has on the properties of the pasta within the framework of CMD. In particular, we study pasta properties with values of λ ranging from λ = 0 to 20 fm with densities ρ = 0.005, 0.03, 0.05, 0.08 fm⁻³, and λ = 30 fm and λ = 50 fm with ρ = 0.005 fm⁻³, where “gnocchi” are formed. The cut-off length was chosen at r_c=λ.

3.2.1 Critical screening length

As observed in previous works [41, 52], in absence of any Coulomb interaction (equivalent to λ = 0), pasta-like structures exist, although only one per cell, which indicates that the structure is limited by the periodic boundary conditions imposed on the box. When there is Coulomb interaction, however, the competition between opposing interactions gives rise to a characteristic length which shapes the pasta structures. By increasing the value of λ, starting from 0 fm, we determine the critical λ_c at which the pasta structures go from the artificial one structure per cell, to the more realistic case of more than one structure per cell.

A way to determine λ_c is by inspecting the nucleons pressure, i.e., the mechanical stability of the structures formed. The pressure is computed by the virial Eq. (A.3) with \langle \rho T\rangle =N{k}_{B}T/V, applying it only to the nucleons and not to the electron gas.

Negative pressure indicates unstable structures kept in place by the periodic replicas. In this case the overall effective interaction is mostly attractive and periodic boundary conditions play a major role in shaping the structure. The pastas formed under negative pressure are artificial and can only exist under periodic boundary conditions (see Refs. [41, 53]). Figure 40 shows that for all λ<10 fm the pressure is negative.

Positive nucleon pressure, on the other hand, indicates that the structures form due to the physical interactions and not to the boundary conditions. Indeed in the positive-pressure regime the density fluctuations are of smaller sizes than the size of the cell, and the morphology of the structures changes drastically. Figure 40 shows that, for those conditions, the pressure becomes positive for λ>10 fm.

λ_c can also be determined by the effect of λ on the morphology of the pastas, as quantified by means of the Minkowski functionals. The plot in Fig. 41 shows the surface, mean breadth and Euler number for ground-state structures (T = 0) as a function of λ. Between λ = 7 fm and λ = 10 fm all three of the Minkowski functionals change drastically before reaching well defined values.

The changes of the Minkowski functionals indicate changes in the pasta structures. As stated in Table 6, the lasagna and spaghetti have Euler characteristics χ ≈ 0, the gnocchi {\chi }_{gn}\approx 2, and a system with N_gn gnocchi \chi \approx 2\cdot N_{gn}, and the surface changes correspondingly as well. The mean breadth, on the other hand, should be positive for spaghetti and gnocchi, zero for lasagna, and negative for tunnels. The structures obtained with λ<10 fm do not exhibit the full range of values the Minkowski functionals should have.

3.2.2 Evolution of the pastas as a function of λ

Figures 42–44 show ground-state structures for λ = 0, 10 and 20 fm, and three different densities. While the structures at λ = 0 and 20 fm are recognizable pasta structures, those obtained at λ = 10 fm are more exotic, probably due to the very rough energy landscape at this value of λ.

More information about the λ_c can be obtained from the size of the gnocchi obtained for various values of λ. Figure 45 shows that at ρ = 0.005 fm⁻³, there is a single gnocco for λ<10 fm, but several for λ between 15 fm to 20 fm. Likewise, Fig. 46 indicates that the average size of the gnocchi gets reduced while its number increases reaching a stable value for λ≥20 fm. This larger number of clusters explains the changes of the Minkowski functionals.

Figure 47 examines the λ = 0 and 20 fm cases at ρ = 0.05 fm⁻³ with the radial distribution function, g(r). It is clear that with Coulomb (λ = 20 fm) the pastas have a more ordered neighboring structures than without Coulomb (λ = 0).

3.3 Summarizing the electron gas

The effect of the electron gas on symmetric and neutron-rich matter was studied at low densities and temperatures by varying the Coulomb interaction strength and screening lengths. Its effect on the fragment size multiplicity, the inter-particle distance, the isospin content of the clusters, the nucleon mobility, and on the modification of the topological shape was studied.

The most general result is the existence of the nuclear pasta structures even without the presence of the electron gas. As seen in Section 2.2, the existence of the pastas, which are usually associated with the presence of the competing interactions of the long range Coulomb potential and the short range nuclear force, can exist even in the absence of the electron gas (i.e., when \alpha \to 0). These pseudo-pastas are the result of the competition between the attractive V_{np}, and the repulsive V_{nn} interactions; this same effect has been seen in other potentials [50] as was reported in Ref. [41]. It must be remembered that, if not calculated with the proper screening lengths, the pseudo-pastas contain effects due to the cell size.

The main effect of the strength of the Coulomb interaction of an electron gas is to allow the system to display true pastas structures, and distribute matter more and form less-compact objects. At a microscopic scale, Fig. 35 shows that varying the strength of the Coulomb interaction does not change the inter-particle distance, but certainly decreases a bit the x content of the fragments as can be seen in Fig. 36. In agreement with this, the increment of nucleon mobility produced by Coulomb (cf. Fig. 37) gets reflected in a reduction of the persistence (cf. Fig. 38), and of the curvature and Euler characteristic of the structures, “compactifying” symmetric matter into spaghetti and lasagna type structures as shown in Fig. 39, although the effect is less pronounced in for x = 0.3 (bottom panel).

Likewise, our study of the effect of the screening length of the Coulomb interaction determined that there is a critical screening length λ_c at which the structures change drastically. For , the Coulomb interaction is barely acting, a negative pressure produces only one structure per cell indicating that the structures formed are due to finite size effects. For \lambda >{\lambda }_c, on the other hand, the pressure becomes positive and the systems present density fluctuations smaller than the cell, and the morphology of the structures stabilize and cease to depend on λ. For the Pandharipande potential λ_c lies between 10 fm and 15 fm depending on the density.

The previous findings have a practical consequence for size of the systems studied computationally. How large should the system studied be to avoid finite size effects? The scales determining the phenomena studied are the range of the potential and the size of the cell used. In our case the range is set by the Coulomb potential which is given by the screening length, λ_c. As summarized before, for the densities studied λ_c lies between 10 fm and 15 fm; once this limit is surpassed, systems with larger number of particles will merely reproduce the structures found in smaller systems. This can be observed in the computationally study of Horowitz using 50 000 particles [54].

4 Neutron star matter

In this Section we extend the study of nuclear matter of Section 2.2 to the realm of neutron star matter. We use CMD to study the existence of pasta-like structures and their possible phase changes. We focus the investigation at saturation and sub-saturation densities, low temperatures, and proton fractions in the range of 10% to 50%, and we also study the behavior of E_{Sym} in the different structures found.

We continue using CMD as presented in Appendix A.2.3. Although the parameters μ_r, μ_a, μ₀ and V_r, V_a, V₀ were first set by Pandharipande for cold nuclear matter [42], a recent improvement [45], here named New Medium, reproduces the cold nuclear matter binding energies more accurately and, thus, is used in this study of NSM. The corresponding values are summarized in Table 2. Figures 79 and 80 contrast these potentials with those of Pandharipande Medium potentials.

As explained in Section 3, the neutron star matter contains protons and neutrons embedded in an electron gas, which introduces a screening effect on the Coulomb potential of the protons. The screened Coulomb potential is included in these calculations using the Thomas–Fermi approximation as explained in Appendix A.1.2. It might be appropriate at this point to emphasize that the exponential cut-off renders the Coulomb effective interaction, short ranged. In this way the energy and the entropy is additive (i.e., energy scales with the number of particles).

4.1 Symmetric neutron star matter

We first study the case of symmetric neutron star matter, i.e., with x = z/A = 0.5, focusing on the caloric curve and the binding energy. In most of the cases presented in this section the total number of nucleons in the primary cell was N = 4000, interacting through the New Medium model (cf. Appendix A.1.1).

4.1.1 The caloric curve

For starters we calculate the internal energy of symmetric neutron star matter. This involves averaging the kinetic and potential energy of each nucleon in the system, cf. Appendix A.2.2. For full details see [55].

As explained in Section 2.2, the caloric curve can be used to detect phase changes. The caloric curve is the relationship between the internal energy and the temperature, and changes in the slope of the E–T curve can be used as indicators of phase transitions. It is convenient to remember that, as seen in Section 2.2, nuclear matter maintains a liquid-like structure for temperatures larger than about 2.0 MeV, and around T = 0.5 MeV it transforms into a pseudo-pasta with nucleons frozen much like in crystalline structures. The results for neutron star matter are somewhat different.

Figure 48(a) presents the caloric curve for the case of ρ = 0.05 fm⁻³, and Fig. 48(b) shows the caloric curves for densities ρ = 0.02, 0.05, 0.10 and 0.16 fm⁻³ in the extended temperature range of 0.2 to 4 MeV. As it can be seen more clearly in Fig. 48(a), the slope of the internal energy exhibits a change at around T\simeq 1.5\text{MeV}, and a sharper change at T\simeq 0.5\text{MeV}. Although these jumps are not as pronounced as those found for NM, they happen at around the same temperatures and can be taken as the same type of transitions.

We claim that the smooth change at T\simeq 1.5\text{MeV} signals the onset of the topological phase transition, (i.e., a bubble appears). A difference with NM is that this change in slope appears more pronounced for smaller densities (ρ≤0.05 fm⁻³), and tends to disappear for larger densities. As this result is different than the NM case, we believe it is due to the presence of the electron gas.

On the other hand, the sharper jump observed at T\simeq 0.5\text{MeV} both in Figs. 48(a) and (b) is maintained throughout the density range studied and, in fact, becomes more pronounced for larger densities. We believe this discontinuity in the derivative of the E–T curve signals the change from amorphous pasta to crystalline pasta already identified in Section 2.2 for NM with the Pandharipande Medium model.

The noticeable differences between the caloric curve of nuclear matter and that of neutron star matter are shown in Fig. 49. The figure compares the NSM results (continuous lines) against the NM curves (dashed lines) for the cases of ρ = 0.05, 0.06, 0.07, and 0.085 fm⁻³ and temperatures from 0.2 to 4 MeV. The NSM lines remain separated and parallel throughout the range of explored temperatures, whereas the NM curves merge with one another at low temperatures. Noteworthy, the NM highest density curve (ρ = 0.085 fm⁻³) corresponds to the lowest energy, while the opposite is true for NSM where the higher density curve has the higher energy of all the curves shown.

Figure 50 shows more examples of the caloric curve for several densities. Each of these densities exhibit a discontinuity in the energy at certain temperatures, a signal of a first order phase transition. This transition can be confirmed and further characterized as a solid–liquid phase transition by looking at the Lindemann coefficient. The Lindemann coefficient for ρ = 0.05 fm⁻³ as a function of temperature, can be seen in Fig. 15(b), along the energy. This figure shows that the discontinuities in Lindemann coefficient and in energy are at the same temperature. These two factors are, effectively, the signature of a solid–liquid phase transition.

4.1.2 Energy–density

As seen in the case of NM (cf. Section 2), the E–ρ curve can signal the saturation density (i.e., the minimum of the curve), bound and unbound states (positive or negative energies), the range of densities of the liquid phase (span of the “∪” shape), and changes of phase, among others. For the case of NSM, however, the situation is different.

Figure 51 shows the energy as a function of the average density at T = 0.1, 0.2, 0.5, 1.0, and 2.0 MeV. At a difference from the usual ∪ shape of symmetric nuclear matter, the NSM curves do not exhibit absolute minima, indicating that NSM at these temperatures and densities does not have an equilibrium point nor a saturation density. Furthermore, at most densities the NSM system appears to be unbound, except for T≤1.0 MeV and ρ≤0.05 fm⁻³. The difference between NM and NSM is the effect of the cloud of electrons which, with its exponential cut-off of the order of 20 fm, introduces a short range Coulomb screening which produces smaller aggregates of nucleons, as seen in Section 3.2.

4.1.3 The Minkowski functionals

As before, the morphology of the pasta structures will be studied with the Minkowski functionals, namely the volume, surface area, Euler characteristic χ, and integral mean curvature B. The calculation of the Minkowski functionals requires the binning of nucleons into “voxels”, as described in the Appendix G.5.

Figure 52 shows the Euler characteristic as a function of the temperature. χ has a clear change of behavior at T\simeq 1\text{MeV}. Although this happens at a lower temperature than its NM counterpart (see Fig. 26), it also appears to be associated to the early stage of the pasta formation. Furthermore, the χ values for the examined densities almost join into a single pattern for T\le 0.5\text{MeV}, much like those of NM at the same temperature but less pronounced; this could also be related to the change of slope found in the caloric curve at T\simeq 0.5\text{MeV}.

A look at the sign of χ can yield information about the morphology of the structure. Figure 52 indicates that the lower density systems attain negative values of χ at, say, T>1 MeV. Since values of χ<0 corresponds to cases where the number of tunnels overcome the number of voids and isolated regions, see Eq. (G.12), it seems that low density configurations tend to be more cavity-like (more tunnels), but as density increases tunnels fill up yielding more compact structures.

When we look at the Minkowski functionals, particularly the Euler characteristic and the mean breadth, we can see that there is again a critical temperature at which both the Euler characteristic and the mean breadth show a sharp transition. We show, as an example, these magnitudes as a function of temperature for density ρ = 0.05 fm⁻³ in Fig. 53. As this transition is signaled by morphological observables, we conclude that this transition is morphological.

These signals of a solid–liquid phase transition (energy and Lindemann’s coefficient discontinuity) and morphological transition (Minkowski functionals discontinuity) point at the same transition temperature, as can be seen in the phase diagram of Fig. 54. This means that as the systems are cooled down at fixed volume, they undergo a thermodynamical and a morphological phase transition, and they do so at the same temperature.

In summary, and comparing to NM, the introduction of the coulomb screened potential appears to smooth out the phase transitions. The Euler functional χ experiences a somewhat smooth change along 1–2 MeV, in correspondence with the energy changes mentioned in Section 4.1. Thus, the pasta forming process may be located at this temperature range, in a similar fashion as in nuclear matter systems. At temperatures below T\simeq 0.5\text{MeV} the (inner) nucleons in the pasta structure freeze into a solid state.

4.1.4 The radial distribution function

Here we use the radial distribution function to explore the phase transformations detected by the caloric curve. Figure 55 shows g(r) for systems with ρ = 0.02 and 0.085 fm⁻³; the lower density case shows more pronounced nearest-neighbor peaks than the high density one, and it corresponds to a more crystalline phase. Looking at Fig. 49, we see that the highest internal energy is associated to a more regular distributions of nucleons within the pasta regime.

In Fig. 56 we show the radial distribution function for three different densities: ρ = 0.03 fm⁻³ (spaghetti), ρ = 0.05 fm⁻³ (lasagna) and ρ = 0.08 fm⁻³ (tunnels), just above and below the transition temperature, as well as a snapshot of the system at the high temperature phase. Since the first peaks (corresponding to the nearest neighbors) are at the same position regardless of the temperature, we conclude that the short range order is present both above and below the transition. However, the peaks for third and higher order neighbors, distinctive of solid phases, disappear as the temperature is increased through the transition.

On top of the disappearance of the long range order characteristic of solids, another feature becomes evident from Fig. 56. As temperature increases through the solid–liquid transition, a very-long range modulation in the pair correlation function survives. This very long range ordering is characteristic of the pasta phases. In Fig. 57, a visual representation of the spatial configuration for ρ = 0.05 fm⁻³ is shown, for temperatures both below and above the transition. In it, we show that not only the solid phase has the usual pasta shape, but the liquid phase preserves it. Below the transition, we have frozen pasta. Just above it, nucleons may flow but confined to a certain pasta or pasta-like structure. In Fig. 58 the usual lasagna is observed, but also intertwined lasagnas and other structures that are not of the usual pasta type.

Figure 59 shows the proton arrangements for three density situations. Figure 59(c) indicates that the ρ = 0.085 fm⁻³ situation is highly homogeneous inside the occupied regions. These occupied regions split into smaller pieces, according to Figs. 59(b) and (a), forming lasagnas and gnocchis, respectively. Thus, as the density diminishes, the structures break into smaller pieces, and the internal energy decreases (see Fig. 51).

Likewise, at ρ = 0.02 fm⁻³, Fig. 55 shows that g(r) has strong neighbor correlations at r≤6 fm⁻³ and almost no correlation [g(r)≈0] for r≥8 fm⁻³. That is, strong correlations exist between neighbors belonging to the same gnocco, and no correlation at large separation distances between gnocchis [see Fig. 59(a)]. For the larger density, ρ = 0.085 fm⁻³, g(r) tends to 1 as expected for homogeneous systems [cf. Fig. 59(c)].

4.1.5 Summary of Section 4.1

We may summarize the section as follows. The internal energy for symmetric neutron star matter increases monotonically for either increasing temperatures and densities; it does not appear to have a saturation point. These effects do not occur for nuclear matter, and they are due to the Coulomb interation of the electron gas as can be seen clearly by comparing Fig. 48(b) to the nuclear matter case. Likewise, the introduction of the Coulomb potential affects the system morphology; as the density diminishes, the pasta structures split into smaller structures; these are true pasta structures, as opposed to the one-structure-per-cell pseudo-pastas observed in NM.

4.2 Non-symmetric neutron star matter

We now turn to a study of NSM with proton fractions more like those in in neutron crusts. All simulations in this Section consist of systems with N = 4000 nucleons under periodic boundary conditions and interacting through the New Medium potential and, whenever used, the binning distance is d = 2.35 fm (see Appendices A and G for details).

4.2.1 The internal energy

We start the study of non-isospin-symmetric NSM by exploring the dependence of the internal energy on the isospin content. Figures 60(a) and (b) show the isothermal energies as a function of x for systems with densities ρ = 0.04 fm⁻³ and ρ = 0.085 fm⁻³. Each of the curves exhibits a ∪ shape indicating the existence of a minimum of the energy at a certain value of x. At ρ = 0.04 fm⁻³ the observed minima appear to shift from x = 0.4 at T = 0.2 MeV to x = 0.3 at T = 2.0 MeV, while at ρ = 0.085 fm⁻³ the minima remain at x = 0.3 at all temperatures. This finding indicates that systems with the freedom of exchanging their content of neutrons and protons would favor a specific isospin ratio x depending on the local temperature and density of the system.

Complementary information can be obtained from the density dependence of the energy. Figure 61 shows the corresponding curves for x = 0.2, 0.4 and 0.5 at T = 0.2 MeV and 1.0 MeV. It is interesting to notice that while the x = 0.4 and 0.5 curves appear to have similar monotonically decreasing behaviors at all densities, the x = 0.2 case deviates at low densities; as explained in Ref. [45], this is due to the ∪ shape of the x dependence of the energy seen in Figs. 60(a) and (b).

4.2.2 The associated morphology

The Euler characteristic χ and the radial distribution function g(r) are now used to study the evolution of the pasta structures.

As explained in Appendix G and Section 2.2 for NM, the Euler characteristic χ can indicate changes of structure. Figure 62 plots χ as a function of temperature for four proton fractions and two densities. Changes of behavior can be seen at around T ≈ 1 MeV. Figure 62(a) shows that at ρ = 0.04 fm⁻³ and T≥2 MeV, the system appears dominated by voids and tunnels independent of the isospin content, while at T≤1 MeV more compact objects appear at all x, except for x = 0.2 which disperses even more. Figure 62(b) shows that at the higher density (ρ = 0.085 fm⁻³) the system attains a more compact structure at all isospin contents and for T≥1 MeV, becoming less pronounced at lower temperatures (T≤1 MeV).

Figure 63 shows the position of the protons for the cases of x = 0.2, 0.4, and 0.5 at T = 0.2 MeV and ρ = 0.085 fm⁻³. These structures correspond to the T = 0.2 MeV points of the χ curves in Fig. 62(b). Notice that no major changes in the percentage of void to filled volumes are observed as it is difficult to appreciate such changes in compact structures.

The situation is different at lower densities. Figure 64 shows the protons for the cases of x = 0.1, 0.2, 0.3, 0.4, 0.45, and 0.5 at T = 0.2 MeV and ρ = 0.04 fm⁻³. Comparing Figs. 64(b), (d) and (f) to the T = 0.2 MeV values of χ in Fig. 62, one can see the how the voids and tunnels decrease as x goes from 0.2 to 0.5.

The radial distribution function is shown in Fig. 65 for two representative proton fractions and densities. No qualitative differences can be distinguished between the x = 0.2 and the x = 0.4 cases. Neither can be distinguished between these and the symmetric case shown in Fig. 55 for ρ = 0.085. Thus, the overall nearest-neighbor distances are not affected by the proton fractions.

The lack of change of g(r) in Fig. 65(b) can be understood by looking at the structures formed. Figure 63 shows the location of protons for three formations at ρ = 0.085 fm⁻³; the hollow spaces correspond to regions occupied by neutrons and voids. The number of hollow regions increases as x decreases, but maintaining the average inter-particle distances without much change.

The radial distribution functions for the lower density ρ = 0.02 fm⁻³ are shown in Fig. 65(a). Both profiles correspond to gnocchi-like structures, similar to those in Fig. 59(a) for symmetric matter (not shown). Figure 65(a) and Fig. 55(a) indicate that as x diminishes, the profiles become somewhat smoother.

Figure 66 splits Fig. 65(a) into the corresponding g(r) for protons and neutrons. Comparing Figs. 66(a) and (b) it is clear that the smoothing of Fig. 65(a) for decreasing x is related to the distribution of nucleons within the gnocchi structure (at, say, r≤6 fm). Specifically, the proton–neutron and neutron-neutron distributions smoothen when the fraction x diminishes.

We may further examine the simulation cell for the intermediate density ρ = 0.04 fm⁻³. Figure 64 shows the corresponding proton structures (neutron not shown) for a sequence of fractions x<0.5. As x decreases, the pasta splits into smaller pieces. The structures for nearly symmetric neutron star matter (x\simeq 0.5) are lasagna-like structures, while the low proton fractions (x<0.3) form gnocchi-like structures. The spaghetti-like structures appear in between (x\simeq 0.4).

Comparing Fig. 60 and Fig. 64 confirm that the changing morphology is due to the ∪-shape pattern for the internal energy as a function of x. The spaghetti-like structure achieves the minimum energy, while the lasagnas and the gnocchis correspond to higher energies on either branch of the ∪-shape.

In summary, two major effects appear when departing from symmetry. A morphological re-arrangement of the pasta structures occurs, attaining some kind of fragmentation as the proton fraction diminishes. Structures with the lowest x attain the gnocchi structure, but, the gnocchis themselves experience (inner) topological changes during its formation process. The resulting energy level surpasses the one at x~0.3 in agreement with the observations made in reference to the energy ∪-shape dependence on x.

4.3 The symmetry energy

In Sections 4.2 and 4.2.2 we learned that morphological changes taking place for decreasing x are related to the ∪-shape dependence of the energy on x shown in Fig. 60. These changes are bound to modify the symmetry energy too. In this section we calculate E_{Sym} according to the procedure detailed in Section 4.3 and Appendix E.

Figure 67 shows the fittings of the internal energy, Eq. (E.1), as a function of x; in this case the term of \mathcal{O}({\alpha }^4) was neglected for simplicity, see Appendix E for details. Data appears to be fitted appropriately by a quadratic profile, however, the low temperature curves shown in Fig. 67(a) depart from x~0.4 data. This is somewhat an \mathcal{O}({\alpha }^4) discrepancy since a noticeable improvement can be obtained if an α⁴ term is added into the fitting procedure (not shown).

Figure 67(b) shows that the fittings for ρ = 0.085 fm⁻³ match better the CMD data than in the case of ρ = 0.04 fm⁻³ [Fig. 67(a)]. Recall that the higher density corresponds to more compact proton structures (Fig. 63) than at ρ = 0.04 fm⁻³ (Fig. 64). Thus, as already noticed in Sections 4.2 and 4.2.2, more “fragmented” structures (say, spaghettis instead of lasagnas) flatten the right-branch of the U-shape pattern in Fig. 67(a), departing from a seemingly quadratic profile.

We now proceed to compute E_{Sym} for neutron star matter. For comparison, Fig. 68 shows E_{Sym} for nuclear matter (similar to in Fig. 31) and for neutron star matter. It is reassuring that the values of E_{Sym} for NM and NSM are approximately equal at T\gtrsim 2\text{MeV}.

Figure 68(b) shows a change in the behavior of E_{Sym} for NSM as compared to NM at temperatures below 1.5 MeV. An important observation is that for T<1.5 MeV, E_{Sym} changes slope at the pasta regime, and decreases with lower values of T; this is the opposite behavior of the NM case. In broad terms, this is due to the changes in pasta structure and increase of gnocchi multiplicity that occur at low densities. See [45] for complete details.

More recently, analytical expression of the symmetry energy, as well as the connection of the pseudo-pastas with the Minkoswki functionals were studied using machine learning technology [56].

4.4 Neutrino transport properties

A property of the NSM pasta of interest for neutron stars is the so-called neutrino opacity, which refers to the neutrino-pasta scattering. Neutrino opacity is of crucial importance for the evolution of neutron stars as they cool down by means of neutrino emission. Here we study the dynamics of the neutrino opacity of the heterogeneous matter at different thermodynamic conditions. For different densities, proton fractions and temperature, we calculate the very long range opacity and the cluster distribution.

4.4.1 Structure function

As can be expected, the absorption of neutrinos by the neutron star crust depends on the structure of the nuclear pasta existing in the crust. The cross section for neutrino scattering by a structured medium is related to the cross section for scattering by a uniform medium by

\begin{array}{}& \frac{\text{d}\sigma }{\text{d}\mathrm{\Omega }}={\left(\frac{\text{d}\sigma }{\text{d}\mathrm{\Omega }}\right)}_{\text{uniform}}\times S(q)\mathrm{,}\end{array}

where S(q) is the static structure factor of nuclear pasta [15], and q is the relative momentum. The structure factor of the system is directly related to the radial correlation function g(r) via a Fourier transform. Here we calculate S(q) according the prescription developed for systems in periodic cells in Refs. [10, 57].

In Fig. 69 we plot the height of the low momenta peak (\lambda \gtrsim 13\text{fm}) as a function of the temperature for the four structures of Fig. 56, all at \rho =0.05{\text{fm}}^{-3}. Different behaviors can be seen at T>0.5 MeV in which all curves are practically identical, and at T<0.5 MeV where the structures have different absorption peaks for low momenta. This is in agreement with Section 4.1, where Fig. 56 showed that g(r) can develop a very-long range ordering characteristic of some pasta phases. This very long-range order is responsible for a peak at very low momentum k (~10 fm wave-length) in S(k).

At high temperatures the nucleons are rather uniformly distributed and no structure is evidenced by S(k): the height of the peaks tend to 1, the value for homogeneous systems. As the temperature is decreased, a peak at low momentum develops. The transition described before manifests in Fig. 69 as the vanishing of fluctuations below the transition temperature T\lesssim 0.5\text{MeV}. Even at temperatures as high as T = 1.0 MeV there is still a recognizable low momentum absorption peak (with height well over 1), but it does not always correspond to a usual pasta (gnocchi, spaghetti or lasagna) in our simulations. At such high temperatures and for most densities, the system is in a “sponge-like” structure which is, nevertheless, ordered enough to produce a recognizable peak in S(k).

Similarly, Fig. 70 plots the height of the low momenta peak as a function of the temperature for the four structures of Fig. 58. These structures correspond to unusual intertwined lasagna and other irregular pasta which tend to appear when the cooling procedure drives the system at temperatures below T~0.7 MeV. In such conditions the system may collapse into several distinct structures, in addition to the usual lasagna which is the ground state at this density. Figure 70 shows the peaks of S(k) in the low-momentum region corresponding to the structures in Fig. 58. When we compare the absorption of each structure, we see that not only the usual lasagna has a peak on the low momenta region, but also unusual pasta phases: two that resemble an intertwined lasagna and another one that does not look like any other pasta (see captions for details).

Despite being different from the usual pasta phases, these shapes have a peak for low momentum in the structure factor. In Fig. 70 we see the corresponding absorption peaks for each structure. The unusual pasta phases show a larger absorption than the usual lasagna in the range of temperatures shown in the figure mentioned above.

4.5 Properties of non-traditional pasta

Usual pasta shapes are ground states (potential energy minima). The nontraditional structures described in the previous section are likely to be local potential energy minima, which abound in frustrated systems like this. The complexity of the energy landscape (many local minima separated by energy barriers) makes it difficult to reach the actual ground state by simple cooling in molecular dynamics simulations. However, since we are working at fixed number of particles, volume and temperature [(N,V,T) ensemble implemented through a thermostat, like Nosè–Hoover], the equilibrium state of the system at finite temperatures is not that which minimizes the internal energy but that which minimizes the Helmholtz free energy, A=E-TS. All of these structures may then be actual equilibrium solutions, as long as they are free energy minima.

An accurate calculation of free energies from MD simulations is computationally very expensive [58], specially at low temperatures when overcoming energy barriers become very improbable events. However, we can easily compute the internal energy distributions over a long evolution at constant temperature. In Fig. 71 we show internal energy histograms constructed from very long thermalized evolutions at T = 0.6 MeV using three of the systems shown in Fig. 58 as initial conditions. We see that, although the histograms clearly differ, they overlap significantly. This fact indicates that the full ensemble of equilibrium configurations at T = 0.6 MeV contains all of these structures, not only lasagna. In light of this we propose that at low but finite temperatures, the state of the system should be described as an ensemble of both traditional and nontraditional structures rather than by a single one.

When we heat up the system to T = 0.8 MeV, these three histograms become indistinguishable, hinting that, for this temperature, the free energy barriers can be surpassed, and the system is more likely to be ergodic.

These observations are relevant because all of these structures show peaks in S(k) at the same wavelenght (within the uncertainty), although of different heights. And more importantly, we find from our calculations that the seemingly amorphous, sponge-like structures can be more efficient in scattering neutrinos of the same momentum that any usual pasta (i.e., have higher peaks), usually invoked as a necessity for coherent neutrino scattering. This result shows that unusual pasta shapes should also be considered when studying the structure of a neutron star’s crust.

4.6 The nucleon thermal conductivity

Neutron stars are expected to reach temperatures as high as T~9 MeV, and be cooled soon after birth. The late-time cooling of neutron stars (e.g., MXB 1659–29) was shown to be consistent with low thermal conductivities [59, 60]. Thermal conductivity is believed to be due to electron flow, and not to energy flux in collisions between nuclear species [61–63].

As the pasta is expected to dominate the inner crust at sub-saturation densities [21, 45], the pasta structures may enhance or hinder the energy transport [64]. Models based on this scenario have attempted to explain the cooling of neutron stars, finding that the conductivity can vary with the alignment of the pasta, reducing it by 37% for randomly oriented pasta slabs [61].

In this section the thermal conductivity of NSM will be studied. We will not consider the energy transport due to electrons, regardless that a screening potential is always present, our concern is with the nucleons’ energy transport. In summary, the phononic thermal conductivity of nucleons embedded in an electron gas will be studied as it undergoes the transition to the pasta regime and to the solid–liquid phase transition. We first present the theoretical background for the thermal conductivity κ in the context of the CMD, followed by the measurements of κ within the pasta scenario divided into isospin symmetric and non-symmetric cases. The procedure used to calculate the thermal conductivity is described in Appendix F.

4.6.1 Procedure

The MD simulations were performed as specified in Appendix A.2.2. We restricted the study to elementary shapes, namely lasagnas or spaghettis, created with 4,000 nucleons at x = 0.5. The procedure was to cool the system from T = 4 MeV down to T = 0.1 MeV where a solid pasta state was formed. The density ranged from 0.03 fm⁻³ to 0.05 fm⁻³.

The lasagnas and spaghettis were produced and aligned with the cell coordinates by means of a transformation. For the thermal conductivity measurements, the temperature was increased from 0.1 MeV to 2.1 MeV, and the configurations were set as the initial conditions for the measurements; the data were collected after steady states were reached.

The nucleonic thermal conductivities were obtained following the Müller–Plathe procedure [65] which requires the binning of the primary cell to compute the temperature gradient across the bins; the number of bins was set to 20 across which a linear temperature profile was established. For each pasta shape the “parallel” (along the pasta) and “transverse” (across the pasta) thermal conductivities were computed. Furthermore, the thermal conductivity was calculated for protons and neutrons separately and for all the nucleons. Computing the conductivity for each specie means that only such specie contributes to the velocity exchange in Eq. (F.3).

4.6.2 κ of symmetric NSM

Figure 72(b) shows the nucleonic thermal conductivities of symmetric neutron star matter at ρ = 0.05 fm⁻³ obtained during heating from 0.1 to 2.1 MeV. Clearly seen are changes at the solid–liquid transition found before at T~0.5 MeV, and the pasta formation at T~1 MeV. The computation was performed for all the nucleons and for the protons only; the figure also includes data reported in literature [64].

Figure 72(b) shows that κ goes from a smooth curve above T\simeq 1.25\text{MeV} to a “decoupling” between the thermal conductivity parallel to the lasagna ({\kappa }_z) and the one orthogonal to this direction ({\kappa }_x). The “decoupling” pattern is essentially the same whether all the nucleons are considered or only the protons.

To understand the decoupling we measured the number of clusters within the simulation cell. The cut-off distance between neighbors belonging to the same cluster was set to r_c=4 fm, in order to get exactly three clusters at T = 0.1 MeV [see Fig. 73(a)]. Weak connections between slabs could be observed at an early stage of the breakdown (say, 1–1.25 MeV), as can be seen in Fig. 73(b). But the coupling of the parallel and transverse conductivities [see Fig. 72(b)] occurred when the connectivity between slabs was established; Fig. 73(c) illustrates this phenomenon.

The reduction of the thermal conductivity across a lasagna (say, for T<1 MeV) is due to the existence of voids between the slabs. The negative slope for the parallel κ (i.e., along the lasagna) is due to the larger conductivity in the solid pasta with respect to the liquid one. Figure 74 shows how these openings spread over the slabs until the lasagna becomes more of an spaghetti-like structure [see Fig. 74(a)]. Figure 75 exhibits the corresponding proton thermal conductivity.

In summary, cold pastas can only conduct heat along the pasta structure, but warming the pastas above T\simeq 1.25\text{MeV} connects regions that were separated allowing heat transfer on any direction, and thus thermal conductivity becomes an (isotropic) value, that may depend on the system density.

4.6.3 κ of non-symmetric NSM

The thermal conductivity for non-symmetric neutron star matter was computed in the same way as in the previous Section. Figure 76 shows the proton thermal conductivity behavior for systems with x = 0.3, evolving from “cold” (solid) temperatures to “warm” ones, along with the corresponding structures (protons only).

Figure 76 is for x = 0.3 as Fig. 72 is for the symmetric case, and they both show qualitatively similar trends. The proton conductivity across slabs vanishes, and the decoupling also occurs at around the same temperature.

The dependence of κ on x is shown in Fig. 77 for the proton thermal conductivity. It is not surprising that on the pasta phase (T\lesssim 1\text{MeV}) the proton κ along the lasagna increases with the x as more protons are available for transmitting heat. The proton κ appears to be independent of x except at the temperature at which the structures melt.

Figure 78 shows the thermal conductivities, {\kappa }_x, {\kappa }_y and {\kappa }_z, of all nucleons and of neutrons only for systems with x = 0.3. The curves are very similar due to the fact that at x = 0.3 most of the conduction is due to the neutrons.

4.6.4 Summary of findings

Studying pasta structures in the form of lasagnas and spaghettis made out of 4000 nucleons at x = 0.3, 0.4, and 0.5, it was determined that the nucleonic thermal conductivity changes whenever the structure of neutron star matter switches from pastas to more isotropic structures at around T ≈ 1 MeV as a fraction of nucleons begin to bridge the slabs and rods.

In the pasta phase the nucleonic thermal conductivity in directions parallel and transverse to the pasta axis separate at T\lesssim 1\text{MeV} due to the creation of void regions. This effect occurs both for isospin symmetric and non-symmetric NSM. The thermal conductivity for non-symmetric matter shares the same qualitative behavior as the symmetric matter, with the differences being explained by the structure of the pasta and the neutrons in excess at different values of x.

4.7 Summary of NSM properties

In this section we have studied neutron star matter, the formation of pastas, phase transitions, symmetry energy, structure function and nucleonic thermal conductivity. The study was performed for isospin symmetric and for asymmetric NSM.

This was done with molecular dynamics studies of systems with 4000 nucleons in a cell with replicas in periodic boundary conditions. The systems were cooled and heated between temperatures in the range of 0.2 MeV≤T≤4 MeV and densities from 0.02 fm⁻³ to 0.085 fm⁻³ for isospin content symmetric x = 0.1, 0.2, 0.3, 0.4 and 0.5.

The caloric curve presented changes in its slope at T ≈ 0.5 MeV and at T ≈ 1.5 MeV. These changes are identified as the onset of the pasta (1.5 MeV) and a liquid-to-solid phase transitions of the nucleons inside the pasta (0.5 MeV). These transitions were characterized with the Lindemann coefficient and with the Minkowski functionals as a function of the temperature. The T ≈ 0.5 MeV phase transition doesn’t alter the typical pasta shape (lasagna and spaghetti).

The morphology of the pastas at different conditions of T, ρ, and x was investigated by means of the Euler characteristic χ. For temperatures above T\simeq 1\text{MeV}, χ went from negative values (tunnel dominated configurations) to positive ones (configurations dominated by cavities) at around ρ = 0.05 fm⁻³.

The structure function (and the neutrino absorption) are also connected to the morphology of the pasta. The maxima of the structure function maintained high values at low temperatures (0\lesssim T\lesssim 0.7\text{MeV}), as the pasta shapes are relatively well ordered in this temperature range.

The symmetry energy and the nucleonic thermal conductivity also presented changes as the phase changes. For temperatures above 2 MeV, E_{Sym} remained basically constant, but it jumped in magnitude at temperatures around which the slope of the caloric curves changed. On the other hand, κ showed two distinct regions above and below T ≈ 1.5 MeV showing again a change of the energy transport due to a change of the pasta morphology.

5 Conclusion

In this review we studied the nuclear pastas with the ultimate goal of understanding properties of neutron stars. With masses between 1 and 3 solar masses and a radius of about 10 km, neutron stars have a crust of about 1 km where β decays yield a neutron-rich environment immersed in a sea of electrons. Densities of neutron stars range from normal nuclear density down to practically zero in the neutron star envelope, and initial temperatures of up to 9 MeV cooling down rapidly by thermal radiation and neutrino emission. Nuclear systems at such conditions form non-uniform structures, known as pastas, which influence the cooling of the star by modifying the neutrino-crust scattering.

Up until recently it was believed that pastas were due solely to the interaction between nuclear and Coulomb forces which constitute a frustrated system, but it is now known that pastas can also form due entirely to the interplay of the attractive-repulsive nuclear forces. We divided the complete study of pastas into one focused on those formed in nuclear matter, in Section 2, followed by a study of the role of the electron gas on the formation of pastas, in Section 3, and finally one on the pastas created in neutron star matter, in Section 4.

5.1 Nuclear matter

The properties of nuclear matter, i.e., in a system of neutrons and protons without a gas of electrons, were studied in Section 2. Because of the different properties NM has at different temperatures, the study was divided into intermediate and low temperatures.

5.1.1 NM at intermediate temperatures

In Section 2.1 the properties of NM at intermediate temperatures (1 MeV \lesssim T\lesssim 15\text{MeV}) were studied. Using molecular dynamics with systems with a ratio of protons of x = Z/A = 0.3, 0.4 and 0.5, temperatures from T = 1 to 15 MeV, and \rho \lesssim 0.2 fm⁻³, the energy per nucleon, pressure saturation densities, compressibility, phases, phase diagram and symmetry energy were investigated.

The density dependence of the internal energy showed characteristic ∪ shapes around saturation densities, which signaled the existence of regions where NM is bound, unbound, and in homogeneous (crystal) or non-homogeneous structures (liquid–gas mixed phase). Although the general shape of the E–ρ curves was maintained for different values of x, the saturation densities and the compressibility varied. The radial distribution function and the mass distribution helped confirm the existence of a liquid phase and a liquid–gas mixture density at sub-saturation densities. The region were such phases exists, i.e., the phase diagram of NM, was determined from the pressure–density isotherms by means of the Maxwell construction; a byproduct was the determination of the critical points for different values of x.

The symmetry energy of NM was obtained at sub-saturation densities and compared to experimental data and other theories. The CMD-based symmetry energies was able to explain the non-zero values of E_{Sym} in the limit of low density, as demanded by experiments.

5.1.2 NM at low temperatures

In Section 2.2 the properties of NM at low temperatures (T\lesssim 1\text{MeV}) were studied. At such temperatures NM was found to produce crystal-like structures around saturation densities, and pasta-like structures at sub-saturation densities. Properties of these phases were studied using molecular dynamics with systems with a ratio of protons from x = Z/A = 0.1, to 0.5, temperatures from T = 0.01 to 4 MeV, and \rho \lesssim 0.2 fm⁻³.

The energy per nucleon as a function of the density also indicated the existence of non-homogeneous structures at lower densities. E(ρ,T) showed three distinct behaviors: pasta for ρ<0.085 fm⁻³, the crystalline structure for ρ>0.14 fm⁻³, and an intermediate region in between the first two. NM at zero temperature develops a simple-cubic crystalline structure at densities around saturation nuclear density. The zero-temperature binding energy found was −16 MeV at a density of 0.16 fm⁻³, and with a compressibility of 283 MeV. At higher temperatures (0.001 MeV≤T≤1.0 MeV), the ∪ shape is maintained still showing crystalline phases at densities \rho \gtrsim 0.13 fm⁻³.

Spherical bubbles, gnocchi, tunnels, lasagnas, and perfect crystals were found at ρ<0.085 fm⁻³. The various shapes were studied with cluster-recognition algorithms, caloric curve, the radial distribution function, the Lindemann coefficient, Kolmogorov statistics, Minkowski functionals.

The caloric curve had a change of slope at T = 0.5 MeV and a less-conspicuous one at T ≈ 2.0 MeV; the changes were observed at all values of x. The Lindemann coefficient, the radial correlation function and p–T curves helped identify the slope-changes as a solid–liquid phase transition within the pasta regime at the lower temperature, and as the pasta onset at the larger one. These results from the caloric curves, g(r), Lindemann coefficient, and the p–T curves, together indicate that the nucleons inside the pasta structures exist in different phases; these results were obtained both in isospin symmetric and non-symmetric nuclear matter.

The onset of the creation of the pasta was further studied with the Kolmogorov statistics. Major discrepancies to homogeneous distributions were detected at T ≈ 2 MeV indicating an early stage of the pasta formation. The changes in pasta morphology were detected at lower temperatures as bubbles widens into tunnels, and into slabs. Apparently, the excess of neutrons inhibits the pasta formation until lower temperatures where protons manage to form the pasta. The released neutrons get distributed along the cell disrupting the pasta structure. Likewise, the Euler characteristic χ helped to understand the changes in morphology. Unfortunately, the χ-curvature classification that worked well for symmetric systems becomes meaningless for the non-symmetric case. In spite of this, a general observation is that the system departs from homogeneity at T~1.5 MeV.

Finally, the symmetry energy of nuclear matter was studied in the pasta region. The E_{sym}, computed through a numerical procedure, showed a connection between the symmetry energy and the morphology of the system. A low value of E_{sym} at T>2 MeV, a larger value in the region 1.5 MeV<T<2 MeV while the liquid pasta is formed, an intermediate value in the range of 0.5 MeV<T<1.5 MeV where the liquid pasta exists, and the highest value at T<0.5 MeV when the liquid-to-solid phase transition happens within the pasta.

In summary, NM at intermediate temperatures can exist in liquid form, gaseous form, and in a mixture between these phases, whereas at lower temperatures crystals develop around normal saturation nuclear density and the pasta is created at sub-saturation densities. All of the properties of the medium reflect these different phases.

5.2 The electron gas

The role the Coulomb interaction strength and screening length of the electron gas have on the formation of the pastas in NM was studied in Section 3. An important observation is the existence of the pasta even without the presence of the electron gas, due to the competition between the attractive and the repulsive nuclear interactions. The strength of the Coulomb interaction emphasizes the pasta structures, makes less-compact objects, decreases the x content of fragments, and increases of nucleon mobility.

Our study determined that there is a minimal screening length of the Coulomb interaction of the electron gas for the CMD to produce structures independent of the cell size. For the Pandharipande potential {\lambda }_c lies between 10 fm and 15 fm depending on the density.

5.3 The pasta in neutron star matter

In Section 4 we studied the formation of pastas in neutron star matter. In particular, phases, phase transitions, symmetry energy, structure function and thermal conductivity were investigated for isospin symmetric and asymmetric NSM. Using molecular dynamics systems of 4000 nucleons per cell in periodic boundary conditions were cooled and heated between 0.2 MeV≤T≤4 MeV, densities from 0.02 fm⁻³ to 0.085 fm⁻³ for isospin content symmetric from x = 0.1 to 0.5.

Again, using the caloric curve, the Lindemann coefficient and the Minkowski functionals, the pasta was found to form at around T ≈ 1.5 MeV, and a liquid-to-solid phase change was detected on the nucleons inside the pasta at T ≈ 0.5 MeV without altering the pasta shape. Similarly, the Euler characteristic χ helped to determine that for T\simeq 1\text{MeV}, the pasta went from being tunnel-dominated to becoming configurations dominated by cavities.

The structure function, related to neutrino absorption, was also found to be dependent on the morphology of the pasta. with the largest neutrino-crust scattering taking place at 0\lesssim T\lesssim 0.7\text{MeV}, where the pastas are well ordered. Likewise, the symmetry energy and the thermal conductivity present values according to the pasta structures. In particular, κ had only two main values, above and below T ≈ 1.5 MeV.

6 Appendices

6.1 Classical molecular dynamics

In this work we use the Classical Molecular Dynamics (CMD) model, which represents nucleons as classical particles interacting through pair potentials and calculates their dynamics by solving their equations of motion numerically. A definite advantage of CMD over other methods is that its dynamics includes all particle correlations at all levels, i.e., 2-body, 3-body, etc. Indeed the method can describe nuclear systems ranging from highly correlated cold nuclei (such as two approaching heavy ions in their ground state), to hot and dense nuclear matter (nuclei fused into an excited blob), to phase transitions (fragment and light particle production), to hydrodynamics flow (after-breakup expansion) and secondary decays (nucleon and light particle emission).

The greatest disadvantage of CMD is the lack of quantum effects, the effect of this deficiency on the formation of pasta, however, is ameliorated by a number of factors. In the region of sub-saturation density and intermediate temperatures the density of available energy states grows much faster that the number of particles per unit volume and thus renders Pauli exclusion principle useless [32], it is not until lower temperatures and higher densities that Pauli blocking has an effect in dictating the nuclear dynamics. For the case of the pasta, at T<1 MeV and at sub-saturation densities, the dynamics is much limited due to the freezing in place of the nucleons. Indeed Pauli blocking at low excitation energies may be disregarded in the context of the pasta formation environment as momentum-transferring collisions cease to be important. Other quantum effects, such as superfluidity and superconductivity, presumed to exist in neutron stars, are not expected to interfere in the macroscopic structure of the pasta, nor on the bulk properties here studied.

For the nuclear case CMD is used with the Pandharipande potentials, which were designed by the Urbana group to reproduce experimental cross sections in nucleon–nucleon collisions of up to 600 MeV [42]. Such potential mimics infinite systems with realistic binding energy, density and compressibility and to produce heavy-ion dynamics comparable to those predicted by the Vlasov–Nordheim equation. This parameter-free model has been successfully used to study nuclear reactions obtaining mass multiplicities, momenta, excitation energies, secondary decay yields, critical phenomena and isoscaling behavior that have been compared to experimental data [26, 66–74]. More recently, and of interest to the present work, the model was used to study infinite nuclear systems at low temperatures [41] and in neutron star crust environments, including the pasta structures that form in NM and NSM [10, 18, 19, 21, 32, 51].

6.1.1 Potentials

Neutron star matter is composed of protons, neutrons, and electrons, CMD uses pp, nn, and np potentials, as well as an screening potential to mimic the effect of the electron gas; these potentials are now described in turn.

6.1.1.1 Nucleon–nucleon potentials

In the CMD model nucleons interact through the Pandharipande (Medium) potentials. These potentials attain a binding energy E({\rho }_0)=-16 MeV/nucleon and a compressibility of about 250 MeV. The corresponding mathematical expressions are

\begin{array}{}\text{(A.1)}& \begin{array}{l}{V}_{np}(r)=\frac{V_r}{r}{\text{e}}^{-{\mu }_{r}r}\text{}\text{}-\text{}\text{}\frac{V_r}{r_c}{\text{e}}^{-{\mu }_r{r}_c}\text{}\text{}-\text{}\text{}\frac{V_a}{r}{\text{e}}^{-{\mu }_{a}r}\text{}\text{}+\text{}\text{}\frac{V_a}{r_c}{\text{e}}^{-{\mu }_a{r}_c}\mathrm{,}\\ V_{nn}(r)=\frac{V_0}{r}{\text{e}}^{-{\mu }_{0}r}-\frac{V_0}{r_c}{\text{e}}^{-{\mu }_0{r}_c}\mathrm{,}\end{array}\end{array}

where r_c is the cutoff radius after which the potentials are set to zero. Although the parameters μ_r, μ_a, μ₀ and V_r, V_a, V₀ were first set by Pandharipande for cold nuclear matter [42], a recent improvement [45], here named New Medium, reproduces the cold nuclear matter binding energies more accurately and, thus, is used in this work. The corresponding values are summarized in Table 2. Figures 79 and 80 contrast these potentials with those of Pandharipande Medium potentials.

It is worth mentioning that these potentials were crafted by Pandharipande [42] to reproduce experimental cross sections in nucleon–nucleon collisions without the use of an explicit Coulomb potential for the proton–proton interactions; the Coulomb interactions among protons can be assumed to be embedded in the Pandharipande potentials. Similarly, as the potentials do not include an explicit momentum dependence, it is not possible to extract the momentum dependence of the properties here studied, such as compressibility, saturation density, symmetry energy, etc.

6.1.1.2 The Coulomb potential

As in neutron star crusts electrons filled the space between protons and neutrons, it is necessary to include the Coulomb effect of an electron gas. The main effect of the electron cloud is to introduce an screening effect on the Coulomb potential of the protons. Given the infinite range of Coulomb interaction some approximation is needed when simulating infinite systems.

In the CMD model protons and neutrons are then assumed to be immersed in a uniform gas of non-interacting electrons, and the effect of the electron gas is implemented either by the Thomas–Fermi screening method, or by the Ewald Summations. In the present work we use first one and refer the reader to [18, 58, 75] for details on the second method.

The implementation of the Thomas–Fermi Screening considers the electron gas as an ideal Fermi gas at the same number density as protons. This electron gas, being uniform, does not exert any force on protons but becomes polarized in their presence effectively “Screening” the proton’s charge. After solving the corresponding Poisson equation, the field generated by a screened proton takes the form:

\begin{array}{}\text{(A.2)}& V_C=\frac{q^2}{r}{\text{e}}^{-r/\lambda }\mathrm{,}\end{array}

where λ is a screening length. Thus, the effective Coulomb potential becomes finite ranged.

The relativistic Thomas–Fermi screening length is given by

\begin{array}{}\text{(A.2)}& \lambda =\frac{{\pi }^2}{2\text{e}}{\left[k_F\sqrt{(k_F^2+m_e^2)}\right]}^{-\frac{1}{2}}\mathrm{,}\end{array}

where m_e is the electron mass and k_F is the electron Fermi momentum defined as k_F={\left(3{\pi }^2{\rho }_e\right)}^{1/3}, where {\rho }_e is the electron gas number density (taken equal to that of the proton’s).

To avoid finite size effects, λ should be significantly smaller than the size L={\left(A/\rho \right)}^{\frac{1}{3}} of the simulation cell. Since the screening length λ depends on the density of the system, it is always possible to satisfy this condition by increasing the simulation box size along with the number of particles. However, this can lead to prohibitively large systems for our current computation capabilities. Following a prescription given in Ref. [15] we set λ = 10 fm, and set V_C^{(Scr)}=0 at a cutoff distance of 20 fm; these values are long enough to reproduce the density fluctuations in the cell size used and, in the low temperature case, they ensure that the properties of the resulting pasta remain essentially constant; see [76]. This implementation has been used before [3, 18, 20].

Figures 79 and 80 show the interaction potentials between nucleons without the existence of the surrounding electron gas, while Fig. 81 shows proton–proton complete potential (including the Coulomb screening). The exponential cut-off renders the Coulomb effective interaction short ranged. In this way the energy and the entropy is additive (i.e., energy scales with the number of particles).

Figure 82 shows an example of the pasta structures for nuclear matter with 6000 nucleons in the simulating cell (and periodic boundary conditions), x = 0.5 at T = 0.2 MeV, and densities ρ = 0.05, 0.06, 0.07, and 0.085 fm⁻³, respectively.

6.1.2 Implementation

6.1.2.1 Intermediate temperature

To study the bulk properties of nuclear matter systems must be created with specific values of density and temperature. For this purpose it is convenient to accelerate the reaching of equilibrium by adding a heat reservoir [77]. The addition of the heat flow variable to the classical equations of motion results in the Nosé–Hoover equations of motion which can be integrated by Størmer finite differences. In principle this approach corresponds to a canonical ensemble and does not conserve energy which is added or removed by the heat reservoir; configurations in thermal equilibrium, however, can be achieved faster than with the usual microcanonical formalism of newtonian mechanics and an Andersen’s thermostat [78].

To mimic an infinite system A = 2000 nucleons were placed in cubic cells under periodic boundary conditions. We focus on systems with isospin content of x = Z/A = 0.3, 0.4, and 0.5, where Z is the number of protons. The number densities were enforced by placing the nucleons in cubical boxes with sizes selected to adjust the density. The temperatures of the systems studied are T = 1, 2, 3, 4, and 5 MeV, and their densities were selected to be around and below the corresponding saturation densities values, which vary with isospin content and temperature.

The procedure followed is straightforward: the nucleons are placed at random within the cell avoiding overlaps (i.e., interparticle distances smaller than 0.1 fm) and endowed with a Maxwell–Boltzmann velocity distribution corresponding to the desired temperature. The system then is rapidly evolved until the temperature is maintained within 1%. After reaching thermal equilibrium, the system continues evolving and its information at selected time steps (nucleon positions and momenta, energy per nucleon, pressure, temperature, density, etc.) is stored for future analysis. Each data point represents the average of 200 thermodynamically independent configurations, the average of the standard deviations is 0.036 MeV.

For the radial potentials used the pressure and the energy per nucleon can be calculated through

where the angular brackets indicate an average over the appropriate ensemble, the number density is ρ = N/V, N is the number of nucleons, T is the temperature (in MeV), and r and F are the position and force vectors between pairs of nucleons, respectively, V_i is the potential energy of each configuration included in the ensemble average, and v_i and m_i are the velocity and masses of each of the nucleons. The first term of the pressure is the ideal gas contribution and the second one is the contribution arising from the inter-nucleon potentials.

6.1.2.2 Low temperature

To study infinite nuclear matter systems we resort to the LAMMPS code [101]. We study the properties of a system of 6000 particles (with periodic boundary conditions) using the Pandharipande and screened Coulomb potentials. The total number of particles is divided into protons (P) and neutrons (N) according to values of x = P/(N+P)=0.3, 0.4, and 0.5. The nuclear system is cooled down from a relatively high temperature (T≥4.0 MeV) to a desired cool temperature in small temperature steps (ΔT = 0.01 MeV) with the Nosé Hoover thermostat [79], and assuring that the energy, temperature, and their fluctuations are stable.

Figure 83 shows an example of the energy per nucleon versus the density for systems with 2000 particles at x = 0.5 and T = 1.5, 1.0 and 0.5 MeV. Clearly visible are the homogeneous phase (i.e., those under the “∪” part of the energy–density curve), and the loss of homogeneity at lower densities.

From Fig. 82 and Fig. 83 it is possible to distinguish three distinct regions to be analyzed. The first one, that goes from very low densities up to approximately 0.08 fm⁻³ in which the system displays pasta structures, a crystal-like region at densities above approximately 0.14 fm⁻³, and a transition region between these two [21].

6.1.2.3 CMD for neutron star matter

In order to mimic (asymmetric) neutron star matter, the LAMMPS CMD code [101] was fitted with the potentials mentioned in Appendices A and A.1.2. We tracked the evolution of systems with A = 4000 nucleons situated in a cubic cell under periodic boundary conditions. The simulated isospin content was x = z/A = 0.1, 0.2, 0.3, 0.4, and 0.5. The densities were varied between 0.02 fm⁻³ to 0.085 fm⁻³. The temperature was controlled with a Nosé–Hoover thermostat slowly varying from T = 4 MeV down to 0.2 MeV (ΔT<0.1%). After placing the nucleons at random, but with a minimum inter-particle distance of 0.01 fm, the nucleons were endowed with velocities according to a Maxwell–Boltzmann distribution to correspond to a desired temperature, and the equations of motion were solved to mimic the evolution of the system. The nucleon positions, momenta, energy per nucleon, pressure, temperature, and density, were stored at fixed time-steps.

6.2 The Maxwell construction

Within the coexistence region, the liquid and gaseous phases are in thermal and chemical equilibrium with one another and the energy needed to add or subtract a nucleon from either phase is zero. As such energy is the Gibbs free energy per particle, g=ϵ-Ts+pV, its infinitesimal changes are given by \text{d}g=-s\text{d}T+V\text{d}p and, for isothermic processes (where dT = 0), it reduces to dg = Vdp. The determination of the boundary of the coexistence region is then reduced to finding the points where g={\displaystyle \int }V\text{d}p=0, which can be done through an interesting geometrical technique known as the Maxwell construction.

[Notice that the present calculation does not consider the isospin chemical potential in its conditions for equilibrium. That is, the Maxwell construction calculated only satisfies that the two phases have the same baryon chemical potential, but not the same isospin chemical potential. This is a limitation of the model used and its effect on the results has yet to be bounded.]

To implement the Maxwell construction, the isothermal pressure–density curves must be first turned into pressure volume curves (simply through V = 1/ρ), and then inverted to yield isothermal V(p) curves; the edges of the coexistence region for such temperature will be the limits of the integral, V_{Gas} and V_{Liq}, at which g={\displaystyle {\int }_{V_{Liq}}^{V_{Gas}}}\text{}V\text{d}p=0.

To carry out this method one can construct a complete pressure–density isotherms by using the curves obtained from CMD [such as those in Fig. 2(a)], connected at very low densities with the p–ρ isotherms of a free nucleon gas by means of an interpolation.

At low densities, the pressure of a mildly-interacting gas of nucleons approaches that of a free Fermi gas, p(\rho \mathrm{,}T)=\frac{2}{3}\rho {\epsilon }_F(\rho \mathrm{,}T) which, using the parametrization of the Fermi energy introduced in Ref. [28], can be approximated by

\begin{array}{}\text{(B.1)}& p(\rho \mathrm{,}T)=\frac{\rho }{3}{\displaystyle \sum _{i=2}^5}i{a}_i{(\rho /{\rho }_0)}^{i/3+1}+\frac{2\rho }{3}{\displaystyle \sum _{i=0}^2}{\epsilon }_i(T){\rho }^i\mathrm{.}\end{array}

where {\epsilon }_i(T)={\displaystyle {\sum }_{j=1}^2}{\epsilon }_{ij}{T}^j, and the coefficients ε_ij and a_i are listed in Table 3.

The first term of Eq. (B.1) is a simple approximation to 2ρ/3 times the Fermi energy of a cold nucleon gas, and the second term is an approximation to 2ρ/3 times the T-dependent Fermi gas energy. Notice that, as expected at low densities, Eq. (B.1) does not depend on the isospin content x, as the system is highly noninteracting.

Using this approximation for the low density region (\rho \lesssim {\rho }_0/6), and the CMD pressure density curves for the liquid phase (from \rho \gtrsim {\rho }_L), it is possible to obtain a complete pressure–density curve by least-squares cubic interpolation constrained to match the three segments and the first derivatives of the curve at the two matching points. Figure 84 shows one example of a resulting pressure–density curve; the values used for the matching liquid density {\rho }_L vary depending on the values of T and x.

To implement the Maxwell construction, the pressure–density curves, p(ρ,T), must be inverted to yield V(p,T) as shown in Fig. 85 for the T = 9 MeV, x = 0.4 pressure isotherm. In this curve the gaseous phase is the region that includes the points from A to D, the unstable (negative compressibility) region goes from point D to point F, and the liquid phase is from F to I and beyond. The boundary of the coexistence region are the points C and G, such that the area to the right of the dotted line C–G and the curve C–D–E equals that to the left of the dotted line and the curve E–F–G, i.e., g={\displaystyle {\int }_{V_G}^{V_C}}\text{}V\text{d}p=0.

To determine the boundary points of the coexistence region the method outlined in the previous paragraph was implemented for all the cases studies with CMD. Figure 86 shows three pressure isotherms plotted versus the volume for the cases T = 1, 10, and 15 MeV and for x = 0.35. The continuous curves of the left show a fit to the CMD results (points), and the dashed curves of the right are the least square fits inverted to appear as a function of the volume. The coexistence region is obtained by repeatedly integrating the area between a given pressure value and the pressure isotherm, until the pressure at which the integral is zero is found. The resulting pressures for the cases T = 1 and 10 MeV are indicated with the horizontal lines; for the case of T = 15 MeV no pressure satisfied the condition of zero integral indicating that for x = 0.35 the coexistence region ends at a lower temperature. For further details see [80].

6.3 Nuclear symmetry energy from CMD at intermediate temperatures

The effect of the excess of neutrons to protons in the nuclear equation of state (EOS) is characterized by the symmetry energy, E_{sym}(T\mathrm{,}\rho ), and its importance in phenomena ranging from nuclear structure to astrophysical processes has prompted intense investigations [81, 82].

Some of the latest experimental and theoretical studies of the symmetry energy have been at subsaturation densities and warm temperatures [33, 34]. Experimental reactions [34–36], for example, have shown that E_{sym} is affected by the formation of clusters. A recent calculation of the symmetry energy at clustering densities and temperatures [33] obtained good agreement with experimental data, corroborating the Natowitz conjecture [34], namely that the asymptotic limit of E_{sym} would not tend to zero at small densities as predicted by mean-field theories.

Formally, the symmetry energy is the correction that must be added or subtracted to the binding energy in systems with different numbers of protons and neutrons. Defining the symmetry energy as

\begin{array}{}\text{(C.1)}& E_{Sym}(\rho \mathrm{,}T)=\frac{1}{2\mathrm{！}}{\left[{\partial }^{2}E(\rho \mathrm{,}T\mathrm{,}\alpha )/\partial {\alpha }^2\right]}_{\alpha =0}\mathrm{,}\end{array}

with \alpha =(N-Z)/(N+Z)=1-2x. The evaluation of the symmetry energy follows the procedure introduced before [32, 33].

It is possible to use the CMD results of the internal energy E(\rho \mathrm{,}T\mathrm{,}x) to construct a continuous function by means of fitting techniques, and use it to obtain the internal symmetry energy. This can be done by fitting the CMD results of E(T\mathrm{,}\rho \mathrm{,}\alpha ) for each T and x with an expression of the type

\begin{array}{}\text{(C.2)}& \begin{array}{ll}E(T\mathrm{,}\rho \mathrm{,}\alpha )=\hfill & E_0(T\mathrm{,}\alpha )+E_1(T\mathrm{,}\alpha )\rho \hfill \\ \hfill & +E_2(T\mathrm{,}\alpha ){\rho }^2+E_3(T\mathrm{,}\alpha ){\rho }^3\mathrm{.}\hfill \end{array}\end{array}

Figure 87 shows the behavior of the energy as a function of the density at T = 1 MeV and x = 0.5, 0.4, and 0.35 (or α = 0, 0.2 and 0.3, respectively). The liquid phase of the system is easily identifiable by its “∪” shape, as well as the liquid saturation density by the minima of the energy curves. Furthermore, the liquid–gas coexistence region can be identified as the lower density region where the CMD data separates from the “∪” curve; this is clearly noticeable for x = 0.4 and 0.5, but not so much for x = 0.35, which appears to be in the continuous phase up to very low densities, perhaps up to 0.02 fm⁻³.

The α dependence of the coefficients E_0(T\mathrm{,}\alpha ), E_1(T\mathrm{,}\alpha ), E_2(T\mathrm{,}\alpha ), E_3(T\mathrm{,}\alpha ) can be easily extracted from the CMD data assuming an α dependence of the type

\begin{array}{}\text{(C.3)}& E_i(T\mathrm{,}\alpha )=E_{i0}(T)+E_{i2}(T){\alpha }^2+E_{i4}(T){\alpha }^4\end{array}

for i = 0, 1, 2, and 3; odd terms in α are not included to respect the isospin symmetry of the strong force (without the Coulomb potential). The smooth behavior of these coefficients with respect to α can be inspected in Figs. 88 and 89; the curves on these figures are the results of least squares fits.

With this, the symmetry energy is then given by

\begin{array}{}\text{(C.4)}& \begin{array}{ll}{E}_{Sym}(T\mathrm{,}\rho )=\hfill & E_{02}(T)+E_{12}(T)\rho \hfill \\ \hfill & +E_{22}(T){\rho }^2+E_{32}(T){\rho }^3\mathrm{,}\hfill \end{array}\end{array}

with the coefficients E_{i2}(T) given in Table 4. Figure 9(a) shows the symmetry energy in the low density (i.e., in the liquid–gas mixture) region and its comparison to experimental data [35, 36].

Notice that the E_{Sym}(T\mathrm{,}\rho ) obtained from the CMD values of E(T\mathrm{,}\rho \mathrm{,}\alpha ) shows a smooth dependence on the density and temperature, and lie on the appropriate range of temperatures for the experimental data, which are known to vary between T = 3 to 11 MeV [34–36, 83]; thus corroborating the conjecture of Natowitz, et al. namely that E_{Sym}(T\mathrm{,}\rho ) should approach a constant value as the density approaches zero.

It is worth mentioning that the mechanism to extract E_{Sym} from CMD data is sensitive to the type of data used. For instance, the study of Ref. [32] was performed solely using data from the liquid phase, and produced results comparable to other theories at saturation densities, but incorrect results at low densities. Indeed, since such calculation did not use the E–ρ information about the liquid–gas mixed phase, it yielded the incorrect result of E_{Sym}\to 0 as \rho \to 0, as all other theories that do not include any clusterization; the validity of the results of [32] is thus limited to near saturation densities. On the other hand, the results presented here for the liquid–gas mixed phase contain clustarization data and, thus, yield the correct non-zero value for E_{Sym} as \rho \to 0 [33].

In summary, the applications of the method must be adapted to the type of data used. The fits used for the “∪” data of E(T\mathrm{,}\rho \mathrm{,}\alpha ) in the liquid phase, for instance, should be different than those used for the data of the liquid–gas mixed phase; the latter needing an extra term in the ρ-expansion, Eq. (C.2). As we will see in Appendix D, the computation of E_{Sym} at low temperatures will require further adaption for the method to work with frozen nuclear matter.

6.3.1 Comparison to other theories

Comparison of symmetry energies obtained with different models is difficult due to a series of factors. As explained in detail elsewhere [84], there are complications in the various methodologies used both in experimental and theoretical works to extract E_{Sym}. This becomes more pronounced when dealing with systems in the liquid phase as opposed to those in the liquid–gas coexistence region, and in dealing with warm systems instead of cold ones. The seriousness of this problem was quantified in the comparison of different theories through the symmetry energy parameters, which managed to approve but a handful of cases out of hundreds of models [85, 86]. We start this section with a brief outline of some of these complications.

To begin with, there is no unique procedure used to obtain the symmetry energy. Although most experimental studies use the method of finite differences instead of that of derivatives, as Eq. (C.1), there are some others that use isoscaling [83]. These methods have different ranges of validity, finite differences, for instance, assumes a quadratic behavior for the dependence on α thus limiting its validity to near symmetric homogeneous matter, and excluding cases with larger isotopic asymmetry. Indeed substantial differences in symmetry energies calculated with the derivative definition and with finite differences have been reported in recent studies [84].

Secondly, these methodological differences become more pronounced at subsaturation densities. At low densities phase transitions produce coexisting regions of low and high densities which can have significant effects on the symmetry energy. This is important as it is known that the liquid and gaseous phases have different isotopic asymmetry: the high-density phase is more isospin symmetric and the low-density phase more asymmetric; a phenomenon called isospin distillation or fractionation [87]. [In fact the present results (finite values of E_{Sym} at low densities) are believed to be due more to the separation of phases than to the process of cluster formation which is not taken into account in our calculations in thermal equilibrium [84].]

Another source of discrepancies surfaces when the concept of symmetry energy is applied to finite systems (as in experiments) or to infinite matter, as well as at zero or finite temperature. Experimental constraints obtained from laboratory experiments, which should apply to finite cold-to-warm systems, are fundamentally different than those coming from stellar systems which obey considerably distinct equations of state.

Yet another complication in comparing finite temperature theories to zero temperature calculations is that their thermodynamics are different. At finite temperatures, and depending on the thermodynamic conditions, the appropriate quantity that needs to be studied may be the free symmetry energy instead of the internal symmetry energy; although the two energies are equal to each other at T = 0 MeV, they can be substantially different at large temperatures, as exemplified in a recent calculation [84].

In summary, the symmetry energy is sensitive both to the methodology used for its calculation, and to the specific thermodynamic conditions of the system. A comparison of symmetry energies from different theories and experiments demands a careful consideration of possible discrepancies, especially if comparing cluster-forming calculations with theories that are unable to undergo phase transitions.

6.4 Nuclear symmetry energy from CMD at low temperatures

The fitting procedure outlined in Section 11 corresponds to a two step-fitting method. The coefficients E_i (i=\mathrm{0,}\cdots \mathrm{,3}) from Eq. (C.2) are computed at the first stage of the procedure. These coefficients, however, depend on the fraction α, although not directly on the fraction x (recall that α = 1−2x). A more suitable (low order) expression for the E_i’s instead of Eq. (C.2) is

\begin{array}{}\text{(D.1)}& \begin{array}{ll}{E}_i(T\mathrm{,}x)\hfill & =E_{i0}+E_{i1}\stackrel{\alpha }{\overbrace{(1-2x)}}+E_{i2}\stackrel{{\alpha }^2}{\overbrace{(1-4x+4x^2)}}+\mathcal{O}({\alpha }^4)\hfill \\ \hfill & \simeq (E_{i0}+E_{i1}+E_{i2})-2(E_{i1}+2E_{i2})x+4E_{i2}{x}^2\mathrm{.}\hfill \end{array}\end{array}

The above expression corresponds to the lowest non-trivial expansion E_i(T\mathrm{,}x)\simeq {\tilde{E}}_{i0}+{\tilde{E}}_{i1}x+{\tilde{E}}_{i2}{x}^2. Thus, the {\tilde{E}}_{ij}\text{’s}(j = 0,1,2) are related to the E_{i0}, E_{i2} coefficients according to the following matrix relation:

\begin{array}{}\text{(D.2)}& \left(\begin{array}{rrr}\hfill 1& \hfill 1& \hfill 1\\ \hfill & \hfill & \hfill \\ \hfill 0& \hfill -2& \hfill -4\\ \hfill & \hfill & \hfill \\ \hfill 0& \hfill 0& \hfill 4\\ \hfill & \hfill & \hfill \end{array}\right)\left(\begin{array}{c}{E}_{i0}\\ \\ E_{i1}\\ \\ E_{i2}\\ \end{array}\right)=\left(\begin{array}{c}{\tilde{E}}_{i0}\\ \\ {\tilde{E}}_{i1}\\ \\ {\tilde{E}}_{i2}\\ \end{array}\right)\mathrm{.}\end{array}

Notice that {\tilde{E}}_{i0}, {\tilde{E}}_{i1}, {\tilde{E}}_{i2} actually correspond to the fitting parameters at the second stage of the procedure outlined in Section C (that is, after the E_i’s were obtained). The E_{i0}, E_{i1} and E_{i2} parameters, though, are now computed indirectly from the linear system (D.2). It is straightforward that E_{i2}={\tilde{E}}_{i2}/4 at this (low order) approach.

The linear system (D.2) attains a single solutions for the three unknowns E_{i0}, E_{i1} and E_{i2}. However, any additional constrain (say, E_{i1}=0 for nuclear matter) would drive the system to an overdetermined state. The E_{i2} values might still be computed (or estimated) in that case, although E_{i0} or E_{i1} would become meaningless.

Our numerical computations for E_{i2} at the examined temperatures and densities indicated in Section 4.3 are detailed in Table 5.

In closing, formally speaking the symmetry energy is approximated from a second order expansion of the energy with respect to the parameter α, and evaluated at α = 0 and at saturation density. Our E_{sym}, however, is a generalization of such expansion that makes it dependent on the density and temperature. Since the use of the liquid drop mass formula that inspired such approximation of the symmetry energy is limited to finite systems at T = 0 and near saturation density, the E_{sym} here obtained cannot be simply inserted into such formula. A more complete equation of state of nuclear matter, i.e., one for infinite systems dependent on density, temperature and isospin content, could benefit from the analytical fit obtained here for non-zero temperature infinite systems at sub-saturation densities.

6.5 Symmetry energy for neutron star matter

The evaluation of the symmetry energy for neutron star matter follows the procedure introduced in Appendices C and D. The symmetry energy is defined as in Eq. (C.1), again using an E(T\mathrm{,}\rho \mathrm{,}\alpha ) as given by Eq. (C.2), but the α dependence modified as follows.

The α dependence of the coefficients E_i(T,α) can be extracted from the CMD data calculated at various values of α, and assuming an α dependence of the type

\begin{array}{}\text{(E.1)}& \begin{array}{ll}{E}_i(T\mathrm{,}\alpha )=\hfill & E_{i0}(T)+E_{i1}(T)\alpha +E_{i2}(T){\alpha }^2\hfill \\ \hfill & +E_{i3}(T){\alpha }^3+E_{i4}(T){\alpha }^4\mathrm{.}\hfill \end{array}\end{array}

Notice that no constrains on the powers of α are introduced in Eq. (E.1), as opposed for nuclear matter. This is a major difference between NM and NSM, since the latter disregards the isospin symmetry of the strong force due to the presence the Coulomb screening potential. The symmetry energy is then given by

\begin{array}{}\text{(E.2)}& \begin{array}{ll}{E}_{sym}(T\mathrm{,}\rho )=\hfill & E_{02}(T)+E_{12}(T)\rho \hfill \\ \hfill & +E_{22}(T){\rho }^2+E_{32}(T){\rho }^3\hfill \end{array}\end{array}

with the coefficients E_ij(T) obtained from the fit of the CMD data.

6.6 Thermal conductivity

The thermal conductivity κ relates the energy flux J to the temperature gradient \nabla T, through the Fourier law

\begin{array}{}\text{(F.1)}& \mathit{J}(t)=-\kappa \nabla T\mathrm{.}\end{array}

The energy flux J is obtained from the mean flux density \langle j\rangle transported across a small volume \mathcal{V}, that is \mathit{J}=\langle j\rangle \mathrm{.}\mathcal{V}. The calculation of thermal transport properties from atomistic simulations is well established [65, 88–90]. For the case of J and \nabla T being collinear (say, for example, along the \widehat{z} axis) a non-equilibrium method for computing the thermal conductivity {\kappa }_z is proposed by Müller–Plathe [65] from the average heat flux and temperature gradient.

\begin{array}{}\text{(F.2)}& {\kappa }_z=-\underset{t\to \infty }{{\displaystyle \text{lim}}}\frac{\langle J_z\rangle }{\langle \partial T/\partial z\rangle }\mathrm{.}\end{array}

If the medium is isotropic, common practice sets the mean thermal conductivity as ({\kappa }_x+{\kappa }_y+{\kappa }_z)/3. Notice that the linear nature of Eqs. (F.1) and (F.2) requires relatively small temperature gradients.

In a nutshell, the Müller–Plathe procedure [65] generates a heat flux of known magnitude and the temperature gradient is obtained as local averages of the kinetic energy. The system is divided in thin bins along the heat flux direction (see Fig. 90 for details); the first slab is labeled as the “cold” slab, while the slab in the middle is labeled as “hot”. A heat flux is generated by exchanging the velocities of two particles (with the same mass), the hottest particle in the “cold” bin and the coldest one in the “hot” bin (see Fig. 90). Thus, the system is (artificially) driven out of equilibrium, and a heat flux J develops through the system of interest in the opposite direction for the equilibrium restoration. This flux is expected to reach the stationary state if the exchanging rate is held regularly for a long time.

In order to generate an external heat flux, particle velocity exchanges are performed periodically during the MD simulation. Recall that the species themselves are not exchanged, but only the velocities. Thus, the “pumping” process only transports kinetic energy (for particles with the same mass). This procedure conserves total energy and linear momentum.

The heat flux introduced by the velocity exchange is hard to compute from dynamical magnitudes. The computation from the net transported (kinetic) energy is somewhat easier since

\begin{array}{}\text{(F.3)}& \langle J_z\rangle +\frac{1}{2A_{xy}}\left[\frac{1}{\tau }{\displaystyle \sum _{n=1}^T}\frac{1}{2}m(v_h^2-v_c^2)\right]=\mathrm{0,}\end{array}

where the expression between the square brackets represents the mean (kinetic) energy exchanged during the time period τ. v_h and v_c refers to the velocities of the hot particle and cold particle, respectively. The factor 2A_xy corresponds to the cross section of the slabs (two faces).

The temperature profile is obtained by computing the local (kinetic) temperature for each slab. Once the steady state is reached, the temperature profile is expected to be linear away from the cold and hot bins where velocities are exchanged, provided the heat flux remains small. Further details can be found in Ref. [65].

We stress the fact that the balance condition (F.3) links the heat flux J to the (artificial) kinetic energy transportation introduced by the Müller–Plathe procedure. The velocity exchange is not restricted to pairs of similar particles, but also across species (with the same mass). Therefore, the procedure enables the computation of the thermal conductivity for the set of all the nucleons, or for the set of protons and neutrons separately. The meaning of either coefficients, though, will be quite different.

6.7 Analytical tools

6.7.1 The radial distribution function

A rigorous definition for the radial distribution function g(r) starts from the following distance distribution

\begin{array}{}\text{(G.1)}& g(\mathit{r})=\frac{1}{{\rho }_0}\langle \frac{1}{N}{\displaystyle \sum _{i=1}^N}{\displaystyle \sum _{j\ne i}^N}\delta (\mathit{r}-{\mathit{r}}_{ij})\rangle \mathrm{,}\end{array}

where ρ₀=N/V is the (mean) density in the simulation cell of volume V (or equivalently L³). r_ij=r_j−r_i is the distance vector between the particle i and the particle j. The δ(·) function corresponds to the Dirac delta. The mean value indicated by 〈·〉 corresponds to the average operation over successive time-steps.

6.7.1.1 g(r) for intermediate temperatures

For cases where the existence of phases is suspected, it is instructive to compare the local density to the average one, g(r)=ρ(r)/ρ; peaks in g(r) will indicate the existence (or lack of) neighbors, thus signaling transitions from one phase to another, say from liquid to gas.

For intermediate temperatures, say T>1 MeV, Eq. (G.1) is applied cumulatively to a large number of systems (say, 200) at the same temperature and densities. The strengths of the nearest-neighbor peaks will indicate higher correlation at higher densities (continuous liquid phase) than at lower densities (gaseous phase composed of a mixture of free nucleons and small clusters). The difference in correlation is due to the larger nucleon mobility in the liquid–gas mixture compared to the homogeneous liquid phase. Similar effect is also observed in the second-neighbor peaks which show stronger correlations in homogeneous phases than in mixed ones. Figure 3 shows the radial distribution functions of a system at T = 1 MeV.

6.7.1.2 g(r) for pasta stuctures

To calculate g(r) for a pasta structure, we assume that particles form a homogeneous slab (that is, a lasagna-like structure) bounded between the \pm z_0 (horizontal) planes. The slab is quasi static, meaning that no averaging over successive time-steps is required.

We first consider a small region Ω(r) with volume δx δy δz in the r domain. Thus, we may evaluate the Eq. (G.1) as follows

\begin{array}{}\text{(G.2)}& {\displaystyle {\int }_{\mathrm{\Omega }(\mathit{r})}}g(\mathit{r}){\text{d}}^{3}x=\frac{1}{{\rho }_0}\frac{1}{N}{\displaystyle \sum _{i=1}^N}{\displaystyle \sum _{j\ne i}^N}{\displaystyle {\int }_{\mathrm{\Omega }(\mathit{r})}}\delta (\mathit{r}-{\mathit{r}}_{ij}){\text{d}}^{3}x\mathrm{,}\end{array}

where the integral on the right-hand side equals unity at the positions r = r_ij within Ω(r).

We now proceed to evaluate the sum over the j≠i neighbors of the particle i. The counting of neighbors is proportional to δx δy δz since the volume Ω(r) is very small. Thus, the tally is

\begin{array}{}\text{(G.3)}& {\displaystyle \sum _{j\ne i}^N}{\displaystyle {\int }_{\Omega (\mathit{r})}}\delta (\mathit{r}+{\mathit{r}}_i-{\mathit{r}}_j){\text{d}}^{3}x\approx \rho (\mathit{r}+{\mathit{r}}_i)\delta x\delta y\delta z\mathrm{.}\end{array}

The magnitude ρ(r+r_i) refers to the density within the small region Ω(r+r_i), its value remains constant inside the slab and vanishes outside. In fact, ρ can be expressed as {\rho }_s\mathrm{\Theta }(z_0-|z+z_i|), with \mathrm{\Theta }(\cdot ) representing the Heaviside function, and {\rho }_s being the slab density ({\rho }_s>{\rho }_0).

The sum over the i particles is evaluated through the integral of the infinitesimal pieces ρ_s d³x＇. Thus,

\begin{array}{}\text{(G.4)}& {\displaystyle \sum _{i=1}^N}\rho (\mathit{r}+{\mathit{r}}_i)\approx {\displaystyle {\int }_{\text{slab}}}\text{}\text{}{\rho }_s^2\mathrm{\Theta }(z_0-|z+z^{\prime }|){\text{d}}^3{x}^{\prime }\mathrm{,}\end{array}

where the slab domain corresponds to the bounded region \left|z^{\prime }\right|\le z_0. Notice that this integral vanishes for \left|z\right|\ge 2z_0 and equals {\rho }_s^2{L}^2(2z_0-|z|) otherwise. The g(r) function is then

assuming that g(r) is somewhat fixed inside the small domain Ω(r) in (G.2). We further replaced N by {\rho }_{s}2z_0{L}^2 in (G.2).

Eq. (G.5) is not exactly the radial distribution function yet because of the angular dependency of g(r). This dependency may be eliminated by integrating g(r) over a spherical shell of radius r, which introduces the normalization factor 4πr². The expression for the radial distribution then reads

\begin{array}{}\text{(G.6)}& g(\mathit{r})=\frac{1}{4\pi r^2}{\displaystyle {\int }_S}\frac{{\rho }_s}{{\rho }_0}(1-\frac{r\left|\text{cos}\theta \right|}{2z_0})r^2\text{sin}\theta \text{d}\theta \text{d}\phi \mathrm{,}\end{array}

where z = r cosθ. Notice that this expression is valid along the interval whenever , but it is constrained to the natural bound \left|\cos \theta \right|\le 1 if 2z_0\ge r. The integration of Eq. (G.6) finally yields

Notice that as r\to 0, g(r) correctly goes to the limit of g(r)→ρ_s/ρ₀. Similarly, for larger values of r, and up to , g(r) decreases linearly, as observed in Fig. 18 for different values of {\rho }_0.

Eq. (G.7) also indicates that g(2z_0)\approx 1 in the case that {\rho }_s/{\rho }_0\approx 2, in agreement with Fig. 18 where g(r) tends to unity at 20 fm (for simulation cells of L~45 fm); this figure, however, does not show the behavior beyond 20 fm as the statistics that can be collected for distances above L/2 are very poor.

From literature references, g(r) is supposed to converge to unity as r\to \infty. But, according to (G.7), the radial distribution vanishes at this limit. The disagreement corresponds to the fact that the condition g(\infty )=1 is only valid for homogeneous systems. The expression (G.7) can actually meet this condition if the slab occupies all the simulation cell, since {\rho }_s\to {\rho }_0 and z_0\to \infty (for periodic boundary conditions).

6.7.2 Lindemann coefficient

The Lindemann coefficient [91] provides an estimation of the root mean square displacement of the particles respect to their equilibrium position in a crystal state, and it serves as an indicator of the phase where the particles are in, as well as of transitions from one phase to another. Formally,

\begin{array}{}\text{(G.8)}& {\Delta }_L=\frac{1}{a}\sqrt{{\displaystyle \sum _{i=1}^N}\langle \frac{\Delta r_i^2}{N}\rangle }\mathrm{,}\end{array}

where \Delta r_i^2=(r_i-\langle r_i\rangle {)}^2, N is the number of particles, and a is the crystal lattice constant; for the nuclear case we use the volume per particle to set the length scale through a={(V/N)}^{1/3}.

6.7.3 Kolmogorov statistic

The Kolmogorov statistic measures the difference between a sampled (cumulative) distribution F_n and a theoretical distribution F. The statistic, as defined by Kolmogorov [92], applies to univariate distributions (1D) as follows:

\begin{array}{}\text{(G.9)}& D_N={\text{sup}}_{\left\{x\right\}}|F_N(x)-F(x)|\mathrm{,}\end{array}

where “sup” means the supremum value of the argument along the x domain, and N is the total number of samples. This definition is proven to represent univocally the greatest absolute discrepancy between both distributions.

An extension of the Kolmogorov statistic to multivariate distributions, however, is not straight forward and researchers moved in different directions for introducing an achievable statistic [93]. The Franceschini’s version seems to be “well-behaved enough” to ensure that the computed supremum varies in the same way as the “true” supremum. It also appears to be sufficiently distribution-free for practical purposes [94].

The three dimensional (3D) Franceschini’s extension of the Kolmogorov statistic computes the supremum for the octants

for any sample (X_i\mathrm{,}{Y}_i\mathrm{,}{Z}_i), i denoting each of the N particles, and chooses the supremum from this set of eight values. The method assumes that the variables X_i, Y_i, and Z_i are not highly correlated.

In the nuclear case, the Kolmogorov 3D (that is, the Franceschini’s version) quantifies the discrepancy in the nucleons positions compared to the homogeneous case.

It is worth mentioning that the reliability of the 3D Kolmogorov statistic has been questioned in recent years [95]. The arguments, however, focus on the correct confidence intervals when applying the 3D Kolmogorov statistic to the null-hypothesis tests. Our investigation does not require computing these intervals, and thus, the questionings are irrelevant to the matter.

6.7.4 Minkowski functionals

Complex geometry may be decomposed into simple structures connected each other side by side. The basic 3D structure is the voxel. Thus, the complexity of any body may be associated to the way the voxels are connected between each other.

The connectivity pattern across voxels (required to represent a complex body) may be characterized by counting the number of unique vertices (n_v), edges (n_e), faces (n_f) and volume (n_c). However, a more practical set of indices corresponds to the set of the Minkowski functionals, defined as follows [96]:

\begin{array}{}\text{(G.11)}& \begin{array}{l}V=n_c\mathrm{,}S=-6n_c+2n_f\mathrm{,2}B=3n_c-2n_f+n_e\mathrm{,}\\ \chi =-n_c+n_f-n_e+n_v\mathrm{.}\end{array}\end{array}

V and S correspond to the body volume and surface, respectively. The B functional represents the mean breadth (or integral mean curvature) of the body. The positive values of B have been associated to spherical or cylindrical bodies (with no holes inside), while negative ones have been associated to bubble-like bodies or hollow cylinders [44, 97]. The B = 0 value (with χ = 0) appears to be associated to slab-like bodies.

An alternative expression for the Euler characteristic χ reads as follows:

\begin{array}{}\text{(G.12)}& \chi =\text{isolated}\text{regions}+\text{cavities}-\text{tunnels}\mathrm{.}\end{array}

This means that solid (including isolated voids inside) may be associated to positive values of χ. But sponge-like or jungle-gym structures may be associted to negative values χ values, since the number of tunnels exceeds the number of isolated regions.

A more detailed set of pastas associated to positive or negative signs of B,χ can be seen in Table 6 (more information in Ref. [21]). As mentioned above, the lasagnas (that is, slab-like structures) correspond to B~χ~0. The gnocchi (spherical-like) or spaghetti (cylinder-like) attain positive values of B and χ≥0. The jungle-gym (sponge-like) structures attain χ<0. The corresponding opposite structures in Table 6 are named with the prefix anti-.

6.7.5 The Minkowski voxels

The Minkowski functionals require the binning of space into “voxels” (that is, tridimensional “pixels”). Each voxel is supposed to include (approximately) a single nucleon. But this is somehow difficult to achieve if the system is not completely homogeneous (and regular).

We start with a simple cubic arrangement of 50% protons and 50% neutrons as shown in Fig. 91. The system is at the saturation density ρ₀=0.16 fm⁻³ (see caption for details).

Notice from Fig. 91 that the nearest distance between protons and neutrons is 1.85 fm. Likewise, the nearest distance between nucleons of the same species is approximately 2.65 fm (that is, the position of the second peak). It can be verified that the latter is approximately \sqrt{2}\times 1.85\text{fm}, as expected for the simple cubic arrangement (within the g(r) binning errors).

Figure 92(a) reproduces the same pattern for the distribution g(r) over all the nucleons. It further shows the Euler functional χ as a function of the voxel’s width (see caption for details). Both curves share the same abscissa for comparison reasons. For small values of d the functional χ is negative (not shown), meaning that the voxels are so small that tunnels prevail in the (discretized) system. At d = 1.65 fm this functional arrives to a maximum where cavities or isolated regions prevail. Some cavities may be “true” empty regions but other may simply be fake voids. The most probable ones, however, correspond to fake voids since the g(r) pattern actually presents a maximum at r>1.65 fm. Thus, increasing the voxel’s size will most probably cancel the fake cavities.

The particles located at the first maximum of g(r) (at the saturation density) may be envisaged as touching each other in a regular (simple cubic) array. Therefore, the mean radius for a nucleon should be 1.85/2 fm. This means, as a first thought, that binning the space into voxels of width 1.85 fm will include a single nucleon per voxels. This is, however, not completely true since approximately half of the first neighbors exceeds the 1.85 fm (see first peak in Fig. 91). Many voxels will be empty, and thus, a relevant probability of finding fake voids exists. Figure 92(a) illustrates this situation.

Notice that whenever a fake empty voxel exists, the contiguous one will perhaps host two nucleons. This is because the nucleon that exceeds the 1.85 fm distance to the neighbor, say on the left, may have shorten the distance to the neighbor on the right. Thus, the number of occupied voxels will probably not match the number of nucleons. Figure 92(b) shows a decrease in the number of occupied voxels (i.e., the volume) for r>1.6 fm.

The space binning should be done wider in order to avoid fake empty voxels. But, too wide voxels may include second neighbors. The most reasonable binning distance appears to be around the first minimum of g(r) (see Fig. 91). That is, at some point between 2.15 fm and 2.35 fm.

A reasonable criterion for the space binning may raise from the Euler functional: the right binning distance should drive the χ functional to unity, that is, to a single compact region. This occurs at d = 2.35 fm for sure, as can be seen in Fig. 93(b). It can be further checked from Fig. 92(b) that this binning allows the hosting of approximately two nucleons per voxels, meaning that the number of fake voids is negligible.

The binning width d = 2.35 fm may be further compared with other literature values. The nucleon radius used in Ref. [52] is 1.5 fm, although these authors realize that this value seems rather large. Indeed, we were able to associate a nucleon radius of 0.925 fm at the saturation density for a simple cubic arrangement. As mentioned above [and exhibited in Fig. 92(b)] the distance 1.5 fm corresponds (roughly) to the maximum distance that can host one nucleon per voxel (at the saturation density). Therefore, this binning width is right if maximum “contrast” is required, as in Ref. [52]. Our criterion, however, is somewhat more conservative since it requires \chi \simeq 1 (at the saturation density), as explained above.

6.7.6 Cluster recognition

The nucleon positions and momenta are used to identify the fragment structure of the system by means of the “Minimum Spanning Tree” (MST) cluster-detection algorithm of [98]. In summary, MST looks for correlations in configuration space: a particle i belongs to a cluster C if there is another particle j that belongs to C and |r_i-r_j|\le r_{cl}, where r_{cl} is a clusterization radius which, for the present study, was set to r_{cl}=3.0\text{fm}.

The main drawback of MST is that, since only correlations in r-space are used, it neglects completely the effect of momentum giving incorrect information for dense systems and for highly dynamical systems such as those formed in colliding nuclei. The accuracy of MST is somehow restricted to cold systems and small momentum transfers; whenever these conditions are not met, or the evolution of the process is expected to be highly dynamical, the MST will not provide useful information, and such systems must be analyzed by means of more computationally expensive algorithms, such as the Minimum Spanning Tree Energy (MSTE) algorithm or the Early Cluster Recognition Algorithm (ECRA) [99, 100]. However, in the case of relatively cold systems, such as nuclear crusts, the MST is sufficient. In our case of periodic boundary conditions, the MST detection of fragments has been modified to take into account the image cells and recognize fragments that extend into adjacent cells.

Figure 94 shows an example of the size distribution of the clusters obtained for a case with 3333 nucleons, x = 0.3, and at T = 0.3 MeV and ρ = 0.009 fm⁻³. The inset shows a projection of the position of the nucleons within the cell. The shown structure was obtained with the “screened Coulomb” treatment that is described in Section 3.1.

CHANGE HISTORY

A Correction to this paper has been published:
https://doi.org

The Retraction Note was updated to reflect reception of statements of disagreement from all authors.

This article has been retracted.Please see the Retraction Notice for more detail:
https://doi.org

References

[1] D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Structure of matter below nuclear saturation density, Phys. Rev. Lett. 50(26), 2066 (1983)

[2] K. Oyamatsu, Nuclear shapes in the inner crust of a neutron star, Nucl. Phys. A 561(3), 431 (1993)

[3] T. Maruyama, K. Niita, K. Oyamatsu, T. Maruyama, S. Chiba, and A. Iwamoto, Quantum molecular dynamics approach to the nuclear matter below the saturation density, Phys. Rev. C 57(2), 655 (1998)

[4] C. P. Lorenz, D. G. Ravenhall, and C. J. Pethick, Neutron star crusts, Phys. Rev. Lett. 70(4), 379 (1993)

[5] K. S. Cheng, C. C. Yao, and Z. G. Dai, Properties of nuclei in the inner crusts of neutron stars in the relativistic mean-field theory, Phys. Rev. C 55(4), 2092 (1997)

[6] T. Kido, T. Maruyama, K. Niita, and S. Chiba, MD simulation study for nuclear matter, Nucl. Phys. A 663-664, 877c (2000)

[7] G. Watanabe, K. Iida, and K. Sato, Thermodynamic properties of nuclear “pasta” in neutron star crusts, Nucl. Phys. A 676(1–4), 455 (2000)

[8] R. D. Williams and S. E. Koonin, Sub-saturation phases of nuclear matter, Nucl. Phys. A 435(3–4), 844 (1985)

[9] M. Hashimoto, H. Seki, and M. Yamada, Shape of Nuclei in the Crust of Neutron Star, Prog. Theor. Phys. 71(2), 320 (1984)

[10] P. N. Alcain and C. O. Dorso, The neutrino opacity of neutron rich matter, Nucl. Phys. A 961, 183 (2017)

[11] D. Page, J. M. Lattimer, M. Prakash and A. W. Steiner, Minimal Cooling of Neutron Stars: A New Paradigm, Astrophys. J. Suppl. 155, 623 (2004)

[12] B. Schuetrumpf, G. Martínez-Pinedo, M. Afibuzzaman, and H. M. Aktulga, Survey of nuclear pasta in the intermediate-density regime: Shapes and energies, Phys. Rev. C 100(4), 045806 (2019)

[13] B. Schuetrumpf, G. Martínez-Pinedo, and P. G. Reinhard, Survey of nuclear pasta in the intermediate-density regime: Structure functions for neutrino scattering, Phys. Rev. C 101(5), 055804 (2020)

[14] G. Watanabe, K. Sato, K. Yasuoka, and T. Ebisuzaki, Microscopic study of slablike and rodlike nuclei: Quantum molecular dynamics approach, Phys. Rev. C 66(1), 012801 (2002)

[15] C. J. Horowitz, M. A. Perez-García, and J. Piekarewicz, Neutrino-“pasta” scattering: The opacity of nonuniform neutron-rich matter, Phys. Rev. C 69(4), 045804 (2004)

[16] B. Schuetrumpf and W. Nazarewicz, Twist-averaged boundary conditions for nuclear pasta Hartree–Fock calculations, Phys. Rev. C 92(4), 045806 (2015)

[17] F. J. Fattoyev, C. J. Horowitz, and B. Schuetrumpf, Quantum nuclear pasta and nuclear symmetry energy, Phys. Rev. C 95(5), 055804 (2017)

[18] C. O. Dorso, P. A. Giménez Molinelli, and J. A. López, in: “Neutron Star Crust”, Eds. C. A. Bertulani and J. Piekarewicz, Nova Science Publishers, ISBN 978-1620819029 (2012)

[19] P. N. Alcain, P. A. Giménez Molinelli, and C. O. Dorso, Beyond nuclear “pasta”: Phase transitions and neutrino opacity of new “pasta” phases, Phys. Rev. C 90(6), 065803 (2014)

[20] C. J. Horowitz, M. A. Pérez-García, J. Carriere, D. K. Berry, and J. Piekarewicz, Nonuniform neutron-rich matter and coherent neutrino scattering, Phys. Rev. C 70(6), 065806 (2004)

[21] C. O. Dorso, P. A. Giménez Molinelli, and J. A. López, Topological characterization of neutron star crusts, Phys. Rev. C 86(5), 055805 (2012)

[22] I. Tanihata, Preprint RIKEN-AF-NP-229, 1996; P. W. Zhao, Z. P. Li, J. M. Yao, and J. Meng, New parametrization for the nuclear covariant energy density functional with a point-coupling interaction, Phys. Rev. C 82(5), 054319 (2010); https://doi.org/10.1103/PhysRevC.85.035201 M. Dutra, O. Lourenço, J. S. Sá Martins, A. Delfino, J. R. Stone, and P. D. Stevenson, Skyrme interaction and nuclear matter constraints, Phys. Rev. C 85(3), 035201 (2012)

[23] S. Kumar and Y. G. Ma, Investigation of compressibilities using neutron-rich projectileʼs fragmentation at intermediate energy, Nucl. Phys. A 898, 59 (2013)

[24] P. Danielewicz, R. Lacey, and W. G. Lynch, Determination of the equation of state of dense matter, Science 298(5598), 1592 (2002)

[25] W. D. Myers and W. J. Swiatecki, The nuclear Thomas–Fermi model, Acta Phys. Pol. B 26, 111 (1995)

[26] A. Barrañón, J. Escamilla Roa, and J. A. López, Entropy in the nuclear caloric curve, Phys. Rev. C 69(1), 014601 (2004)

[27] P. J. Siemens, Liquid–gas phase transition in nuclear matter, Nature 305(5933), 410 (1983); https://doi.org/10.1038/336110a0 P. J. Siemens, Macroscopic behaviour of nuclear matter, Nature 336(6195), 110 (1988)

[28] J. A. López and C. O. Dorso, Lecture Notes on Phase Transitions in Nuclear Matter, World Scientific, 2000

[29] H. Müller and B. Serot, Phase transitions in warm, asymmetric nuclear matter, Phys. Rev. C 52(4), 2072 (1995)

[30] J. A. López, A. Gaytán Terrazas, and S. Terrazas Porras, Isospin-dependent phase diagram of nuclear matter, Nucl. Phys. A 994, 121664 (2020)

[31] See, e.g., https://www.ks.uiuc.edu/Research/vmd/current/ug/node73.html, retrieved Sep. 2, 2019

[32] J. A. López, E. Ramírez-Homs, R. González, and R. Ravelo, Isospin-asymmetric nuclear matter, Phys. Rev. C 89(2), 024611 (2014)

[33] J. A. López and S. Terrazas Porras, Symmetry energy in the liquid–gas mixture, Nucl. Phys. A 957, 312 (2017)

[34] K. Hagel, J. B. Natowitz, and G. Röpke, The equation of state and symmetry energy of low-density nuclear matter, Eur. Phys. J. A 50(2), 39 (2014)

[35] S. Kowalski, J. B. Natowitz, S. Shlomo, R. Wada, K. Hagel, J. Wang, T. Materna, Z. Chen, Y. G. Ma, L. Qin, A. S. Botvina, D. Fabris, M. Lunardon, S. Moretto, G. Nebbia, S. Pesente, V. Rizzi, G. Viesti, M. Cinausero, G. Prete, T. Keutgen, Y. E. Masri, Z. Majka, and A. Ono, Experimental determination of the symmetry energy of a low density nuclear gas, Phys. Rev. C 75(1), 014601 (2007)

[36] R. Wada, K. Hagel, L. Qin, J. B. Natowitz, Y. G. Ma, G. Röpke, S. Shlomo, A. Bonasera, S. Typel, Z. Chen, M. Huang, J. Wang, H. Zheng, S. Kowalski, C. Bottosso, M. Barbui, M. R. D. Rodrigues, K. Schmidt, D. Fabris, M. Lunardon, S. Moretto, G. Nebbia, S. Pesente, V. Rizzi, G. Viesti, M. Cinausero, G. Prete, T. Keutgen, Y. El Masri, and Z. Majka, Nuclear matter symmetry energy at 0.03≤ρ/ρ₀≤0.2, Phys. Rev. C 85(6), 064618 (2012)

[37] L. W. Chen, C. M. Ko, and B. A. Li, Isospin-dependent properties of asymmetric nuclear matter in relativistic mean field models, Phys. Rev. C 76(5), 054316 (2007)

[38] E. L. Medeiros and J. Randrup, Thermostatic properties of Seyler–Blanchard nuclei, Phys. Rev. C 45(1), 372 (1992)

[39] C. J. Horowitz and A. Schwenk, Cluster formation and the virial equation of state of low-density nuclear matter, Nucl. Phys. A 776(1–2), 55 (2006)

[40] J. Xu, L. W. Chen, B. A. Li, and H. R. Ma, Temperature effects on the nuclear symmetry energy and symmetry free energy with an isospin and momentum dependent interaction, Phys. Rev. C 75(1), 014607 (2007)

[41] P. A. Giménez Molinelli, J. I. Nichols, J. A. López, and C. O. Dorso, Simulations of cold nuclear matter at sub-saturation densities, Nucl. Phys. A 923, 31 (2014)

[42] A. Vicentini, G. Jacucci, and V. R. Pandharipande, Fragmentation of hot classical drops, Phys. Rev. C 31(5), 1783 (1985); https://doi.org/10.1103/PhysRevC.34.177 R. J. Lenk and V. R. Pandharipande, Disassembly of hot classical charged drops, Phys. Rev. C 34(1), 177 (1986); https://doi.org/10.1103/PhysRevC.42.372 R. J. Lenk, T. J. Schlagel, and V. R. Pandharipande, Accuracy of the Vlasov–Nordheim approximation in the classical limit, Phys. Rev. C 42(1), 372 (1990)

[43] G. Raciti, R. Bassini, M. Begemann-Blaich, S. Fritz, S. J. Gaff, N. Giudice, C. Gross, G. Immé, I. Iori, U. Kleinevoss, G. J. Kunde, W. D. Kunze, U. Lynen, M. Mahi, T. Möhlenkamp, W. F. J. Müller, B. Ocker, T. Odeh, J. Pochodzalla, G. Riccobene, F. P. Romano, A. Sajia, M. Schnittker, A. Schüttauf, C. Schwarz, W. Seidel, V. Serfling, C. Sfienti, W. Trautmann, A. Trzcinski, G. Verde, A. Wörner, H. Xi, and B. Zwieglinski, A systematic study of the nuclear caloric curve, Nuovo Cim. 111(8-9), 987 (1998)

[44] H. Sonoda, G. Watanabe, K. Sato, K. Yasuoka, and T. Ebisuzaki, Phase diagram of nuclear “pasta” and its uncertainties in supernova cores, Phys. Rev. C 77(3), 035806 (2008)

[45] C. O. Dorso, G. Frank, and J. A. López, Phase transitions and symmetry energy in nuclear pasta, Nucl. Phys. A 978, 35 (2018)

[46] C. J. Horowitz, Links between heavy ion and astrophysics, Eur. Phys. J. A 30(1), 303 (2006)

[47] G. Watanabe and K. Iida, Electron screening in the liquid–gas mixed phases of nuclear matter, Phys. Rev. C 68(4), 045801 (2003)

[48] T. Maruyama, T. Tatsumi, D. N. Voskresensky, T. Tanigawa, and S. Chiba, Nuclear “pasta” structures and the charge screening effect, Phys. Rev. C 72(1), 015802 (2005)

[49] C. J. Horowitz, M. A. Perez-Garcia, D. K. Berry, and J. Piekarewicz, Dynamical response of the nuclear “pasta” in neutron star crusts, Phys. Rev. C 72(3), 035801 (2005)

[50] J. Piekarewicz and G. T. Sánchez, Proton fraction in the inner neutron-star crust, Phys. Rev. C 85(1), 015807 (2012)

[51] J.A. López and E. Ramírez-Homs, Effect of an electron gas on a neutron-rich nuclear pasta, Nuc. Sci. Tech. 26, S20502 (2015)

[52] A. S. Schneider, C. J. Horowitz, J. Hughto, and D. K. Berry, Nuclear “pasta” formation, Phys. Rev. C 88(6), 065807 (2013)

[53] K. Binder, B. J. Block, P. Virnau, and A. Tröster, Beyond the van der Waals loop: What can be learned from simulating Lennard–Jones fluids inside the region of phase coexistence, Am. J. Phys. 80(12), 1099 (2012)

[54] C. J. Horowitz, D. K. Berry, C. M. Briggs, M. E. Caplan, A. Cumming, and A. S. Schneider, Disordered nuclear pasta, magnetic field decay, and crust cooling in neutron stars, Phys. Rev. Lett. 114(3), 031102 (2015)

[55] C. Dorso, G. Frank, and J. A. López, Symmetry energy in neutron star matter, Nucl. Phys. A 984, 77 (2019)

[56] J. A. López, J. A. Muñoz, C. O. Dorso, and G. Frank, Machine learning Minkoswki functionals of neutron star crusts, J. Phys. Conf. Ser. (2019); J. A. López and J. A. Muñoz, Analytical expression and neural network study of the symmetry energy, CERN Proc. 1, 29 (2019)

[57] P. N. Alcain, Dependencia en el isospín de la ecuación de estado de la materia nuclear, Ph.D. Thesis, Universidad de Buenos Aires, 2019

[58] D. Frenkel y B. Smit, Understanding Molecular Simulations, 2nd Ed., Academic Press, 2002

[59] A. Deibel, A. Cumming, E. F. Brown, and S. Reddy, Late-time cooling of neutron star transients and the physics of the inner crust, Astrophys. J. 839(2), 95 (2017)

[60] E. F. Brown, A. Cumming, F. J. Fattoyev, C. J. Horowitz, D. Page, and S. Reddy, Rapid neutrino cooling in the neutron star MXB 1659-29, Phys. Rev. Lett. 120(18), 182701 (2018)

[61] A. S. Schneider, D. K. Berry, M. E. Caplan, C. J. Horowitz, and Z. Lin, Effect of topological defects on “nuclear pasta” observables, Phys. Rev. C 93(6), 065806 (2016)

[62] R. Nandi and S. Schramm, Transport properties of the nuclear pasta phase with quantum molecular dynamics, Astrophys. J. 852(2), 135 (2018)

[63] C. J. Horowitz, and D. K. Berry, Shear viscosity and thermal conductivity of nuclear “pasta”, Phys. Rev. C 78(3), 035806 (2008)

[64] J. M. Dunn, Nanoscale phonon thermal conductivity via molecular dynamics, Ph.D. Thesis, Purdue University, 2016

[65] F. Müller–Plathe, A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity, J. Chem. Phys. 106(14), 6082 (1997)

[66] A. Barrañón, C. O. Dorso, J. A. López, and J. Morales, LATINO: A semi-classical model to study nuclear fragmentation, Rev. Mex. Fis. 45(suppl. 2), 110 (1999)

[67] A. Chernomoretz, L. Gingras, Y. Larochelle, L. Beaulieu, R. Roy, C. St-Pierre, and C. O. Dorso, Quasiclassical model of intermediate velocity particle production in asymmetric heavy ion reactions, Phys. Rev. C 65(5), 054613 (2002)

[68] A. Barrañón, C. O. Dorso, and J. A. López, Searching for criticality in nuclear fragmentation, Rev. Mex. Fís. 47(sup. 2), 93 (2001)

[69] A. Barrañón, C. O. Dorso, and J. A. López, Time dependence of isotopic temperatures, Nucl. Phys. A 791(1–2), 222 (2007)

[70] A. Barrañón, R. Cárdenas, C. O. Dorso, and J.A. López, The critical exponent of nuclear fragmentation, Acta Physica Hungarica A: Heavy Ion Phys. 17(1), 59 (2003)

[71] C. O. Dorso and J. A. López, Selection of critical events in nuclear fragmentation, Phys. Rev. C 64(2), 027602 (2001)

[72] A. Barrañón, J. Escamilla Roa, and J. A. López, The transition temperature of the nuclear caloric curve, Braz. J. Phys. 34(3A), 904 (2004)

[73] C. O. Dorso, C. R. Escudero, M. Ison, and J. A. López, Dynamical aspects of isoscaling, Phys. Rev. C 73(4), 044601 (2006)

[74] C. A. Dorso, P. A. G. Molinelli, and J. A. López, Isoscaling and the nuclear EoS, J. Phys. G 38(11), 115101 (2011); C. O. Dorso, P. A. G. Molinelli, and J. A. López, Searching for the origin of isoscaling: Confinement and expansion, Rev. Mex. Phys. S57 (1), 14 (2011)

[75] T. M. Nymand and P. Linse, Ewald summation and reaction field methods for potentials with atomic charges, dipoles, and polarizabilities, J. Chem. Phys. 112, 6152 (2000)

[76] P. N. Alcain, P. A. Giménez Molinelli, J. I. Nichols, and C. O. Dorso, Effect of Coulomb screening length on nuclear “pasta” simulations, Phys. Rev. C 89(5), 055801 (2014)

[77] B. L. Holian, A. F. Voter, and R. Ravelo, Thermostatted molecular dynamics: How to avoid the Toda demon hidden in Nosé–Hoover dynamics, Phys. Rev. E 52(3), 2338 (1995)

[78] H. C. Andersen, Molecular dynamics simulations at constant pressure and/or temperature, J. Chem. Phys. 72(4), 2384 (1980)

[79] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81(1), 511 (1984)

[80] J. A. López, S. Terrazas Porras, and A. Rodríguez Gutiérrez, Thermodynamics of neutron-rich nuclear matter, AIP Conf. Proc. 1753, 050001 (2016)

[81] B. A. Li, L. W. Chen, and C. M. Ko, Recent progress and new challenges in isospin physics with heavy-ion reactions, Phys. Rep. 464(4–6), 113 (2008)

[82] B. A. Li, A. Ramos, G. Verde, and I. Vidana, Topical issue on Nuclear Symmetry Energy, Eur. Phys. J. A 50(2), 9 (2014)

[83] J. B. Natowitz, G. Röpke, S. Typel, D. Blaschke, A. Bonasera, K. Hagel, T. Klähn, S. Kowalski, L. Qin, S. Shlomo, R. Wada, and H. H. Wolter, Symmetry energy of dilute warm nuclear matter, Phys. Rev. Lett. 104(20), 202501 (2010)

[84] S. Typel, H. H. Wolter, G. Röpke, and D. Blaschke, Effects of the liquid–gas phase transition and cluster formation on the symmetry energy, Eur. Phys. J. A 50(2), 17 (2014)

[85] M. Dutra, O. Lourenço, J. S. Sá Martins, A. Delfino, J. R. Stone, and P. D. Stevenson, Skyrme interaction and nuclear matter constraints, Phys. Rev. C 85(3), 035201 (2012)

[86] M. Dutra, O. Lourenço, S. S. Avancini, B. V. Carlson, A. Delfino, D. P. Menezes, C. Providência, S. Typel, and J. R. Stone, Relativistic mean-field hadronic models under nuclear matter constraints, Phys. Rev. C 90(5), 055203 (2014)

[87] M. Colonna, V. Baran, M. D. Toro, and H. H. Wolter, Isospin distillation with radial flow: A test of the nuclear symmetry energy, Phys. Rev. C 78(6), 064618 (2008)

[88] Y. Zhou, B. Anglin, and A. Strachan, Phonon thermal conductivity in nanolaminated composite metals via molecular dynamics, J. Chem. Phys. 127(18), 184702 (2007)

[89] J. Dunn, E. Antillon, J. Maassen, M. Lundstrom, and A. Strachan, Role of energy distribution in contacts on thermal transport in Si: A molecular dynamics study, J. Appl. Phys. 120(22), 225112 (2016)

[90] K. H. Lin and A. Strachan, Thermal transport in SiGe superlattice thin films and nanowires: Effects of specimen and periodic lengths, Phys. Rev. B 87(11), 115302 (2013)

[91] F. A. Lindemann, The calculation of molecular vibration frequencies, Phys. Z. 11, 609 (1910)

[92] Z. W. Birnbaum, Numerical tabulation of the distribution of Kolmogorov’s statistic for finite sample size, J. Am. Stat. Assoc. 47(259), 425 (1952)

[93] E. Gosset, A three-dimensional extended Kolmogorov–Smirnov test as a useful tool in astronomy, Astron. Astrophys. 188, 258 (1987)

[94] G. Fasano and A. Franceschini, A multidimensional version of the Kolmogorov–Smirnov test, Mon. Not. R. Astron. Soc. 225(1), 155 (1987)

[95] G. J. Babu and E. D. Feigelson, Astronomical Data Anal ysis Software and Systems XV, Eds. C. Gabriel, et al., ASP Conference Series, 351, 127 (2006)

[96] K. Michielsen and H. De Raedt, Integral-geometry morphological image analysis, Phys. Rep. 347(6), 461 (2001)

[97] B. Schuetrumpf, M. A. Klatt, K. Iida, J. A. Maruhn, K. Mecke, and P. G. Reinhard, Time-dependent Hartree–Fock approach to nuclear “pasta” at finite temperature, Phys. Rev. C 87(5), 055805 (2013)

[98] A. Strachan and C. O. Dorso, Time scales in fragmentation, Phys. Rev. C 55(2), 775 (1997); https://doi.org/10.1103/PhysRevC.56.995 A. Strachan and C. O. Dorso, Fragment recognition in molecular dynamics, Phys. Rev. C 56(2), 995 (1997)

[99] C. O. Dorso and J. Randrup, Early recognition of clusters in molecular dynamics, Phys. Lett. B 301(4), 328 (1993)

[100] P. N. Alcain and C. O. Dorso, Dynamics of fragment formation in neutron-rich matter, Phys. Rev. C 97(1), 015803 (2018)

[101] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117(1), 1 (1995)

Rights & permissions

Higher Education Press 2021