Effect of an electron gas on a neutron-rich nuclear pasta

Special Section on International Workshop on Nuclear Dynamics in Heavy-Ion Reactions (IWND 2014)

Effect of an electron gas on a neutron-rich nuclear pasta

Alberto López Jorge，

Ramírez-Homs Enrique

Nuclear Science and Techniques

Vol.26, No.2

Article number S20502

Published in print 20 Apr 2015

Available online 20 Apr 2014

DOI：10.13538/j.1001-8042/nst.26.S20502

59300

We studied the role that electron gas has on the formation of nuclear structures at subsaturation densities and low temperatures (T < 1 MeV). Using a classical molecular dynamics model we studied isospin symmetric and asymmetric matter at subsaturation densities and low temperatures varying the Coulomb interaction strength. The effect of such variation was quantified on the fragment size multiplicity, the inter-particle distance, the isospin content of the clusters, the nucleon mobility and cluster persistence, and on the nuclear structure shapes. We found that the presence of an electron gas distributes matter more evenly, disrupts the formation of larger objects, reduces the isospin content, and modifies the nucleon average displacement, but does not affect the inter-nucleon distance in clusters. The nuclear structures are also found to change shapes by different degrees depending on their isospin content, temperature and density.

Electron gasMolecular dynamicsNeutron rich matterNuclear pasta

I. INTRODUCTION

Understanding nuclear matter at subcritical temperatures and subsaturation densities is of importance to understand phenomena ranging from heavy-ion collisions to the stability of stellar objects. This region has been investigated through nuclear reactions leading to the fragmentation of the nuclei, and it is now known that excited nuclear matter can exhibit a first-order phase transition [1], spinodal decomposition [2], critical behavior [3] and other phenomenon associated with liquid-gas systems. More recently, such studies have been extended to neutron-rich nuclei [4-7].

At lower temperatures (T 1 MeV) and subsaturation densities, theoretical studies predict that nuclear matter forms non-homogeneous structures known as nuclear pasta. Synoptically, static models [8-13] (e.g. liquid-drop model, Thomas-Fermi, and Hartree-Fock) find such structures from a balance between the different components of the nuclear and Coulomb energies, while dynamical models [14-23] (e.g. quantum and classical molecular dynamics) find them as equilibrium states of dynamical evolutions of such systems.

The non-homogeneous phases are expected to exist on the outer crust of neutron stars and play a role in the emission of neutrinos [20], starquakes and pulsar glitches [24]. In such stellar conditions an unbalance between the weak decay $n \to p + e^{-} + {\bar{ν}}_{e}$ and its opposite process, p + e n+ν is expected to yield an accumulation of neutrons and produce a neutralizing electron gas. Thus researchers have studied such nuclear structures while embedded in an electron gas. Although most calculations have shown that the structures are formed thanks to the interplay of nuclear forces with the Coulomb forces of the electron gas, none have dissected this effect quantitatively. The question we aim to answer is, can the nuclear pasta structures form without the embedding electron gas?

This investigation will focus on the role of the Coulomb forces on the creation of the pasta using a classical molecular dynamics which was shown to be effective in these type of studies [21]. By using a microscopic view of the nucleons and a combination of topological tools the different structures formed can be characterized into recognizable patterns, opening the door to a quantifiable analysis of the effect the Coulomb forces of the electron gas have on the pasta formation.

In 2003 [23] a static liquid-drop model was used to study the effect of an electron gas on cold nuclear structures, and it was found that its main effect was to extend the range of densities where bubbles and clusters appear. Later, in a 2005 study [25] 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. The authors remarked the importance of extending such study to finite temperatures and with dynamical models.

Recently, Monte Carlo simulations, which approximated the Coulomb interaction via an Ewald summation [20, 26-28], determined the lowest proton fractions that are compatible with β equilibrium in neutron star crust environments (at T= 1.0 MeV). The authors recommended that further investigation be performed with molecular dynamics simulations to study dynamical observables and to examine the role of the temperature in the study of nuclear matter under subsaturation density and low temperature conditions.

The goal of this work is to study the effect that electron gas has on the formation of the nuclear pasta using a dynamical model, at non-zero temperatures, and in a systematic way that takes advantage of the tools developed in a previous work [21]. It must be remarked that it is not the intention of this work to study either the structure or other properties of neutron star crusts nor to propose the present molecular dynamics method for such purpose; the only goal of the present exercise is –again– to understand the role the electron gas plays in determining the pasta shapes. This will be achieved through the calculation of the pasta structures with and without the effect of the electron gas as well as with softened interactions which –of course– do not exist in nature. Reiterating, we are not hoping to mimic neutron star crusts nor proposing the existence of nuclear structures with weakened electron gas interactions.

In the next section the model used (classical molecular dynamics) is briefly introduced along with the analysis tools employed. Section III presents a series of results of the simulations and analyses performed. We summarize some of our findings in Section IV and present a quick look at work in progress.

II. MODEL USED

The model used here was developed to study nuclear reactions from a microscopic point of view, and it is composed of a molecular dynamics code retrofitted with cluster recognition algorithms and a plethora of analysis tools. The model was first used for infinite nuclear matter elsewhere [7, 21, 22], here we simply mention some basic ingredients of the model.

The classical molecular dynamics model (CMD), first introduced in [29], has been used in heavy-ion reaction studies to explain experimental data [30], identify phase transitions and other critical phenomena [31-34], explore the caloric curve of nuclear matter [35, 36] and isoscaling [37, 38]. CMD uses two two-body potentials to describe the motion of "nucleons" by solving their classical equations of motion. The potentials were developed phenomenologically by Pandharipande [39], we directinterested readers to this original work for a review of the basic properties of nuclear matter and finite nuclei that can be described by this model. The potential used is:

\begin{array}{l} V_{n p} (r) = & V_{r} [\exp (- μ_{r} r) / r - \exp (- μ_{r} r_{c}) / r_{c}] \\ - V_{a} [\exp (- μ_{a} r) / r - \exp (- μ_{a} r_{c}) / r_{c}], \\ V_{N N} (r) & = V_{0} [\exp (- μ_{0} r) / r - \exp (- μ_{0} r_{c}) / r_{c}], \end{array}

(1)

where Vnp is the potential between a neutron and a proton, and VNN is the repulsive interaction between either nn or pp. The cutoff radius is rc=5.4 fm and for r>r_c both potentials are set to zero. The Yukawa parameters μr, μa and μ₀ were determined by Pandharipande to yield an equilibrium density of ρ₀=0.16 fm^-3, a binding energy E(ρ₀)=16 MeV/nucleon and a selectable compressibility of 250 MeV for the "Medium" model, and 535 MeV for the "Stiff" ’ [39]. We remark that, since the potential VNN does not permit bound states between identical nucleons, pure neutron matter is unbound and, differing from potentials used by other models [20], the Pandharipande potentials have a hard core.

The main advantage of the CMD model is the possibility of knowing the position and momentum of all particles at all times, as this allows the study of the structure of the nuclear medium from a microscopic point of view. The outcome of CMD, namely, the time evolution of the particles in (r,p), can be used to study the resulting pasta structures.

On the other hand, the CMD model fails to comply with several quantum effects. Luckily, the model can still be applied for the nucleon dynamics since, as explained elsewhere [21, 40, 41], the mean inter-particle distance (V/N)^1/3 in nuclear matter in the region of interest (ρ 0.1 fm^-3, T 1 MeV) is larger than the nucleon thermal wavelength [7]. Other flaws of the CMD are to be discussed after the results are presented.

To simulate an infinite medium, systems with thousands of nucleons were constructed using CMD under periodic boundary conditions. Cases symmetric in isospin (i.e. with x=Z/A=0.5, 1000 protons and 1000 neutrons), and cases with asymmetric isospin (with x=0.3 comprised of 1000 protons and 2000 neutrons) were constructed in cubical boxes with sizes adjusted to have densities between 0.01 fm^-3 ρ ρ₀; it must be remarked that by fixing the mass and size of each cell, the scale of the density fluctuations is also determined. Our previous study, [21], presented an smörgåsboard of structures obtained through this method.

A. Coulomb interaction

To study nuclear matter embedded in a degenerate electron gas, the system is first taken as overall neutral (β-equilibrated), and then a proper approximation is used to take into account the Coulomb interaction (which is of infinite range). The two most common approaches are the Thomas-Fermi screened Coulomb potential (used, e.g., in QMD [15]) and the Ewald summation procedure [42]; CMD has been used with both approximations (see [21] for a comparison of methods under CMD) but for computational reasons we opt for the use of a screened Coulomb potential.

We approximate the electron gas as a uniform ideal Fermi gas with the same number density as the protons. Its inclusion in the CMD is through the use of a screened Coulomb potential obtained from the Poisson equation, $V_{C}^{(S c r)} (r) = (e^{2} / r) \exp (- r / λ)$ , where the screening length, λ, is related to the electron mass, Fermi momentum, and number density. Here, pragmatically, and after testing various lengths, we used the prescription used before [15, 20] and set λ=10 fm, which satisfies the requirement of being sufficiently smaller than the size of the simulation cell, L=(A/ρ)^1/3.

B. Simulation procedure

The trajectories of the nucleons are then governed by the Pandharipande and the screened Coulomb potentials. The equations of motion are solved using a symplectic Verlet integration with energy conservation of 𝒪(0.01 %). The nuclear system is force-heated or cooled using isothermal molecular dynamics with the Andersen thermostat procedure [43] with small temperature steps. We are interested in the range of 0.1 MeV T 1.0 MeV and densities 0.01 fm^-3 ρ ρ₀. In comparison to heavy-ion reactions, these conditions correspond to practically semi-frozen states.

To obtain reliable statistics we sample each T, ρ and x configuration 200 times in a stationary ergodic process. That is, after reaching equilibrium at a given set of values of T, ρ and x, the (r,p) values of all nucleons are stored for further processing while the system continues its evolution for a "long" time (4000 time steps) until it reaches a state independent of the previous one, point at which the procedure is repeated for a total of 200 times. The recorded nucleon positions and momenta are then used later to identify clusters and to characterize the structure as explained next.

Given that, for the densities and temperatures of interest, the nucleons do not have large mobilities, the task of identifying self-bound cumuli is sufficiently simple as to be achieved with a basic "Minimum Spanning Tree" (MST) algorithm. The MST finds cluster membership simply by identifying those particles that are closer to each other than some clusterization radius, which here is set to 3.0 fm based on the range of the nuclear potential and densities investigated. The generic method was modified to recognize fragments that extend into any of the adjacent periodic cells.

C. Analysis tools

Once the fragments have been identified, their global properties are studied as follows. MST yields information about the fragment multiplicity. In Ref. [21] systems at x=0.3 and 0.5 and at various densities (0.01 fm^-3 ρ 0.015 fm^-3) and temperatures (0.1 MeV T 1.0 MeV) exhibited increases in multiplicity with temperature (in opposition to the rise-fall that takes place in intermediate-energy heavy-ion reactions), as well as the formation of percolating structures that encompass most of the nucleons in the simulation.

At a microscopic level, the dynamics of the nucleons in such systems can be quantified through their average displacement as a function of "time", i.e. through the time steps of the simulation; for the cases mentioned before (cf. [21]) the mobility between configurations was found to be between 2 and 20 fm. Likewise, the microscopic stability of the clusters can be gauged through the "persistency" [40] which measures the tendency of members of a given cluster to remain in the same cluster. For the cases previously studied, colder clusters (T 0.5 MeV) exchange less than half of their nucleons during 200 ergodic configurations, while those at higher temperatures (T 0.8 MeV) replace up to 90 Another interesting descriptor is the isospin content of each cluster produced, i.e. their x values. In the previous study it was found that in the case of x=0.3, small clusters (A≲10) tend to have less protons than the global average, while such effect was not present in the case of x=0.5. Another global characterization tool is the pair correlation function, g(r), which gives information about the spatial ordering of the nuclear medium by comparing the average local density to the global density. In the previous study g(r) showed that nucleons in clusters have an interparticle distance of about 1.7fm at all studied densities, and was useful in determining the onset of lasagna-like structures.

Proceeding beyond global measures, the shapes of nuclear structures can be characterized by their volume, surface area, mean curvature, plus an interesting construct known as the Euler characteristic. These four objects comprise what is known as the "Minkowsky functionals" which completely describe all topological properties of any three-dimensional object. The evaluation of the mean curvature (averaged over the whole surface of the object) and Euler number (a topological invariant connected to the vertices, edges and faces of polyhedra), however, requires the mapping of the nuclear clusters into a polyhedra. This procedure can be accomplished through the Michielsen–De Raed algorithm [18, 44].

In [21] it was shown that generic structures, such as "gnocchi", "spaghetti", "lasagna" and "crossed-lasagnas" or "jungle gym" and their inverse structures (with voids replacing particles and vice versa), all have well defined and distinct values of the mean curvature and Euler numbers. It should be remarked that the assignment of curvature and Euler values to a given structure becomes relative to the overall scale. This becomes important for structures with near-zero Euler number which signal spaghetti, lasagna and their anti-structures. Curled spaghetti rods of a given length, for instance, can easily oscillate between positive and negative near zero Euler numbers by curling or uncurling a segment. This classification, which was respected at all x values, densities and temperatures studied, is shown in Table 1.

Classification Curvature - Euler

	Curvature <0	Curvature ~ 0	Curvature >0
Euler >0	Anti-Gnocchi		Gnocchi
Euler ∼0	Anti-Spaghetti	Lasagna	Spaghetti
Euler <0	Anti-Jungle Gym		Jungle Gym

It is with this model and arsenal of tools that we now embark on a study of the effect that electron gas has on the formation of the nuclear pasta.

III. THE EFFECT OF THE ELECTRON GAS

Along the lines of previous studies [23, 25, 28] and adhering to their advise, in this study we used a molecular dynamics model to study spatial and dynamical effects the electron gas may produce on the nuclear structures at non-zero temperatures. Specifically we extended the approach of one of the previous studies and compared pasta structures obtained with and without the electron gas as well as with varying strengths.

In summary, pasta structures were obtained at subsaturation densities (0.015 fm^-3 ρ 0.072 fm^-3) 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 {ρ,T,x}. 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 $α V_{C}^{(S c r)}$ , where 0 α 1. For a fair comparison, all corresponding cases with different values of α were produced with identical initial conditions of {x,ρ,T} and started off from the same initial random configuration. As an illustration, Fig. 1 shows two corresponding structures obtained with and without the electron gas.

Fig. 1.

(Color online) Characteristic structures obtained with (α=1) and without (α=0) the proton-electron gas interaction for symmetric matter (x=0.5) at density ρ=0.015 fm^-3 and temperature T=0.1 MeV.

The effect of the screened potential can also be observed on the fragment size multiplicity. Fig. 2 shows the distribution of cluster sizes observed in 200 configurations of asymmetric matter (x=0.3) at density ρ=0.015 fm^-3 and temperatures T=0.1 MeV and T= 1.0 MeV. The figures correspond to configurations with and without the Coulomb interaction. The increase in the number of fragments of sizes A 500 as α goes from α=1 to α=0 underline the role the electron gas has on determining the mass distribution.

Fig 2.

Effect of the screened potential on the fragment size multiplicity. Plots show the distribution of cluster sizes observed in 200 configurations of x=0.3 nuclear matter at density ρ=0.015 fm^-3 and temperatures T= 1.0 MeV (top) and T= 0.1 MeV (bottom). The figure on the left panel correspond to configurations with Coulomb interaction, and those on the right to the case without such potential.

The change of the inner structure can be quantified through the use of the radial distribution function, g(r). An example of this effect is shown in Fig. 3, which shows this function calculated with the structures obtained with symmetric nuclear matter (x=0.5) at density ρ=0.072 fm^-3 and temperature T= 0.1 MeV and for varying strengths of the Coulomb potential, namely α=1 (full Coulomb), 0.8, 0.2, and 0 (without Coulomb). As it can be easily observed, the electron gas does not affect the nearest neighbor distances, although the magnitude of the function for the second nearest neighbors increases as Coulomb weakens.

Fig. 3.

Examples of the radial correlation function for varying strengths of the Coulomb potential: α=1 (full Coulomb), 0.8, 0.2, and 0 (without Coulomb).

The isotopic content on the cluster is also affected by the presence of electron gas. Fig. 4 shows the x content of the clusters formed in 200 configurations of asymmetric matter (x=0.3) at density ρ=0.015 fm^-3 and temperature T= 1.0 MeV. Clearly visible is the enhancement up to x 0.4 that would occur without Coulomb. Also noticeable is the increase of the maximum fragment sizes which grow up to A 800 without the electric interaction.

Fig. 4.

x content of the clusters formed in 200 configurations of asymmetric matter (x=0.3) at density ρ=0.015 fm^-3 and temperature T= 1.0 MeV. The top panel shows the case with Coulomb and the bottom one the case without such interaction. Notice that the abscissas have different scales.

As suspected in [25], in this study we found that the interactions between protons and electrons are responsible for the rearrangement of protons, although we found this to be of varying degree. Fig. 5 shows the RMS displacement of nucleons as a function of the evolution of the stationary ergodic process. Each curve is made of 200 points, each of which represents the average displacement (with respect to the original configuration) of the 2000 nucleons used in the x=0.5 cases shown. For each case tracked (ρ=0.015 and 0.72 fm^-3 at T=0.1 and 1.0 MeV), the displacements obtained with and without Coulomb are presented. Curiously, at low temperatures Coulomb seems to enhance the nucleon mobility much more than at high temperatures.

Fig. 5.

RMS displacement of nucleons during the evolution of the stationary ergodic process with Coulomb (continuous line) and without Coulomb (dashed line) as labeled; see text for details. Notice that at low temperatures the mobility with Coulomb is smaller than without Coulomb, while at higher temperatures the effect is the opposite.

A similar effect is observed in the persistence, which is the tendency of nucleons in a given fragment to remain in the same cluster. Fig. 6 shows the persistence 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 seen in the mobility, Coulomb enhances the transfer of nucleons thus decreasing the persistence as shown by the two curves.

Fig. 6.

Persistence with and without Coulomb as it evolves through 200 configurations of the stationary ergodic process.

The effect of electron gas on the shape of the structures formed can be quantified using the Euler Characteristic and Mean Curvature introduced before. Following the procedure explained in more detail elsewhere [21], digitized polyhedra were constructed for each of the 200 nuclear structures obtained with and without Coulomb at each value of {x,ρ,T} considered. The corresponding values of the mean curvature and Euler numbers were calculated and, although the configurations corresponding to the same set of conditions were not identical, their curvature and Euler numbers showed robust average values.

The actual values of the curvature in all the studied configurations varied approximately from -6,000 to 12,000 and those of the Euler numbers from -700 to 1,000. As these values depend on the characteristics of the polyhedra (i.e. number of sides, corners, etc.) they yield only relative information. Table 2 shows the curvature-Euler numbers as percentages of the maximum values obtained, and Fig. 7 shows their location on the curvature-Euler plane as well as their changes in position as the Coulomb strength is diminished. The standard deviations of the points are, roughly speaking, of the size of the points used in the plots.

Fig. 7.

Average values of the Curvature and Euler numbers of the structures listed in Table 2; circles correspond to structures with Coulomb and squares to structures without Coulomb, arrows indicate the average displacement of the structures as α goes from 1 to 0.

Classification Curvature - Euler

	ρ (fm^-3)	T (MeV)	α=1		α=0
	ρ (fm^-3)	T (MeV)	Curvature	Euler	Curvature	Euler
A	0.072	1.0	-14.5	8.9	-0.2	0.47
B	0.072	0.1	-3.18	-1.58	0.245	-0.19
C	0.015	0.1	14.76	0.89	4.28	0.35
D	0.015	1.0	18.43	-0.17	2.28	0.86
E	0.015	1.0	75.75	-61.2	91.76	-24.3
F	0.015	0.1	90.03	45.4	95.3	15.9
G	0.072	1.0	-24.74	44.7	-37.6	60.2
H	0.072	0.1	-37.38	79.8	-49.28	93.3

A. Analysis

The results obtained in this study allow us to understand the effect that electron gas has on the formation of the nuclear pasta. Perhaps the most interesting result of this investigation is the fact that the characteristic topological structures of the nuclear pasta –which are usually associated with the presence of competing interactions of the long range Coulomb potential and short range nuclear force– exist even when there is no Coulomb (i.e. when α 0). We believe this to be the result of the competition between attractive and repulsive nuclear forces at very low temperatures. In our case this is due to the attraction of the Vnp and the repulsion of the Vnn and Vpp which operates at medium ranges. Incidentally, we find the same effect also in other potentials [28] as was reported elsewhere [22].

Thus the role of the Coulomb forces introduced by the electron gas is to emphasize such competition between attractive-repulsive forces. The examples of structures formed with and without Coulomb (Fig. 1) suggest that the electric forces tend to distribute matter more and form less-compact objects. This behavior is in agreement with the reduction of the fragment size multiplicity as illustrated in Fig. 2 and 4.

At a microscopic scale, Fig. 3 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. 4. In agreement with this, the increment of nucleon mobility produced by Coulomb (cf. Fig. 5) gets reflected in a reduction of the persistence (cf. Fig. 6).

With respect to the modification of the shapes due to the electric interaction, we find that Coulomb tends to reduce both the curvature and Euler number of the structures. Indeed, Coulomb appears to "compactify" symmetric matter into spaghetti and lasagna type structures as shown clearly in Fig. 7 for symmetric matter (top panel). Of course, the effect becomes blurrier in systems with smaller charge densities, such as the asymmetric matter case of x=0.3 (bottom panel).

After having seen the results of the present study, it is now easier to gauge the effect of those quantum features not included in our classical model. CMD is at best an independent-particle approximation to a quantum system in which all higher order correlations are ignored. This excludes all effects ascribed to the closing of shells (magic numbers, shell polarization and deformations) and to any collective phenomena (pairing rotations and vibrations and their effect on fission decays, formation of Cooper pairs and their superfluidity effects). Although some of these effects are expected to be noticeable in the low density low temperature region of our study, it is easy to see that few of them –if any– have a bearing on the macroscopic structure of the clusters formed. Likewise, as all of the missing quantum effects are independent of the existence (or lack) of an embedding electron gas, it stands to reason that they would not have any effect on the results. In other words, had the comparison between pasta shapes with and without Coulomb interactions been performed with a model that included all quantum effects, the results would have been the same.

IV. CONCLUDING REMARKS

The effect of electron gas on symmetric and neutron-rich matter was studied at low densities and temperatures by varying the Coulomb interaction strength. 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 due, perhaps, to the interplay between the attractive and repulsive parts of the nuclear forces. In general, the presence of the electron gas tends to distribute matter more, form less-compact objects, increase nucleon mobility, reduce the persistence and the fragment size multiplicity, decrease the x content of the clusters, but does not change the inter-particle distances. All of this gets reflected on a modification of the shapes of the nuclear structures formed by a reduction of their curvature and Euler numbers. This effect is more pronounced in symmetric than asymmetric matter.

References

[1]

J Pochadzalla, T Mohlenkamp, T Rubehn, et al.

Probing the nuclear liquid-gas phase transition

. Phys Rev Lett, 1995, 75: 1040-1043. DOI: 10.1103/PhysRevLett.75.1040